Tight Cramér-Rao type bounds for multiparameter quantum metrology through conic programming

In the quest to unlock the maximum potential of quantum sensors, it is of paramount importance to have practical measurement strategies that can estimate incompatible parameters with best precisions possible. However, it is still not known how to ﬁnd practical measurements with optimal precisions, even for uncorrelated measurements over probe states. Here, we give a concrete way to ﬁnd uncorrelated measurement strategies with optimal precisions. We solve this fundamental problem by introducing a framework of conic programming that uniﬁes the theory of precision bounds for multiparameter estimates for uncorrelated and correlated measurement strategies under a common umbrella. Namely, we give precision bounds that arise from linear programs on various cones deﬁned on a tensor product space of matrices, including a particular cone of separable matrices. Subsequently, our theory allows us to develop an efﬁcient algorithm that calculates both upper and lower bounds for the ultimate precision bound for uncorrelated measurement strategies, where these bounds can be tight. In particular, the uncorrelated measurement strategy that arises from our theory saturates the upper bound to the ultimate precision bound. Also, we show numerically that there is a strict gap between the previous efﬁciently computable bounds and the ultimate precision bound.


Introduction
Quantum sensors, by employing quantum resources, promise to estimate physical parameters with unprecedented precision beyond what is possible using classical resources. Quantum metrology is a research field that studies quantum sensors. Quantum metrology schemes require the ability to both prepare parameter-dependent quantum probe states and perform quantum measurements on these states. Armed with the statistics of the measurement outcomes, one can thereafter estimate the underlying parameters. A central question in quantum metrology is to find measurement strategies with the ultimate precision for these multiparameter estimates. In the simplest scenario of estimating a single-parameter, the ultimate precision, given by the quantum (CR) Cramér-Rao bound [1,2,3,4,5], along with the corresponding optimal measurement strategy, are efficient to compute with knowledge of both the probe state and its dependence on the single parameter.
The theory of multiparameter quantum metrology 1 is considerably richer than the single parameter setting. For example, parameters can be fundamentally incompatible, as is often the case in quantum systems. If we want to design the best quantum sensors that can simultaneously estimate incompatible parameters, we must find optimal practical measurement strategies for multiparameter quantum metrology. However, even after decades of research, the question of how to determine these optimal measurement strategies remains unanswered and unknown.
In lieu of determining these optimal measurement strategies, the field has focused on determining bounds for the ultimate precision of quantum sensors that estimate incompatible parameters simultaneously. Since most optimizations for precision bounds [3,5,8,9,10,11,12] are not based directly on measurement strategies, even if we can calculate the best precision bounds from these optimizations, we still will not know what the optimal measurement strategies are.
Regarding the theory of precision bounds, the theory of single parameter estimation differs substantially from the multiparameter case. Namely, while the CR bound is tight for both correlated and uncorrelated measurement strategies in single-parameter estimation, this is not the case for multiple parameters. This is because the SLD Cramér Rao (SLD) bound does not give the tight bound in the multiple-parameter case nevertheless it gives the tight bound in the single parameter case. That is, for multiparameter quantum metrology, the Holevo-Nagaoka (HN) bound [3,4,5,8,9] is efficient to compute and always tight for correlated measurement strategies across multiple probe states. Although it is often called the Holevo Cramér Rao bound, it is called Holevo-Nagaoka bound in this paper and its reason is explained later. Such correlated measurement strategies however require a large quantum device across multiple probe states, this bound is not practical. To accomplish state estimation in a practical way, we need to design uncorrelated measurement strategies across multiple probe states, using only measurement devices that access individual probe states. (See Fig. 1.) The Nagaoka-Hayashi (NH) [4,10,11,12] bound is not only efficient to compute, but also addresses such uncorrelated measurement strategies. However, the NH bound might not be always tight for uncorrelated measurement strategies. Hence, an efficient way to determine the ultimate precision bound for uncorrelated strategies remains unknown.
The ultimate precision bound for uncorrelated strategies in the multiparameter setting, or simply the tight bound, was formulated more than two decades ago as the optimal value of an infinite-dimensional optimization program with linear objective and constraint functions on a topological vector space [13]. The advantage of this formulation of tight bound is that it is a direct optimization over measurement strategies in contrast to all other multiparameter precision bounds.
There however remain many outstanding questions pertaining to the tight bound.
(1) How to efficiently determine the tight bound?
(2) How to determine the optimal uncorrelated measurement strategy for multiparameter quantum metrology that saturates the tight bound?
(3) What is the relationship between the SLD bound, the HN bound, the NH bound, and the tight bound?
(a) (b) Figure 1: Measurement strategies in parameter estimation using quantum probe states can be either correlated (a), or uncorrelated (b) across identical copies of the probe states. In (a), the measurement device collectively measures the input states. That is, the measurement device needs to access a large quantum system. In (b) the measurement device individually measures the input states, but classical feedback is allowed to improve the measurement. This strategy is composed of measurement device to access a single system.

(4)
Is there a gap between the tight bound and the NH bound?
With regards to question (3), we unify the theory of the SLD, HN, NH and tight bounds under a common umbrella. Remarkably, these bounds can alternatively be formulated as conic programs with the same linear objective and constraint functions. The only difference between these programs is the different choices of their cones. Namely, the cone for the tight bound is a strict subset of the cone for the NH bound, and the cone for the NH bound is a strict subset of the cone for the HN and SLD bounds. This solves the open problem (3).
Regarding question (2), using our reformulation of the tight bound, we also construct an efficient algorithm that calculates the tight bound. From the dual program of our reformulation of the tight bound, we construct an associated semidefinite program (SDP), and show how to use its optimal solution to approximately solve the tight bound. Namely, we provide an efficient algorithm that calculates both upper and lower bounds to the tight bound, thereby solving open problem (4). Regarding question (4), using our algorithm, we numerically demonstrate that the tight bound can be strictly tighter than the NH bound.
We also solve questions (1) and (2), where we show concretely how we can efficiently compute optimal uncorrelated measurement strategies that saturate the tight bound along with the tight bound.
Our numerically tight estimates to the tight bounds are useful beyond multiparameter estimation theory. More abstractly, we develop an approach to optimize over the separable cone for bipartite systems. Hence we expect our theory can pave the way towards gaining new insights into optimizations over separable bipartite states, which could be useful for entanglement theory [14,15] and more general quantum resource theories [16,17,18]. Now we sketch the organization of our paper. We also give the structure of our paper visually in Fig. 2. In Table 2, we list the important notations we used in our paper along with their meanings. In Section 2 we review various CR-type bounds in multiparameter quantum metrology and define the tight bound. In Section 3, we answer question (3), where we unify various CR-type bounds via conic linear programming with various cones on a space that is the tensor product of real symmetric matrices and complex Hermitian matrices. Here, we reformulate the tight bound, (i.e.Theorem 2), the HN bound (Theorem 4) and the NH bound (Theorem 4). In Section 4, we develop our theory of how to calculate upper and lower bounds on the tight bound by applying an SDP with constraints labeled by unit vectors from a real vector space. We also derive an SDP with optimal value equal to the upper bound, and which is directly optimized over uncorrelated measurement strategies. From the solution for this SDP, we derive a concrete uncorrelated measurement strategy that has its precision given by our upper bound to the tight bound. Since our two-sided bounds to the bound can be tight, we therefore have derived a near-optimal measurement strategy We thereby can compute optimal uncorrelated measurement strategies in the asymptotic limit, and this answers question (2). In Section 5, we construct a strategic subset of the above unit vectors from design theory. In Section 6, based on the results of Section 5, we give the calculation complexity of the tight bound within an additive error of ϵ. This thereby answers question (1). In Section 7, we answer question (4), where we describe how we obtain our numerical lower bounds for the tight bound and illustrate its results. In Section 8, we explain how quantum multiparameter estimation theory is applicable in learning parameters of Hamiltonian models, and also in 3D-field sensing. In Section 9, we discuss our results and its implications in more detail.
2 Formulation and review of existing results

Various lower bounds for tight CR bound
We describe the formulation of quantum state estimation, and briefly review existing results on quantum parameter estimation. We refer readers to references [2,3,19,20,21,22] for more details.
A quantum system is represented by a Hilbert space H. In the following, when H is finite-dimensional, we can ignore the word "bounded" and "trace class", and can replace "self-adjoint operator" by Hermitian matrix. Let B sa (H) be the set of bounded self-adjoint operators on H, which is a real vector space. Let T sa (H) be the real vector space composed of set of trace class self-adjoint operators on H. A quantum state ρ is a positive semidefinite matrix on H with unit trace. The set of all quantum states on H is denoted by S(H) := {ρ | ρ ≥ 0, Trρ = 1} ⊂ T sa (H). Maximization of Tra + TrS subject to (23) for all x ∈ R d (a * * , S * * ) the optimal solution of D0 S(D0) = Tra * * + TrS * * the optimal value of D0,

Symbol
Nagaoka-Hayashi (NH) bound A measurement Π can be represented mathematically as a positive operator-valued measurement (POVM) which is a set of positive semidefinite matrices Π = {Π x } x∈X that satisfies a completeness condition. When the set of measurement outcomes is finite, X is a finite set and the completeness condition is x∈X Π x = I. When the set of measurement outcomes is a continuum, such as when X = R, the completeness condition is X Π x dx = I. For notational simplicity, for any POVM Π, we employ the notation for discrete-valued POVMs, and we always use x∈X Π x = I to represent the completeness condition.
When one performs a POVM Π on ρ, the Born rule gives the probability of getting an outcome x ∈ X . We are interested in a model of quantum parameter estimation given by a parametric family of quantum states on H: where Θ ⊂ R d denotes the set of parameters. Here, d denotes the number of parameters that we want to simultaneously estimate. To avoid mathematical subtleties, we impose regularity conditions; we require ρ θ to be differentiable sufficiently many times, assume ∂ρ θ /∂θ i to be linearly independent, and for the sake of clarity only consider full-rank states here 2 . Now given a measurement Π and an estimatorθ, we denoteΠ = (Π,θ) as an estimator. We define the mean-square error (MSE) matrix for the estimatorΠ as where E θ [f (X)|Π] denotes the expectation of a random variable f (X) with respect to the probability distribution p ρ θ (x|Π) = Trρ θ Π x obtained from the Born rule. In multiparameter quantum metrology, the objective is to find an optimal estimatorΠ = (Π,θ) that in some sense minimizes the MSE matrix.
Since the minimization of an MSE matrix is not properly defined, we seek precision bounds where we minimize Tr[GV θ [Π]], which is the weighted trace of the MSE matrix according to a given positive matrix G. Here, G is a weight matrix and quantifies the trade-off between estimating different vector components of the parameter θ. For instance, when G is the size d identity matrix I d , minimizing TrGV θ [Π] corresponds to minimizing the average variance of estimators.
A problem fundamental to quantum metrology is that of finding the ultimate precision bound under reasonable assumptions on the estimators we use. We say that an estimator Π is unbiased if for all θ = (θ 1 , . . . , θ d ) ∈ Θ, we have However, such an unbiased estimator invariably does not exist. Hence one often relaxes this unbiasedness condition to a locally unbiased condition in the neighborhood of a chosen point θ. To have the locally unbiased condition to hold, we require that for all parameter indices i, j ∈ {1, 2, . . . , d}, we have the equations Note that we can derive this condition by applying the Taylor expansion to the usual unbiasedness condition at a point θ to first order. Then, we introduce the fundamental precision limit by where the minimization is carried out for all possible estimators under the locally unbiasedness condition, which is indicated by l.u. at θ. In this paper, any lower bound for the weighted trace of the MSE matrix V θ [Π] is referred to as the CR type bound. When a CR type bound equals to the fundamental precision limit C θ [G] as in (6), it is called the tight CR bound in our discussion. That is, C θ [G] is called the tight CR bound. In the following, we discuss some CR type and tight CR bounds.
In fact, the set of MSE matrix V θ [Π] under the local unbiasedness condition is characterized as follows.
where J θ (Π) is the Fisher information matrix of the distribution family {P θ,Π } θ and P θ,Π (x) := Tr[ρ θ Π x ] [25, Exercise 6.44]. When multiple copies of the unknown state are prepared, only individual measurements for each copy are allowed, and the error is measured by weighted sum of mean square error with the weight matrix G, it is impossible to realize estimation precisions exceeding C θ [G]. Although the measurement to achieve this bound depends on the true parameter θ, when classical adaptive improvement for the choice of measurement is allowed, it is possible to achieve the bound C θ [G] [26,27,28,29]. Therefore, we can consider that the tight CR bound C θ [G] expresses the ultimate precision of the optimal estimator in the asymptotic limit of infinitely many probe states when only adaptive individual measurements for each copy are allowed. This setting corresponds to the strategy A2 in Section 3.2 of [22]. Furthermore, when n copies of the unknown state are given, even when any separable measurement over the n-fold system is allowed, it is impossible to overcome the precision error C θ [G] although there exists a separable measurement that requires quantum correlation over the n-fold system [25,Exercise 6.42].
We denote the minimum of (6) when our measurement is limited to measurement with discrete value and the number of elements in X is m by C [G, m].
In addition, since our interest is the minimization (6), without loss of generality, we can assume that θ is zero.
To get a CR type bound, we often focus on the SLD L i , which is defined as any Hermitian matrix that satisfies The SLD Fisher information matrix J = (J i,j ) is defined as Here, when ρ is strictly positive, the choice of Hermitian matrix L i is unique. Otherwise, it is not unique. However, the definition of the SLD Fisher information matrix J in (9) does not depend on the choice of Hermitian matrix L i under the condition (8). Under the locally unbiasedness condition at θ, we have SLD CR inequality [2] For the proof, see [2,3], [ (10) can be achieved by a local unbiased estimator constructed by their simultaneous spectral decomposition. In the choice of SLDs L i , extending the Hilbert space is allowed. However, when ρ is a strictly positive density matrix, it is sufficient to check for the commutativity of SLDs L i without extending the Hilbert space. In general, there is a possibility that the equality in (10) be achieved only with an extending Hilbert space. Taking a weighted trace in (10), we obtain the following bound.
• The SLD CR bound, which is the tight CR bound for any one-parameter model [2]: where J denotes the SLD Fisher information matrix about the model M.
To get a tighter bound than SLD bound, we focus on the vector of self-adjoint operators ⃗ Z = (Z 1 , . . . , Z d ) that satisfies the condition Then, we define the Hermitian matrix Z( ⃗ Z) whose (i, j) component is TrρZ i Z j . When an estimatorΠ = (Π,θ) satisfies the condition: we have [3, (6.7.73)] Using this relation, we have [4] TrGV θ [Π] ≥ TrGRe Z( ⃗ Z) + Tr|G where the operator |X| is defined as √ X † X. Therefore, we obtain the following bound; • The Holevo-Nagaoka (HN) bound [3,4] (Nagaoka [4] derived this bound by using several formulas obtained in Holevo [3], and called this bound the Holevo bound. Recently, the reference [30] called this bound the Holevo Cramer Rao bound, and cited Nagaoka [4] as its formulation. However, this bound does not appear in Holevo [3]. It should be called Holevo-Nagaoka CR bound. For its detailed history, see the latest review [31].) where the minimization takes the vector of self-adjoint operators ⃗ Z = (Z 1 , . . . , Z d ) to satisfy the condition (12). The HN bound C HN [G] improves the SLD CR bound C S [G], and gives the asymptotic limit of precision of the minimum estimation error when any quantum correlation is allowed in measurement apparatus and multiple copies of unknown states are prepared [33,5,34,35,29].
In addition, the HN bound C HN [G] satisfies the additivity condition, i.e., this value with the m-copy setting equals the value with the one-copy setting divided by m [5]. The reference [8,Eq. (11)] derived a calculation formula based on SDP for the bound C HN (G).
Using (17), we obtain the following bound.
• The Nagaoka bound [4,10], which is given only in the case with d = 2, and is tighter than the HN bound.
where the minimization takes the vector of self-adjoint operators ⃗ Z = (Z 1 , Z 2 ) under the condition (12).
In the qubit case, the Nagaoka bound C N θ [G] equals the tight CR bound C[G] for a twoparameter model (d = 2).
The reference [12, (22)] derived a calculation formula based on SDP for the bound C N H (G). Overall, the HN and NH bounds focused on the relation between the difficulty of joint measurement and the multiparameter estimation. This kind of relation was also discussed in the recent paper [36]. [37,13] To get another form of C[G], the papers [37,13] treated the minimization C[G] as a minimization with respect to a POVMΠ over R d . In this case, a POVMΠ is considered as an element of convex cone. Then, it focused on the following maximization problem.

Another expression of tight CR bound by
For a d × d real matrix a and a self-adjoint operator S on the Hilbert space H, we consider the condition: We consider the maximization: where the maximization takes the pair (a, S) to satisfy the condition (23). Here, D0 shows the maximization problem presented in (24), and the maximum value is denoted by S(D0). Then, the paper [13] showed the following proposition: Proposition 1 ([13, Theorem 6]).
In fact, the maximization (24) can be regarded as the dual problem when the minimization C[G] is considered as a conic linear programming with respect to a POVMΠ over R d . Hence, we call the maximization (24) the dual problem D0.
To see this, we focus on the relation where x = (x 1 , . . . , x d ) takes values in R d . Using the above relation, we have Hence, we can easily see the inequality ≥ in (25).
To show the opposite inequality, one needs to show the non-existence of the duality gap, i.e., the minimum of the primal problem equals the maximum of the dual problem in the framework of a conic linear programming, as explained in Appendix B. Although the non-existence of the duality gap holds for a conic linear program in a finite-dimensional vector space, the set of POVMsΠ over R d is a subset of an infinite-dimensional space, because the number of elements in R d is not finite, even when H is finite-dimensional. Hence, the paper [13] showed the non-existence of the duality gap in this problem setting by discussing a complicated issue related to topological vector space.
In fact, the formula (25) has various merits. Appendix A summarizes its two applications. For example, use of the relation (25) enables us to derive the minimum MSE under the locally unbiasedness condition under the one-parameter case. As another example, using the relation (25), the papers [37,13] derived the tight CR bound C[G] for a three-parameter model (d = 3) in the qubit system.

Conic programming with various cones on tensor product
When H is finite-dimensional, Proposition 1 guarantees that the tight CR bound C[G] is given as the maximum value S(D0) with the variables (a, S) of finite-dimension. However, even with a finite-dimensional space H, the calculation of the maximum value S(D0) is not so easy, To get a more computable form for C[G], we introduce the real vector space R spanned by {|0⟩, |1⟩, . . . , |d⟩}. Let M rs (R) be the set of real symmetric matrices on R. The weight matrix G can be considered as an element of M rs, This section aims to derive various CR bounds as conic linear programming problems on the vector spaces B := M rs (R) ⊗ B sa (H) and T := M rs (R) ⊗ T sa (H). That is, B is defined as an operator in T given by (57) real number depending on a, D j and ρ (see (80)) C 1 real number depending on D j and ρ (see (110)) δ(W R ) covering radius of W R κ related to the covering radius of W R Table 2: Notations from Section 3 onwards. For our conic programming framework, we consider primal conic programs P 1, P 2, P 3, P 4, P 5 and their dual conic programs D1, D2, D3, D4, D5. Primal programs are minimizations over cones, while dual programs are maximizations over the dual cones (also cones). There is no duality gap between the primal programs and the dual programs, so P i has the same optimal value as Di, so that S(P i) = S(Di). The optimizations for the primal programs are over variables of the form X, which are matrices in the tensor product space R C ⊗ H. The primal conic programs all have objective functions Tr(G ⊗ ρ)X and the dual conic programs all have objective functions Tra + TrS, where a is a size d real matrix, and S is a Hermitian matrix. Now S(P 1) and D(P 1) are equal to the tight bound C[G], and can be approximated using the semidefinite program (SDP) [P 1, W R ], where we can take W R to be a collection of points on a norm-1 hypersphere in a real vector space. Namely , denoting the optimal solutions as X * and (a * , S * ) respectively. Using (a * , S * ) we can calculate C 2 (a * ). Then we find that (in Theorem 13) T is similarly defined. Since any real symmetric matrix can be regarded as a Hermitian matrix on the complexified space R C of R, any element of B can be considered as a self-adjoint operator on R C ⊗ H. Then, we define several cones in the space B as follows. Let B sa,+ (H) be the set of bounded positive semi-definite operators on H. Then, we define the cone S SEP in B as S SEP := conv(M rs,+ (R) ⊗ B sa,+ (H)), where conv is the convex hull. Then, we define the cone S * SEP in B as the dual cone of S SEP . That is, In the context of entanglement theory, S * SEP is often called entanglement witnesses [38,39,40].
Also, we define the cone S P as Then, we have Hence, when H is finite-dimensional, we have Relaxing the condition of B we extend the space B as We define the set S(R C ⊗ H) P (S(R C ⊗ H) P P T ) of positive semi-definite (positive partial transpose) self-adjoint operators on R C ⊗ H. Since any component of an element X of S P is a self-adjoint operator, the element X satisfies the conditions of B ′′ . Also, the operator transposed on the system R C of X is positive. Hence, the element X satisfies the conditions of B ′′ and S(R C ⊗ H) P P T . Therefore, since the set S SEP is generated by separable pure states in R C ⊗ H, we find the relations To derive various CR bounds as conic linear programming problems on the vector spaces, we discuss the relation between our estimator and S SEP . For any estimator Although the references [41, (2011)] considered a matrix with the separable form from a POVM, our matrix X(Π,θ) with the separable form is different from their matrix with the separable form in the following point. Their matrix with the separable form is composed of the tensor product of the POVM element and the density matrix of their guess. We consider the tensor product of the POVM element and the one-dimensional operator given by the superposition of |0⟩ and the parameter of our estimate. The component |0⟩⟨0| enables us to check the condition for the POVM.
Since Π satisfies the condition of POVM, X(Π,θ) satisfies the following condition: The condition (5) for a locally unbiased estimator guarantees the following condition: Under the above condition, the objective function TrGV θ [Π] is rewritten with X(Π,θ) as In this notation, the components G 0,0 , G 0,j , and G j,0 are zero. Therefore, we consider the minimization: and we denote the above primal conic linear programming problem by P 1, and the minimum value is denoted by S(P 1). Although the above discussion shows the inequality C[G] ≥ S(P 1), we have the following theorem.

Theorem 2. We have
The above theorem shows the tight CR bound can be calculated by solving the primal conic linear programming P 1. Since the tight CR bound equals the tight bound, the primal conic linear programming P 1 characterizes the tight bound.
When H is finite-dimensional and dim H = n, the dimension of B is d(d+1) 2 · n 2 . Due to Carathéodory's theorem, any element of S SEP can be written as a convex sum of d(d+1)n 2 2 +1 extremal elements of S SEP . Then, we have the following lemma.
When we replace the cone S SEP by another cone in the primal conic linear program P 1, we obtain another primal conic linear program. The relation (34) lists various alternative cones. For instance, when the condition X ∈ S SEP is replaced by X ∈ S P , we denote the primal conic linear programming problem as P 2. In fact, since the cone S P is given as the set of positive semi-definite matrices in B, P 2 can be simply solved by semi-definite programming.
When the condition X ∈ S SEP is replaced by X ∈ S(R C ⊗ H) P P T ∩ B ′′ , (X ∈ S(R C ⊗ H) P ∩ B ′′ ), we denote this primal conic linear programming problem as P 3 (P 4). When we denote the minimum of the primal conic linear programming problem P l by S(P l) for l = 2, 3, 4. The relation (34) implies Then, the NH bound C N H [G] and the SLD bound C S [G] can be calculated by solving S(P 2) and S(P 4), respectively, as follows.

Theorem 4. We have
The above theorem shows that the NH and SLD bounds are the optimal values of the primal conic linear programs P 2 and P 4. Their difference lies in the choice of their cones. The cone of P 2 is composed of positive semi-definite matrices in B = M rs (R) ⊗ B sa (H), and the cone of P 4 is composed of all positive semi-definite matrices on R C ⊗ H that are also in the cone B ′′ .
Although the minimizations P 2, P 3, P 4 are defined in cones different from the cone S(R C ⊗ H) P , these minimizations P 2, P 3, P 4 can be considered as semi-definite programming (SDP) in the cone S(R C ⊗ H) P because their cones are given as linear constraints in the cone S(R C ⊗ H) P . Hence, the NH and SLD bounds C N H [G] and C S [G] can be solved by semi-definite programming while the reference [12, (22)] already gave the same SDP form as S(P 2) for C N H (G).
Since any element in B can be written in the above form due to the definition of B, it is sufficient to discuss S(P 2) to consider the matrix with the form X(X ′ , ⃗ Z). The objective functions of all of these conic programs is the same, and the only difference between them lies in the constraints that their operator variables must satisfy. The operator spaces for each of these conic programs are different; P1 is optimized over a bipartite separable space on the tensor product of R and H. All the other spaces use the complexification of R. The possible choices of X in the programs P 1, P 2, P 4 and P 5 are contained within one another, which we depict visually as nested rectangles. Because of this, we have S(P 1) ≥ S(P 2) ≥ S(P 5) ≥ S(P 4).
Any element of R C ⊗ H can be written as the form |(y, Considering the case with z = − d j=1 Z j y j , as pointed by [12, (22)], we find the equivalence between the following two conditions for X ′ and ⃗ Z.
(ii) The inequality ⟨y|X ′ − Π( ⃗ Z)|y⟩ ≥ 0 holds for any |y⟩ ∈ R ′ C ⊗ H. Since the condition (i) is the conditions for P 2 and the condition (ii) is the conditions for the NH bound C N H [G], the desired statement is obtained.
Next, we proceed to the proof of (45). We choose X ′ as an element of To get a relation with HN bound, C HN (G), we introduce a linear constraint to the operator X ∈ B ′′ as Tr X((|j⟩⟨i| − |i⟩⟨j|) ⊗ ρ) = 0 (48) for i, j = 1, 2, . . . , d. Using this linear constraint, we define the subspace B ′′ ρ of B ′′ as We consider the minimization: Since we have the relation we have the following relations Then, we have the following theorem.

Theorem 5.
We have The minimization P 5 can be considered as an SDP in the cone S(R C ⊗ H) P because the cone is given as linear constraints in the cone S(R C ⊗ H) P in the same way as the minimizations P 2, P 4. Therefore, the tight, NH, SLD, and HN bounds are summarized as Fig. 3. Although the reference [8,Eq. (11)] derived an SDP form of the HN bound, their SDP is different from P 5. They assumed the finite-dimensional system for H in the derivation [8,Eq. (11)]. In contrast, we do not require this assumption to derive the equation (53).
Proof. We show (53) by using the second expression of S(P 5) in (50). Since X(X ′ , ⃗ Z) ∈ S(R C ⊗ H) P ∩ B ′′ , the components of ⃗ Z are self-adjoint. First, we fix ⃗ Z. When the condition (48) holds, the condition X(X ′ , ⃗ Z) ∈ S(R C ⊗ H) P ∩ B ′′ is rewritten as follows.
|i⟩⟨j|⊗I satisfies the condition (48), when X ′ * satisfies the condition (48), when X ′ satisfies the condition (48). The pair of the condition (48) and the condition X(X ′ , ⃗ Z) ∈ S(R C ⊗ H) P ∩ B ′′ is rewritten as the pair of the condition (48) for X ′ * and the condition (54).
Since Tr (G ⊗ ρ)X ′ * = Tr G(Tr H ρX ′ * ) , our minimization with fixed ⃗ Z is rewritten as which equals Tr |G the minimum value of the objective function with fixed Since the condition for ⃗ Z is the same as the HN bound, we obtain (53).
Now, for a d × d real matrix a and a self-adjoint operator S on the Hilbert space H, let us define the operator as an element of T : The dual problem D1 of P 1 is given as the following maximization: Also, the dual problem D2 of P 2 is given as the following maximization: In the same way, the dual problems D3 and D4 of P 3 and P 4 are given as the following maximizations:
For its proof, see Appendix C. When H is a finite-dimensional space, the conic linear programing problems appearing in Theorem 6 are conic linear programing problems on finite-dimensional vector spaces.
Also, as shown in Appendix D, the following theorem holds.

Theorem 7.
We have The combination of Theorems 2, 6, and 7 implies Proposition 1. The combination of Theorems 2, 6, and 7 implies another proof of Proposition 1.

Upper bound
Although the minimization problem P 2, i.e., the NH bound C N H [G] can be solved by semidefinite programming, the cone in the primal conic linear programming P 1 is the separable cone S SEP , which is different from the set of positive semi-definite matrices. Hence, the primal conic linear programming P 1 is still not so easy even with a finite-dimensional space H. In fact, this type of problem appears in the membership problem of the separable cone S SEP , which appears in the decision problem of the existence of entanglement [14,15] and in generalized robustness of entanglement [17]. Also, the reference [42, Proposition 2] showed that the communication value can be calculated as a conic linear programming with the separable cone S SEP .
In this section, we consider an efficient algorithm to solve S(P 1). Although we present our algorithm for the conic linear programming S(P 1), this algorithm can be applied to a general conic linear programming with the separable cone S SEP including the problem presented in [ S(W R ), X s ∈ T sa,+ (H) with a block diagonal matrix m s=1 |s⟩⟨s| ⊗ X s ∈ T (W R ), we have another form of S[P 1, W R ]: where the condition C[P 1, We call the above problem [P 1, W R ]. This problem can be considered as an SDP with d 2 + n 2 constraints with block diagonal matrices, which is a typical case of sparse SDP. We define the number δ(W R ) as We consider a sequence W R,n such that δ(W R,n ) → 0. As shown in Theorem 13, we have Hence, we can use SDP for this problem. Many existing algorithms for SDP work well for sparse positive semi-definite matrices. The paper [43] studied the calculation complexity for generic primal-dual interior-point method for SDP. When we apply their analysis for the calculation complexity to the case with block diagonal matrices, as explained in Remark 8, we find that the calculation complexity of S[P 1, =O(m(n 6 + d 2 n 3 + d 4 n 2 )). (72)

Remark 8.
To get (72), we explain how to apply the analysis by the reference [43]. In the reference [43], m is the number of constraint, and n is the dimension of the vector space to define the positive semi-definite matrix. In their analysis for generic primal-dual interiorpoint method for SDP, the dominant part of the calculation complexity is given as the calculation of B ′ and r ′ in Eq. (8') of [43]. When the matrices are given as block diagonal matrices, the number of blocks is n 1 and size of each block is n 2 , the calculation complexity of B ′ is O(mn 3 2 n 1 + m 2 n 2 2 n 1 ) and the calculation complexity of r ′ is O(n 3 2 n 1 + mn 2 2 n 1 ). Therefore, the total calculation complexity is given as O(mn 3 2 n 1 + m 2 n 2 2 n 1 ). Applying this estimation to our setting, we obtain (72).

Estimator attaining upper bound
Here, we explain how to construct an estimator to attain the upper bound S[P 1, For the optimal solution X * 1 , . . . , X * m to the SDP [P 1, W R ], we can define the positive semi-definite matrix The first condition (68) implies Therefore, the estimator ({M s },θ) satisfies the locally unbiasedness condition. Its weighted mean square error is calculated as Therefore we have the following result.

Lower bound
Since the solution S[P 1, W R ] of the SDP gives an upper bound of S(P 1), we need to derive a lower bound for S(P 1). However, it is not easy to derive an lower bound of S(P 1) from the solution S[P 1, W R ]. In the following, we present a lower bound of S(P 1) under the algorithm presented in the above section. Although the algorithm and the upper bound presented in the above section work with a general conic linear programming with the separable cone S SEP , the following lower bound does not necessarily work with a general conic linear programming with the separable cone S SEP because the lower bound depends on the form of the conic linear programming S(P 1) with the separable cone S SEP .
To derive its lower bound, we denote the optimum a and S of the dual problem [D1, with optimal value S[D1, W R ] := i (a * ) i i + TrS * . Then, we define the real number for a d × d matrix a. Using this value, we define Then, we have the following lemma.

Lemma 10. A pair (a, S) satisfies the condition (23) if and only if the relation
Proof. In this proof, (c j ) j expresses the vector whose components are c 1 , . . . , c d . A pair (a, S) satisfies the condition (23) if and only if any y ∈ R d and |ϕ⟩ ∈ H satisfy The LHS is rewritten as That is, when |ϕ⟩ is fixed, the minimum of the LHS for y is realized when y = 1 2⟨ϕ|ρ|ϕ⟩ Therefore, the range of y ∈ R d in the condition (84) can be limited to the set y ∈ R d ∥y∥ ≤ Since we have using |ϕ ′ ⟩ := ρ 1/2 |ϕ⟩, we have Hence, the range of y ∈ R d in the condition (84) can be limited to the set y ∈ R d ∥y∥ ≤ C 2 (a) . Thus, we obtain the desired statement.
Then, for the lower bound of S(P 1), we have the following theorem.
Theorem 11. We have the following relation.
where n is the dimension of Hilbert space H. In particular, (94) holds with equality in the limit δ → 0.
When we calculate a lower bound of S(P 1) from S[P 1, W R ], we have knowledge for a * and X * in a real numerical calculation. Hence, it is possible to calculate κ by applying Theorem 11.
Proof. It is sufficient to show that the pair (a * , S * −κI) satisfies the condition (23) because in this case, the corresponding objective function of [D1, W R ] is Tra * + Tr(S * − κI). Due to Lemma 10, it is sufficient to show that the pair (a * , S * − κI) satisfies for y ∈ R d with ∥y∥ ≤ C 2 (a * ). For y ∈ R d with ∥y∥ ≤ C 2 (a * ) and v ∈ H with ∥v∥ = 1, we have which implies (88).

Theoretical Evaluation of κ
When we derive a lower bound of S(P 1) in our numerical calculation, we can calculate κ numerically, and can directly use Theorem 11. However, when we estimate the calculation complexity to get the tight CR bound within additive error ϵ, we need to evaluate κ theoretically. For this aim, we define Then, the following lemma gives an upper bound of κ.

Lemma 12.
The quantity δ(W R ) gives an upper bound of κ.
5 Construction of W R

Quantum t-design
To implement the proposed method, we need to construct the subset W R ⊂ R d+1 . For this aim, we like to pack m points uniformly on a d-sphere of radius 1, such that any point on the sphere is as close to some point as possible. As one possible choice, we randomly generate pure states |w 1 ⟩⟨w 1 |, . . . , |w l ⟩⟨w l | on the system H subject to the Haar measure. Another choice for W R is a quantum t-design.
A subset W H of normalized vectors on the finite-dimensional Hilbert space H is called a quantum t-design on H when where µ is the Haar measure on the set of pure states on R.
For a good choice of W R and W H , we can consider a quantum t-design because a quantum t-design has a symmetric property. For example, the paper [46] discusses approximate construction of quantum t-designs.

Spherical t-design
A quantum t design has not been sufficiently studied, but a spherical t-design has been well studied [47]. The references [48,49,50] study a spherical t-design. Fortunately, a quantum t-design can be constructed from a spherical t-design [44].
A subset V of S d−1 is called a spherical t-design when the relation holds any polynomial f with degree t, where µ S d−1 is the Haar measure on S d−1 . Given a spherical 2t-design V ⊂ S 2d−1 , we define a subset V H of pure states on H = C d as follows. For x ∈ S 2d−1 , we define the normalized vector w(x) ∈ H as w(x) j := x 2j−1 + x 2j i. As shown below, the set V H := {w(x)} x∈V is a quantum t-design. For any sequences e 1 , . . . , e t and f 1 , . . . , f t , we have Hence, we have which shows that the set V H := {w(x)} x∈V is a quantum t-design. In the same way, given a spherical t-design V ⊂ S d−1 , the set V is a quantum t-design on R = R d .

Construction reflecting structure of κ
Since the lower bound of S(D1) depends on the value κ, we may construct W by reflecting the structure of κ as follows. We choose a subset S ⊂ S d−1 , which can be a spherical t-design. Then, we define a subset S(ϕ) := {(cos ϕ, sin ϕy)} y∈S . We choose ϕ 0 > 0 as tan ϕ 0 = C 2 (a * ). Using k, we choose W R as ∪ k j=0 S( ϕ 0 j k ). We can modify this construction as follows. That is, to construct S( ϕ 0 j k ), depending on j, we choose a subset S j ⊂ S d−1 , which also may be a spherical t-design. We define a subset S( ϕ 0 j k ) := {(cos ϕ 0 j k , sin ϕ 0 j k y)} y∈S j . We choose ϕ 0 > 0 as tan ϕ 0 = C 2 (a * ). Using k, we choose W R as ∪ k j=0 S( ϕ 0 j k ). Indeed, when d = 2, the choice of S is very easy because we can choose S as {(cos 2πl N , sin 2πl N )} N l=1 . In this case, we find that δ(S) = ∥(1, 0) − (cos 2π/N, sin 2π/N )∥ = 2 − 2 cos 2π/N . Appendix E discusses the theoretical evaluation of κ in these choices while it is better to directly calculate κ from the obtained data a * and X * in a real numerical calculation.

Another construction
First, we construct the real discrete subset D n,d ⊂ S d , where S d is the d-dimensional sphere in R d+1 . We define D n,1 as (105) In general, we have This can be inductively shown as follows. Therefore, Thus, when δ(D n,d ) = δ, we have

Calculation complexity of tight CR bound within additive error ϵ
Now, we evaluate the calculation complexity of tight CR bound within additive error ϵ under the choice of W R given in Section 5.4. For further evaluation, we introduce the quantity: Then, we have the following lemma; Lemma 15.
Proof. We define A j = ρ −1/2 (D j ρ)ρ −1/2 . Then, since ∥a∥ 2 I ≥ a † a, we have Then, we have When the pair a * , S * is the optimal solution of S[D1, W R ] with a weight matrix W , we define ξ := n∥Π(a * , S * )∥(1 + ∥a * ∥ 2 C 2 1 ). Then, Theorem 13 and Lemma 15 guarantee that the error ϵ is upper bounded by ξδ(W R ). Using this fact, we have the following theorem Theorem 16. Suppose that we choose W R as D n,d with ϵξ −1 = δ(D n,d ). Then, the calculation complexity of the tight CR bound for probe states of size n within additive error ϵ is Proof. We consider how |W R | relates to the additive error of the CR bound. From Theorem 13, our lower bound to O(D0) depends on the covering radius of |W R . When the relation ϵ ≥ ξδ holds, Lemma 15 and (94) guarantee that the error of CR bound is upper bounded by ϵ.
Hence, substituting the above value to m in (72), the overall complexity of the algorithm is Since Theorem 16 evaluates the calculation complexity of the tight CR bound within additive error ϵ by using the value x, we need to evaluate the quantity ∥Π(a * * , S * * )∥, where (a * * , S * * ) denotes the optimal solution of D0. The following theorem gives its upper bound.

Theorem 17.
We denote the SLD Fisher information matrix by J. We have and Finally, we evaluate the calculation complexity of the tight CR bound within additive error ϵ by Theorems 16 and 17. Since a * * , S * * are close to a * (W R ), S * (W R ), ∥a * (W R )∥ and ∥Π(a * (W R ), S * (W R ))∥ have the same order as ∥a * * ∥ and ∥Π(a * * , S * * )∥, respectively. We assume that G is the identity matrix. Combining Theorems 16 and 17, we find that Hence, the calculation complexity is upper bounded as For simplicity, we consider the case when J is the identity matrix. The value TrJ −1 is d. The above value is simplified as

Numerical lower bound for S(D1)
Here we numerically find a lower bound for S(D1). Here n denotes the dimension of H, and d is the number of parameters that we estimate. Here, we consider the metrology problem with probe state where I n denotes an identity matrix of size n. We set the number of parameters to estimate to be and d = 2. Now let us define the algorithm MakeRandomDs(n) that generates size n random traceless matrices D 1 and D 2 . We can interpret this as corresponding to a quantum model with true parameters θ 1 and θ 2 , where if x is in the neighborhood of θ = (θ 1 , θ 2 ), we have ρ ((x 1 , x 2 )) ≈ I n /n + D 1 (x 1 − θ 1 ) + D 2 (x 2 − θ 2 ). Physically, this could correspond to quantum parameter estimation on a set of states that are close to maximally mixed, which describes the scenario where noise dominates the quantum system. Now we describe how we make the set W R . First we make the algorithm circlepoints(N ) which takes as input a positive integer N , and returns a set calS. We find that W R has 5750 points.
We implement our approximation of κ in Theorem 11 by the following method. 2. Set κ = −∞.

For x in myrange:
4. For y in myrange:
Otherwise, we set isgap = 0. Since (D 1 , D 2 ) is generated randomly by MakeRandomDs, isgap is a binary random variable. We define f n as the probability that isgap takes the value 1. Since f n is an unknown parameter, we estimate f n from our 50 numeric experiments that independently generates (D 1 , D 2 ), which determines isgap in the above method. Let g n be the number of times isgap = 1 out of the 50 experiments. We plot the values of g n along with their error bars in in Fig. 5.
For fixed n, we plot the largest values of S[D1, W R ]/S(P 2) among our 50 experiments in Fig. 6. Values of S[D1, W R ]/S(P 2) that are strictly larger than 1 illustrate a positive gap between S(P 1) and S(P 2). We also plot the corresponding values of S(P 4)/S(P 2) in Fig. 6.
As expected, when n = 3, S[D1, W R ] is never greater than the NH bound given by S(P 2). However, for n = 4, 5, 6, 7, 8, 9, in the majority of our numerical experiments, S[D1, W R ] exceeds the NH bound.

Learning parameters of Hamiltonian models
When one has a physical system on many qubits, the underlying Hamiltonian can always be written as H = k a k P k where P k are multiqubit Pauli operators. The Pauli operators When g n is close to 1, this indicates that it is very likely that all random derivatives of probe states have a C[G] > S(P 2). Hence, the numerical evidence show that there is very often a gap between the tight bound and the NH bound, demonstrating the benefit of using the tight bound over the NH bound. We calculate the 99% confidence interval of our estimator g n according to the Clopper-Pearson method (using the function binofit in MatLab.)  and S(P 5)]/S(P 2). In Fig. 5, we see that there is very often a gap between the tight bound and the NH bound. Here, this figure shows us what the magnitude of this gap is. Since the gap is quite small, we can see that there is numerical evidence that the NH bound is a good approximation to the tight bound. P k are known, but often, we do not know the precise values of a k . In a physical system with n qubits, there are up to 4 n coefficients for a k to determine. However, in realistic physical systems comprising of qubits, the predominant interactions are between pairs of qubits, and the interaction strength decreases as the separation of the qubits increases. When the interaction between the qubits is dominated by nearest-neighbor interactions, and the qubits are arranged in a 2D or 3D array, the vast majority of the coefficients a k have absolute values close to zero. In fact in this case, the number of unknown coefficients a k that we like to estimate will be polynomial in n. This can arise for instance in physical systems where only spatially adjacent qubits have a non-negligible interaction. For example, we could have Heisenberg model with qubits interacting according to low-dimensional graph [51], and the goal is to find the coefficients of the non-negligible interaction terms.
In general, given the Hamiltonian H = d k=1 a k P k , we like to estimate the parameters a 1 , . . . , a d using a initial n-qubit probe state ρ. After we initialize the physical system as the probe state ρ, we allow the physical system to naturally evolve with time according to the quantum channel N a 1 ,...,a d . In the noiseless setting, the quantum channel N a 1 ,...,a d describes evolution of quantum systems according to the Schrödinger equation, and is given by In more realistic scenarios however, the quantum channel N a 1 ,...,a d that models the evolution of the initial probe state, is not necessarily a unitary channel. In general, N a 1 ,...,a d can be described using Lindblad's formalism [52], where the channels are generated by exponentiating Liouville operators. Such models have been considered for instance in Ref [53]. We can see that in both settings with and without noise, the channel N a 1 ,...,a d is differentiable with respect to the parameters a 1 , . . . , a d . In the language of multiparameter quantum estimation, we can define the model which is differentiable set of quantum states. The task at hand is then to find the optimal quantum estimator that lets us learn the values of a 1 , . . . , a d while also minimizing the sum of the mean squared errors of a 1 , . . . , a d . Matsumoto [32] showed that when the model M is a set of pure states, then the corresponding tight bound is in fact equivalent to the HN bound. Since the noiseless setting corresponds to having N a 1 ,...,a d as a unitary channel, having the initial probe state to be a pure state would make the model M comprise of only pure states. In this case, we can simply evaluate S(P 4), the HN bound, and know that S(P 4) is equal to the tight bound. Now, if we consider the case where N a 1 ,...,a d is a non-unitary channel which is generated by exponentiating Liouville operators, even if the initial probe state ρ is pure, the model M need not be a set of only pure states. In this scenario, we have no guarantee that the tight bound is equal to the HN bound. To find the optimal quantum estimator, we can nonetheless use our method for finding the asymptotically for the tight bound.
As a concrete example, consider a two-qubit Hamiltonian H a,b = aX ⊗ X + bY ⊗ Y , where X and Y are the usual single qubit Pauli operators. Such a Hamiltonian H a,b models an XY Heisenberg interaction between a pair of qubits, and such interactions are ubiquitous in nature. Being able to estimate a and b is essential to characterize such Hamiltonians. We also consider having some dissipative jump operators L j that describe noise in the physical system. In this case, the Liouville operator is given by L a,b , where and γ j are real numbers. The quantum channel is and the corresponding model is M = {exp(L a,b (ρ)) : a, b ∈ R}. Our method can address such problems. As Section 7 illustrates, we are able to calculate tight two-sided bounds on the tight bound for such models with two-parameter estimation. We are also able to obtain corresponding near-optimal quantum estimators for the tight bound. Hence, for this example, we are able to determine the optimal uncorrelated measurement strategy to estimate a and b with the minimum mean-square error, which will be useful in characterizing interactions in spin-systems.

3D field sensing
The mathematical structure of the Hamiltonian learning problem in Section 8.1 also allows us to discuss the multiparameter quantum metrology of field sensing.
In the research area of 'field sensing' [54], we have a classical field that interacts with a physical system comprises of qubits, and the classical field interacts with each qubit in exactly the same way. The classical field is a 3D-vector (x, y, z) with three real components. The interaction Hamiltonian on n qubits is given by where X j , Y j and Z j denote multiqubit Pauli operators that apply the Pauli X, Y and Z respectively on the jth qubit, and apply the identity operator on all other qubits. For example when we have only two qubits, the Hamiltonian is Just as in Section 8.1, we consider the evolution operator to be generated by a Liouville operator L, to take into account the effects of noise in the quantum system. Namely, we consider L of the same form as (124), except that we replace the Hamiltonian H with what we consider in (127). In this case, the model is where ρ is a fixed two-qubit density matrix. We could for example consider ρ as a GHZ state. Then we can numerically calculate upper and lower bounds on the tight bound for such models, and furthermore also also obtain corresponding near-optimal quantum estimators for the tight bound. Hence, it is possible to use our algorithm for determining optimal estimators in multiparameter quantum metrology for the extremely practical problem of field sensing.

Discussion
To summarize our results, we have presented a unified viewpoint to understand existing bounds for the Cramer-Rao bound from the viewpoint of conic linear programming. However in general, existing bounds such as the NH bound and HCR bound, are not the tight Cramer-Rao bound. Therefore, we proposed an algorithm to calculate the tight bound by using conic linear programming. Therefore, the main contribution of our paper is to propose an algorithm to calculate the tight bound. Our method also enables us to find the optimal POVM for the parameter estimation.
In more detail, we have characterized three types of bounds, HN bound, NH bound, and the tight bound as the solution of conic linear program with different cones in the same linear space. Their size relationship is characterized as the inclusion relation among their corresponding cones. Since the conic linear programing corresponding to the HN and NH bounds are given as an SDP, they can be efficiently solved. However, since the cone corresponding to the tight bound is the separable cone S SEP , i.e., the cone composed of positive semi-definite operators with separable form over a bipartite system, it is hard to solve the corresponding conic linear programing. In the second part of this paper, we have tackled the problem, i.e., how to efficiently solve this conic linear programing. For this aim, we have proposed a method to convert the conic linear programing with the separable cone S SEP to an SDP with constraints labeled by unit vectors from a real vector space. We have derived upper and lower bounds on the tight bound by using the solution of the above converted SDP. From the optimal solution of the SDP, we give a concrete way to derive a corresponding uncorrelated measurement strategy which is asymptotically optimal. In addition, we have proposed the method to choose a strategic subset of the above unit vectors from design theory. Then, we have given the calculation complexity of the tight bound within an additive error of ϵ when the above method is applied. Then, we have numerically applied our method to our examples. In these examples, we have numerically shown that the tight bound is strictly larger than the NH bound, which shows the existence of the gap between the tight bound and the NH bound.
In fact, several problems of quantum information can be written with the separable cone S SEP . For example, the communication value can be calculated as a conic linear programming with the separable cone S SEP [42,Proposition 2]. Also, the existence of entanglement can be written as a problem by using the separable cone S SEP [14,15]. We can expect that our method to approximately solve conic linear programing with the separable cone S SEP can be extended to the above problems, and other similar problems [17,18]. Such extensions are interesting for future study.
Further, after completing this research, the reviewer informed us about the method discussed in the reference [55] to calculate the conic linear programing with separable cone. The comparison between our method and this approach is an interesting topic, but beyond the scope of this paper. Also, as another application of our method includes the simultaneous estimation of phase and dephasing which was discussed in [56]. There is a possibility that our presented approach provides a new insight into this problem. This kind of research is another interesting future study.
Another open problem is whether our conic programming framework can also be used to unify multiparameter quantum metrology problems where instead of minimizing the mean square error, one minimizes a different regret function. For instance in the reference [36], minimization of an information regret function was considered. It will be interesting to see if our framework also applies in that setting.
We set a j i = a √ g j δ j i and S = − a 2 4 (I − ρ). The condition (23) can simplified as Hence, the objective value of the dual problem is Hence, we find that ( j √ g j ) 2 is a lower bound for CR bound.
Moreover, we expect that ( j √ g j ) 2 is the solution with a = 2 j √ g j .
This expectation can be checked by constructing a locally unbiasd estimator to realize the equality.
We denote the spectral decomposition of L j as k s j,k E j,k . Define the vector x(j, k) ∈ R d as follows. It has only non-zero element in the j-th entry. Its j-th entry is s j,k Then, we define the locally unbiasd estimator as follows. It takes the outcomes in the set Here, we multiply the probabilistic factor √ g j j ′ √ g j ′ in the POVM element. Hence, to satisfy the locally unbiasedness condition, we need to multiply its inverse j ′ √ g j ′ √ g j in the measurement outcome. Hence, we can check that this POVM satisfies the condition for the locally unbiasedness condition.
In fact, when x = x(j, k), Hence, we find that This equation shows the lower bound ( j √ g j ) 2 is attainable.
Proof of (131). We choose σ i as Without loss of generality, we can assume that ρ = 1 2 (I + α(µ)σ 3 ) with 0 ≤ a ≤ 1. Notice that the SLD Fisher information matrix is identical. The condition (131) is covariant. That is, if the condition (131) under a coordinate, it holds even with another coordinate converted from orthogonal transformation. Hence, it is sufficient to check the condition (131) only under a specific coordinate. Now, we consider the case that In this case, the SLD Fisher information matrix is identical, and the condition (131) holds.

B Linear programming with general cone
We consider conic linear programming (EP) with a general cone. Let X and Z be topological real vector spaces. Let P ⊂ X be a cone. Let A be a linear function from X to Z. Given c ∈ X * and b ∈ Z, we consider the following minimization.
As the duality problem, we consider the following maximization.
When the pair of x ∈ P and w ∈ Z * satisfy the constraints, we have c( holds. The difference EP − EP * is called the duality gap, and its existence/non-existence is not trivial in general. To discuss this problem, we define two sets: Since F ∩ E ⊂ F ∩ E and F ∩ E is a closed set, we have the inclusion relation The opposite inclusion relation is not trivial, and is related to the duality gap as follows.

Theorem 18 ([58, Proposition 4.2]). When
we have When X and Z are finite-dimensional, the relation (142) holds. Hence, there is no duality gap.

C Proof of Theorem 6
This appendix aims to show the non-existence of the duality gap in the respective problems P l for l = 1, 2, 3, 4, 5 by using Theorem 18 based on the notations given in Appendix B. Since the condition (143) holds in the finite-dimensional case, Theorem 6 immediately follows from Theorem 18. This point is in contrast with the proof of Proposition 1 by the paper [13] because the paper [13] had to manage a difficulty in the possibility of the existence of the duality gap even with a finite-dimensional space H. The difficulty in [13] was caused by handling the space of POVMs with Ω = R d as an infinite-dimensional linear vector space. Hence, in our proof, we encounter the difficulty related to the potential existence of the duality gap only with an infinite-dimensional space H.
Step 1: Preparation and structure of this proof.
We define X as the set of bounded operators on R d+1 ⊗ H. We define Z as the tensor product of R d 2 and the set of bounded operators on H. Then, the linear map A is defined as the left sides of (36) and (37). The element b ∈ Z is defined as the right sides of (36) and (37). The element c ∈ X * is defined as (38). Then, for l = 1, we choose P as S SEP . For l = 2, we choose P as S P . For l = 3, we choose P as S(R C ⊗ H) P P T ∩ B ′′ . For l = 4, we choose P as S(R C ⊗ H) P ∩ B ′′ . For l = 5, we choose P as S(R C ⊗ H) P ∩ B ′′ ρ . We consider the operator norm in the space F when H is infinite-dimensional.
We denote the SLD Fisher information matrix by J i,j , and its inverse matrix by J i,j . We define L j := i J i,j L i . Hence, we have Due to Theorem 18, it is sufficient to show (143). In the following, we show (143). Given an element s ∈ R with the condition (s, b) ∈ F ∩ E, we choose a sequence {X n } in X such that {(c(X n ), AX n )} converges to (s, b). It is sufficient to show that there exists a sequence {Ẋ n } in P such that AẊ n = b and c(Ẋ n ) converges to s, which shows that F ∩ E ⊂ F ∩ E.
Step 2: Definition of a new elementẊ n ∈ P.
We define the d × d matrix a (n) with component a We define the vectorȧ (n),i ∈ R d asȧ by u x u E u . Then, we define the matrix X n,i : Then, we have We defineẊ where Since X n,i belongs to P,X n,i belongs to P for any l = 1, 2, 3, 4. Hence, T nXn T n belongs to P for any l = 1, 2, 3, 4. Thus,Ẋ n belongs to P for any l = 1, 2, 3, 4.
Step 3: Asymptotic behavior of new elementẊ n ∈ P.

D Proof for Theorem 7
The reference [13, Theorem 6] spends long pages to show Proposition 1 because the conic linear programming P 0 is given as a conic linear programming in an infinite-dimensional space. However, due to Theorems 2 and 6, we can show Proposition 1 by showing the relation S(D0) = S(D1).

Theorem 19. We have
Proof. The condition for D0 can be rewritten as Then, the condition (169) is rewritten as 0 ≤Tr |v(x)⟩⟨v(x)| ⊗ |y⟩⟨y| for y ∈ H and x ∈ R d . This condition is rewritten as 0 ≤Tr |v⟩⟨v| ⊗ |y⟩⟨y| for y ∈ H and v ∈ R to satisfy ⟨v|0⟩ ̸ = 0. When (170) holds for y ∈ H and v ∈ R to satisfy ⟨v|0⟩ ̸ = 0, it holds for y ∈ H and v ∈ R because of the following reason. We choose v ∈ R such that ⟨v|0⟩ = 0. Then, we choose a sequence v n such that v n → v and ⟨v n |0⟩ ̸ = 0. Since (170) holds with v n , by taking the limit, (170) holds even with v.
Therefore, the condition for D0 can be rewritten as the condition that the inequality (170) holds with any y ∈ H and any v ∈ R. This condition is equivalent to (58). Hence, we conclude that the two dual problems D0 and D1 are equivalent. E Evaluation of κ based on W R constructed in Section 5.3 Since the lower bound of S(D1) depends on the value κ, it is better to construct W by reflecting the structure of κ. Hence, we choose W R as follows. We choose a subset S ⊂ S d−1 . We define a subset S(ϕ) := {(cos ϕ, sin ϕy)} y∈S . We choose ϕ 0 > 0 as tan ϕ 0 = C 2 (a * ). Using k, we choose W R as ∪ k j=0 S( ϕ 0 j k ). Indeed, when d = 2, the choice of S is very easy because we can choose S as {(cos 2πl N , sin 2πl N )} N l=1 . In this case, we find that δ(S) = ∥(1, 0) − (cos 2π/N, sin 2π/N )∥ = 2 − 2 cos 2π/N .
By using a vector y = t −1 x ∈ R d , (191) is rewritten as That is, it is rewritten as which implies (189).
In the following, we assume that G is the identity matrix.
Proof. The condition (23) for (a, S(a)) is written as follows.
We take the diagonalization as 1 2 d j=1 a j i L j = k c k |ϕ k ⟩⟨ϕ k |. We choose x i ′ = 0 for i ′ ̸ = i and x i = c k in (195). Then, we have Hence, we have where (a) follows from the fact that c k is the eigenvalue of 1 2 d j=1 a j i L j with the vector ϕ k . Also, (b) follows from (196).
Taking the sum with respect to k, we have − Tr Lemma 25.
Proof. We have where (a), (b), and (c) follow from the inequality ( d j=1 where (a) follows from Lemma 23.
The equality holds when a is a constant times of J −1 . Therefore, we have the relation max a d 4Tr[a T Ja] The combination of (209) and (212) yields (203). Also, (204) is shown as where (b), (c), and (d) follows from the definition of t a , (210), and (212), respectively. In addition, (a) follows from (208) and the relation Tra * = 1.
Proof of (116). When G is the identity matrix, (116) follows from Lemmas 25 and 26.
F.2 Proof of (117) Lemma 27. When G is the identity matrix, a * * is a positive semi-definite matrix.
Proof. Assume that (a, S) satisfies the condition (23). We choose the polar decomposition of a as a = |a|o with an orthogonal matrix o. Then, we have for y ∈ R d . We choose y ′ = oy. Then, we have Hence, (|a|, S) also satisfies the condition (23) and Tr|a| ≥ Tra. In fact, the equality of the above inequality holds if and only if a is a positive semi-definite matrix. Therefore, we can conclude that the optimizer a * * is a positive semi-definite matrix.