Optimal probes and error-correction schemes in multi-parameter quantum metrology

We derive a necessary and suﬃcient condition for the possibility of achieving the Heisenberg scaling in general adaptive multi-parameter estimation schemes in presence of Markovian noise. In situations where the Heisenberg scaling is achievable, we provide a semideﬁnite program to identify the optimal quantum error correcting (QEC) protocol that yields the best estimation precision. We overcome the technical challenges associated with potential incompatibility of the measurement optimally extracting information on diﬀerent parameters by utilizing the Holevo Cramér-Rao (HCR) bound for pure states. We provide examples of signiﬁcant advantages oﬀered by our joint-QEC protocols, that sense all the parameters utilizing a single error-corrected subspace, over separate-QEC protocols where each parameter is eﬀectively sensed in a separate sub-space.

Recently, a general theory providing a necessary and sufficient condition, the HNLS condition (an acronym for "Hamiltonian-Not-in-Lindblad-Span"), for achieving the HS in a finite-dimensional system in the most general adaptive quantum metrological protocols under Markovian noise, has been developed [35,36]. The theory allows for a quick identification of the most promising quantum metrological models and provides a clear recipe for designing the optimal adaptive schemes based on appropriately tailored QEC protocols. However, HNLS is restricted to the single-parameter estimation case, while a lot of relevant metrological problems, like vector field sensing (e.g. magnetic field) [43], imaging [44], multiple-arm interferometry [45,46] or waveform estimation [47,48] are inherently multi-parameter estimation problems. Multi-parameter estimation problems drew a lot of attention in recent years [49][50][51][52][53][54][55][56], yet the fundamental questions regarding the achievability of the HS as well as the form of the optimal metrological protocols in multiple-parameter estimation in presence of noise have not been answered so far. The aim of this paper is to fill this gap. The methods developed for the single parameter estimation case, in particular the semidefinite program (SDP) that allows to identify the optimal QEC protocol [36], are not applicable in the multi-parameter regime. The reasons are threefold.
First, the widely used quantum Cramér-Rao (CR) bound for multiple parameters is not in general saturable, due to the incompatibility of the optimal measurements for different parameters [49,51,52]. Therefore, unlike in the single-parameter case, the quantum Fisher information (QFI) does not provide the full insight into the problem [51,[57][58][59][60]. On the other hand, stronger bounds, such as the HCR bound [59][60][61][62][63][64], have no closed-form expressions (except for specific cases [65]). Moreover, the HCR bound, although solvable via an SDP [66], does not shed light on the corresponding optimal measurements saturating it. In general, the HCR bound is only saturable when collective measurements on all copies of quantum states are allowed, which is a demanding condition in practice [63,64]. As a result, the optimal measurements on the output quantum states are hard to identify.
Second, there is no general recipe to find the optimal input states in multi-parameter estimation even in the noiseless case [52,53], unlike in the single-parameter estimation case where the optimal input state is simply the equally weighted superposition between the eigenstates corresponding to the maximum and minimum eigenvalues of the Hamiltonian.
Finally, in the single-parameter case [36], all valid two-dimensional QEC codes were mapped into a traceless Hermitian matrix representing the difference between logical zero and one codes. In the multi-parameter case, however, it is not clear whether valid QEC codes could be mapped into a convex set when the code dimension is large, which is inevitable in multi-parameter estimation.
In this paper, we generalize the HNLS condition to multi-parameter scenarios, and provide an SDP to find the best possible QEC protocol in the situations when the HNLS condition is satisfied (including all noiseless cases). The solution yields an explicit form of the optimal input state, QEC codes and measurements. No collective measurements are required on the output states. Our protocol goes beyond the typically used QFI-based formalism and overcomes all the challenges related with the multi-parameter aspect of the problem mentioned above. Our work reveals the advantage of QEC protocols in multi-parameter estimation and we expect that the SDP formulation of our problem will also be an inspiration for other research areas in quantum error correction and quantum metrology.

Formulation of the model
We assume the dynamics of a d-dimensional probe system H S is given by a general quantum master equation [67][68][69]: where the parameters to be estimated ω = [ω 1 , . . . , ω P ] enter linearly into the Hamiltonian of the evolution via Hermitian generators G = [G 1 , . . . , G P ] T (where T denotes transpose) so that H = ω · G ≡ P k=1 ω k G k , and L k are operators representing a general Markovian noise. Similar to the previous investigations [35,36] we consider the most general adaptive scheme (see Fig. 1) [27] with an unlimited number of ancillae (denoted jointly as H A ), instantaneous perfect intermediate unitary operations U i and a general POVM on the final state ρ ω . E ω t represents the probe system dynamics integrated over time t, whereas the total probe interrogation time is T . Such schemes are the most general schemes of probing quantum dynamics, assuming the total interrogation time is T , and encompass in particular all QEC procedures.
In single-parameter estimation the optimal protocol is the one that yields the minimum estimation variance. In multi-parameter case the estimator covariance matrix is the key object capturing estimation precision, defined as [58,59]: for i, j = 1, . . . , P , where M ≥ 0, M = 1, are measurement operators ("≥ 0" for matrices means positive semidefinite) andω( ) is an estimator function mapping a measurement result to the parameter space. Diagonal entries of Σ represent variances of estimators of respective parameters while offdiagonal terms represent covariance between different parameters. As a figure of merit one may simply choose Tr(Σ) which will be the sum of all individual parameter variance, or more generally Tr(W Σ), where W is a real positive cost matrix that determines the weight we associate with each parameter in the effective scalar cost function Note that we require strict positivity of W which is equivalent to saying that it is an estimation problem of all P parameters, and not a problem where effectively only a smaller number of parameters are relevant. We assume the measurement-estimation strategy to be locally unbiased at some fixed parameter point ω, i.e.

The necessary and sufficient condition for the HS
We say that the HS in a multi-parameter estimation problem is achieved when there exists an adaptive protocol such that for every W > 0, ∆ 2 Wω ∝ 1/T 2 in the limit T → ∞. This is equivalent to a requirement that all parameters (and any combination of parameters) are estimated with precision that scales like the HS. The following theorem generalize the HNLS condition [35,36] to multi-parameter scenarios.

Theorem 1 (Multi-parameter HNLS). The HS can be achieved in a multi-parameter estimation problem if and only if
in the Hilbert space of Hermitian matrices under the standard Hilbert-Schmidt scalar product, whereas the superscripts H , AH denote the Hermitian and anti-Hermitian part of an operator respectively.
Proof. Let us start with a brief review of the HNLS condition in the single-parameter case, where H = ωG involves only a single generator G. As shown in [35,36], the necessary and sufficient condition to achieve the HS is that G / ∈ S, or in other words that G ⊥ = 0. In particular, following [36] (see the section named "QEC code for HL scaling when HNLS holds"), an explicit construction of the optimal QEC code was provided, where the code space H C ⊆ H S ⊗ H A is defined on the Hilbert space of the probe system H S extended by an ancillary space H A ∼ = H S . The code space satisfies the QEC condition [36,72]: where the operator S acting on H S was tensored with identity on H A and Π H C denotes the projection onto H C . Metrological sensitivity is guaranteed by the fact that G acts non-trivially on H C : As a result we obtain a noiseless unitary evolution generated by G C leading to the HS in the estimation precision of ω.
(Necessity) Suppose (G i ) ⊥ 's are linearly dependent. Then there exists a linear (invertible) transformation on the parameter space A ∈ R P ×P : ω = ωA −1 , (where we also modify accordingly the generators G = AG and the cost matrix W = AW A T , so that H and ∆ 2 Wω remain unchanged), such that (G i ) ⊥ = 0 for some i. Then, from the single-parameter theorem, ω i cannot be estimated with precision better than ∆ 2ω i ∼ 1/T which contradicts the HS requirements. (Sufficiency) Suppose (G i ) ⊥ 's are linearly independent. We assume the ancillary space to be a direct sum of P subspaces H Ai so that the whole Hilbert space is H S ⊗ (H A1 ⊕ · · · ⊕ H A P ) (see Fig. 2a). We may construct separate code spaces for each parameter using orthogonal ancillary subspace H Ci ⊆ H S ⊗ H Ai so that the QEC conditions Eq. (6) are satisfied within each code space H Ci separately. While constructing the code space for the i-th parameter, we include all the remaining generators G j (j = i) in the Lindblad span, so effectively treating them as noise i.e.
As a result thanks to the QEC and hence within a given subspace only one parameter is being sensed via the effective generator G Ci generators act trivially. If |ψ i ∈ H Ci is the optimal state for measuring ω i , the state to be used in order to obtain HS for all parameters which is not affected by noise reads ρ in = 1 H Ai -there is no measurement incompatibility issue because different parameters are encoded on orthogonal subspaces.
Similar to the single-parameter case, it must be admitted that in realistic situations with generic noise, HNLS is often violated [24,35]. Therefore, a more pragmatic approach is required taking into account the fact that in a real experiment the total time of evolution T is always finite. Let us consider a situation where H ∈ S, but where some noise components are weak [36]. Specifically, we divide Lindblad operators in Eq. (1) into two sets-strong noise generators {L k } and weak noise generators {J m } satisfying := m J † m J m where · denotes operator norm. If the HNLS condition is satisfied for the strong noise part, we could choose the code space H C which allows to completely erase the strong noise {L k } and the resulting effective noise rate would be upper bounded by [36]. This means that the distance between state of the error-corrected probe and the state evolving under ideal noiseless evolution will be of the order Θ( T ). Therefore, for sufficiently short evolution, T = o(1/ ), the precision of estimation will still scale quadratically with the total time ∆ 2 Wω ∝ 1 T 2 whereas for larger T , it will gradually approach the standard 1/T scaling.

Optimal probes, error correction schemes and measurements
In Sec. 3, we provided a QEC code where each parameter is sensed separately in different errorcorrected subspaces (see Fig. 2a). Such protocols will be referred as separate-QEC protocols (SEP-QEC). In contrast to this construction, we will now consider QEC strategies which allow for simultaneous estimation of all the parameters in a single coherent protocol by utilizing states within a single protected code space, which we will call the joint-parameter QEC scheme (JNT-QEC). In this section we provide a general method to identify the optimal JNT-QEC, while its potential advantages over the optimal SEP-QEC will be discussed in Sec. 5, as well as in Sec. 6 when studying concrete estimation models. From now on, we assume the multi-parameter HNLS condition is satisfied. Without loss of generality, we assume the generators {G i } P i=1 ⊂ S ⊥ are orthonormal, since the components in S do not contribute and there is always a linear transformation A on parameters leading to orthonormality. The following theorem provides a recipe to find the optimal JNT-QEC.
Theorem 2 (Optimal JNT-QEC). Given a cost matrix W . If the multi-parameter HNLS condition is satisfied with generators {G i } P i=1 ⊂ S ⊥ , the minimum cost ∆ 2 Wω that can be achieved in a JNT-

an orthonormal basis of Hermitian operators acting on H
i , B i are Hermitian (P + 1) × (P + 1) matrices (with matrix indices taking values from 0 to P ), and ν i ∈ R. Γ and K are real P × P matrices (with matrix indices from 1 to P ).
As a semidefinite program (SDP) it could be easily solved numerically, for example, using the Matlab-based package CVX [73]. Before giving a formal proof, let us briefly review some existing bounds in multi-parameter metrology and discuss their saturability.

General bounds in multi-parameter metrology
In this section, we discuss bounds for general parameter estimation problems on fixed quantum states. We use H to denote a general finite-dimensional Hilbert space.
Most commonly, quantum multi-parameter estimation problems are analyzed utilizing the standard quantum CR bound [57][58][59]: where F is a P ×P QFI matrix and Λ i (symmetric logarithmic derivatives) satisfy This bound is not saturable in general, due to potential non-compatibility of the optimal measurements, unless Im[Tr(ρ ω Λ i Λ j )] = 0 [51]. Moreover, a direct minimization of the CR bound over all JNT-QEC, with the saturability constraint imposed, does not necessarily guarantee the identification of the optimal protocol-there is a possibility that the optimal protocol does not meet the saturability constraint for the CR bound.
Therefore, in order to identify truly optimal protocols we need to resort to a stronger HCR bound [59,61,62]: where L(•) denotes the set of all linear operators acting on •, Re and Im denote the real and imaginary part of a matrix (not to be confused with the Hermitian and anti-Hermitian part of a matrix), and Tr[abs(·)] is the sum of absolute values of the eigenvalues of a matrix. When the second term is dropped the HCR bound reduces to the standard CR bound [51,59]. Unlike the CR bound this bound is saturable in general using collective measurements on many copies [63,64]. On the other hand, the HCR is defined via an optimization problem, making it usually more difficult to deal with than the standard quantum CR bound with a closed-form expression.
In the case of pure states ρ ω = |ψ ω ψ ω |, however, we note that the HCR is exactly equivalent to [49]: which we will call the Matsumoto bound. It was also shown in [49] that the bound is always saturable using individual measurements (which relaxes the requirement of collective measurements in the pure state case). However, there are no known efficient numerical algorithms to find the solutions {|x i } of the Mastumoto bound and it is not clear yet whether or not the Mastumoto bound will be useful in finding the optimal QEC protocols in our situation. In the following, we will start the proof of Theorem 2 by first sketching the proof of the Matsumoto bound and then reformulating it in such a way that the final optimization problem becomes an SDP, which will eventually be incorporated into the QEC protocol optimization procedure.

Proof and reformulation of the Matsumoto bound (Eq. (11))
According to the Naimark's theorem [59] One may see that, thanks to the projective nature of measurements {E }, scalar products of vectors |x i yield the covariance matrix of the estimator: Now, instead of minimizing over the measurement {M } on H, we can perform the minimization directly over the vectors |x i ∈ H M , imposing the following constraints: These constraints correspond respectively to the projective nature of the measurement {E } and the unbiasedness conditions as given in Eq. (4). At this point one may wonder how big the space H M should be (as for a general measurement it might be arbitrary large). However, we can always map span{|ψ ω , {∂ i |ψ ω , |x i } P i=1 } ⊆ H M isometrically to a (2P +1)-dimensional space. Therefore when looking for the bound, under the constraint Eq. (14), it is enough to perform the minimization over |x i ∈ span{|ψ ω , ∂ 1 |ψ ω , ..., ∂ P |ψ ω } ⊕ C P , which results in equation Eq. (11).
Finally, we show that indeed for any set of |x i satisfying Eq. (14) there exists a proper projective measurement on H ⊕ C P and a locally unbiased estimator satisfying Eq. (12), and consequently there exists a corresponding general measurement on H. To see this, notice that since ∀ i ψ ω |x i = 0 and ∀ i,j x i |x j ∈ R one may choose an orthonormal basis {|b i } of span{|ψ ω , |x 1 , . . . , |x P } satisfying: ∀ i ψ ω |b i ∈ R\{0} and ∀ i,j x i |b j ∈ R. Then one can define a projective measurement: with the corresponding estimator: which is locally unbiased and satisfies This proves Eq. (11). Specifically, if dim(H) ≥ 2P + 1 we may choose span{|ψ ω , ∂ 1 |ψ ω , ..., ∂ P |ψ ω } ⊕ C P as a subspace of H and optimize over |x i ∈ H. In this case we may also reformulate the Matsumoto bound in a slightly different form. First, note that any vectors {|x i } satisfying Eq. (14) need to be linearly independent. Let {|c i } P i=1 be an orthonormal basis of span{|x 1 , ..., |x P }, satisfying ∀ i,j Im x i |c j = 0 (such a set may be generated using the Gram-Schmidt orthonormalization procedure). The locally unbiased conditions may now be rewritten as: which (after introducing matrices This formulation will be more convenient to use when we will formulate the QEC protocol optimization problem as an SDP.

Optimizing the error-correction codes
Now we apply the reformulated Matsumoto bound to our task of identification of the optimal JNT-QEC. Consider a given input state |ψ in . Let H C be any code subspace of H S ⊗H A containing |ψ in and satisfying the QEC conditions Eq. (6)-in order to be in accordance with the reformulated Matsumoto bound, this space may be required to be at least 2P + 1 dimensional, but as we show in the following it will always be possible to reduce its dimensionality to P + 1 effectively. Using QEC, our goal is to preserve an effective unitary evolution in the encoded space and coherently acquire the sensing signal. Therefore, we are effectively dealing with pure state |ψ ω , which allows us to utilize Eq. (19) as a formula for the minimal cost of sensing multiple parameters. The effective evolution after implementing QEC is given by We focus on the estimation around point ω = [0, . . . , 0] (which can always be achieved by applying inverse Hamiltonian dynamics [52]) and denote |c 0 = |ψ ω=0 for notational simplicity. Then for )|c 0 , and according to Eq. (19) the minimum achievable cost for a fixed code space H C is given by: From the above formulation it is clear that we may always reduce the code space H C to span{|c k } P k=0 without increasing the cost. Hence, the problem of optimization over both probes and errorcorrection protocols is now equivalent to identification of the set {|c k } P k=0 that minimizes the cost with the constraint that H C = span{|c k } P k=0 satisfies the QEC conditions. To solve this problem, it will be convenient to formally extend the Hilbert space H S ⊗ H A by tensoring it with a (P + 1)-dimensional reference space H R = span{|0 , . . . , |P } (see Fig. 2c). This reference space will be representing the effective evolution of the probe state that happens within the code space and it will allow us to encode QEC conditions in a compact and numerically friendly way.
First, we introduce a matrix Q ∈ L(H R ⊗ H S ) that represents a code which, more formally, will be written as . This matrix is proportional to the reduced density matrix of the maximum entangled state between H R and H C . By its construction Q ≥ 0 and contains all relevant information on the code states in H C .
Next, we introduce effective generators G R i acting on H R so that they represent properly the action of the physical generators on the code space The effective evolution generators are related with the Q matrix via: Note that the identity operator here acts on the reference space H R , and not on the ancillary space H A . Taking into account the orthonormality of |c k and the QEC condition Eq. (6), we obtain the following constraints on Q form an orthonormal basis of Hermitian operators in L(H S ) such that S = span R {1 d , (S i ) P i=P +1 }. Any non-negative Q satisfying Eqs. (23)-(24) has the following form: where ν i ∈ R and B i are Hermitian. Conversely, for any nonnegative defined Q ≥ 0, we can consider its purification |Q ∈ H R ⊗ H S ⊗ H A , which when written as |Q = P k=0 |k H R ⊗ |c k H S ⊗H A yields the code states |c k . Note that it implies that the rank of Q corresponding to the dimension of the ancillary space. It is always sufficient to assume the dimension of the ancillary space to be dim H A = (P + 1)d. Therefore {G R i } is an achievable set of effective generators in L(H R ) (satisfying the QEC condition) if and only if there exist such ν i ∈ R and B i , for which Q ≥ 0.
Finally, in order to have an explicit dependence of the cost on the total time parameter T , we introduce a matrix Γ = 1 , and we end up with:

Reduction to an SDP
In order to reformulate Eq. (26) as an SDP, we first show that we may assume without loss of generality that Γ √ W −1 ≥ 0. Note that for any full rank matrix Γ, the polar decomposition theorem implies that there always exists an orthonormal matrix O such that OΓ √ W −1 ≥ 0. Next, as Γ ij = Im[ i|G R j |0 ], multiplication Γ by O is equivalent to rotating the base in the reference space H R . Since, according to Eq. (22) such a rotation cannot change the non-negativity of Q and at the same time it does not affect the figure of merit Tr W (Γ T Γ) −1 , the statement is proven. To put Eq. (26) in the form of an SDP, we introduce a positive matrix K ∈ R P ×P and a positive real number w. Now, using the following two relations, we see that P min w = min Tr(K 2 ) = min Tr W (Γ T Γ) −1 in Eq. (8), making it equivalent to Eq. (26). Hence the problem takes the form of an SDP.

Discussion
It should be remarked that JNT-QEC do not contain SEP-QEC as a subclass. In SEP-QEC, unlike in JNT-QEC, the noises are not fully corrected in the entire space, and the decoherence is only avoided by choosing a properly mixed state input. In general one could combine both these approaches in a unified framework by dividing the set of all parameters into smaller subsets and then applying JNT-QEC for each of these subset separately-in this approach the SEP-QEC case would correspond to the situation where JNT-QEC optimization is applied to single parameter subsets. Such an optimization is in principle doable, but will involve much large numerical effort and it is not clear that it will lead to better protocols. It is also worth noting that, apart from the improved metrological performance provided by QEC protocols when dealing with noisy systems, the above algorithm is also applicable in the noiseless scenario when S = span R {1 d }. In such a situation no QEC is required (for simplicity we may still use H C for span{|c k } P k=0 , but no recovery operation or projection Π H C is needed during evolution), but the condition ω = [0, . . . , 0] (which is achievable by applying inverse Hamiltonian dynamics [52]) is still required, as otherwise the derivatives of the state may not scale linearly with T . In such situations, the solution of JNT-QEC yields an optimally ancilla-assisted sensing protcol under arbitrary system dynamics (Hamiltonitans) that resolves the potential incompatibility issues between sensing of different parameters. It should be stressed that our approach is universal and unlike existing approaches [43,52,53] does not assume any specific structures of the Hamiltonians.

Advantages of JNT-QEC over SEP-QEC
In this section we investigate the potential advantages JNT-QEC provide compared with SEP-QEC. The characteristic feature of the SEP-QEC is the division resources so that each part of the resources is used to measure a given parameter independently of the others. This is reflected in the form of the input probe state ρ in = 1 P P i=1 |ψ i ψ i |. As a consequence we effectively measure each parameter only once in every P repetitions of an experiment (corresponding to the 1/P factor in the ρ in ). Therefore for a fixed total number of measurements, the uncertainty of estimating a given parameter will grow proportionally to P . Contrastingly, in JNT-QEC there is a chance to avoid the division of resources and use a single pure state and a single measurement to estimate all of them simultaneously. If there existed a code space, a state and a measurement which were all simultaneously optimal for all parameters, then the decrease the cost by a factor P would be achievable. This would be the largest possible advantage offered by JNT-QEC. In general we can write (see Appx. A for more formal derivation):

This inequality may be saturated only in special examples.
It is important to note that Theorem 2 gives us a recipe for identification of the optimal JNT-QEC, while the explicit construction of SEP-QEC discussed in the proof of Theorem 1 was only aimed at demonstrating the possibility of the HS and hence the protocol was not optimized. Therefore, to make a fair comparison between the two approaches, we need to compare the performance of the best JNT-QEC with the best SEP-QEC.
Below we present a way to find the optimal SEP-QEC, i.e. the protocol for which all parameters are measured independently on mutually orthogonal subspaces. We will also provide a useful lower bound for the minimal cost achievable in such protocols.

Optimization of SEP-QEC
Adapting the results from [36] (see the section named "Geometrical picture"), we can infer that for a given set of {G i }, the minimal variance achievable in estimation of a single parameter ω i in SEP-QEC is given by: where · denotes operator norm. Let |ψ i be the optimal state for measuring ω i . Therefore, Wii Fi 2 . Next, we may still improve the peromance of SEP-QEC by choosing a different QEC code based on a new set of generators obtained via a linear transformation on the parameters A ∈ R P ×P : so that the cost function remains unchanged. Note also, that rescaling any generator by constant factor has no impact on the results, so we may restrict to linear transformations satisfying (AA T ) ii = 1. Therefore the cost for the optimal SEP-QEC is given by: We want to stress that introducing this (relatively complicated) procedure of SEP-QEC optimization is necessary in order to distniguish the cases where the true advantage is offered by the joint multi-parameter approach compared with the situation where the advatage is only apparent and results from a suboptimal choice of separate protocols. Formula Eq. (32) is rather complicated and hard to compute in general. However, in case W = 1 the following lower bound is valid: where λ * ± are extreme eigenvalues of G * . Since the reasoning leading to Eq. (29) holds for any set of generators, not necessary the ones optimal for SEP-QEC, therefore, the optimal joint-estimation cost, in case W = 1, may be bounded as (see Appx. A for a formal derivation): where λ i± are extreme eigenvalues of G i . In order to provide the reader with some intuition on the concept presented above, we discuss below two extreme cases illustrating the apparent and maximal advantage of joint-estimation protocols over separate ones. These examples deal with noiseless scenarios where JNT-QEC is always optimal and are aimed to prepare the reader for more physical noisy examples discussed in Sec. 6. Still, in order not to introduce additional abbreviations we will still use the acronyms JNT-QEC, SEP-QEC to describe separate and joint estimation schemes, even though the role of QEC in these examples is trivial.

Example: Apparent advantage of JNT-QEC
Consider a noiseless physical system comprising P = 2 r qubits (for technical reasons we require the number of qubits to be the integer power of 2). Each qubit, regarded as a spin 1/2, senses the z-component of a local magnetic field which is assumed to be independent for different qubit locations. The corresponding sensing Hamiltonian for this system reads: (35) and the cost matrix is assumed to be W = 1. The minimum variance of estimating each parameter independently is lower bounded by ∆ 2ω i ≥ 1 T 2 (λ+−λ−) 2 . As the maximum and the minimum eigenvalues of each G i are λ i± = ±1, the minimal variance reads ∆ 2ω i = 1 4T 2 . According to Eq. (34)) this implies that the optimal JNT-QEC cost is lower bounded by P 4T 2 . Moreover, since the state |ψ = 1 ⊗P is simultaneously optimal for all the parameters (and the tensor structure guarantees no measurement incompatibility issue), this bound is saturable and hence On the other hand, the estimation strategy where each parameter ω i is estimated separately leads to the cost equal P 2 4T 2 . Therefore, one may naively think that it is an example of superiority JNT-QEC over SEP-QEC.
However, the separate protocol may be significantly improved here, by estimating different combinations of the parameters. According to Eq. (33), we should be looking for a proper linear combination of G i with the biggest difference of extreme eigenvalues. The most obvious choice is √ P , and therefore ∆ 2 WωSEP ≥ P 4T 2 , which is exactly equal to ∆ 2 WωJNT ! Below we show, that this bound may be saturated. Let the matrix A defining the transformation to the new set of parameters be proportional to a Hadamard matrix of size P × P (i.e. square matrix whose entries are either +1 or −1 and whose rows are mutually orthogonal) defined as follows: where the indices of the matrix are written using the binary representation, i = r−1 k=0 i k 2 k , where i k represent the binary digit of i at position k. Then for G = AG we have: As a result every ω i may be measured separately on subspace span{|λ i+ , |λ i− } with precision 1 4P T 2 . Therefore, for the initial sate ρ in = 1 which is exactly the same as the optimal ∆ 2 WωJNT . We see that the apparent advantage of a JNT-QEC disappears once a proper combinations of parameters are estimated in the SEP-QEC.

Example: Maximal advantage of JNT-QEC
As a contrasting example, here we discuss a situation where the advantage of JNT-QEC is genuine and is the maximal possible. This example also provides an intuition, what relation between the generators G i is responsible for this. Consider a noiseless system H S = span{|i } P i=0 with Hamiltonian: We focus on the estimation around point ω = [0, . . . , 0] with the cost matrix W = 1. As all generators are orthonormal Tr(G i G j ) = δ ij , the difference between extreme eigenvalues of any normalized combination of them, such as G * in Eq. (33), is at most equal to √ 2. Therefore, the optimal SEP-QEC cost is bounded by Importantly, the state |ψ in = |0 is simultaneously optimal for measuring all the parameters. There is also no measurement incompatibility problem as the optimal measurement for all the parameters is the measurement in the basis |i . Therefore, is achievable, and hence the there is a factor P decrease in the cost in case of JNT-QEC compared with the optimal SEP-QEC protocols.

Examples
Here we provide representative examples of a large class of multi-parameter estimation models, where there are unavoidable tradeoffs in determining the optimal states and measurements that arise due to the multi-parameter nature of the problem. Our methods are useful in identifying optimal strategies in all such models, provided the structure of noise admits the achievability of the HS via application of the most general QEC schemes. By the construction of our algorithm (Theorem 2), we have the guarantee that the solutions found are the optimal ones. For an interested reader, a broader discussion and generalizations of the results presented in this section, including proofs and analytical constructions of the codes are provided in Appx. B, Appx. C, Appx. D.

Single qubit case
Consider first the simplest single-qubit case with d = 2. The HS is achievable via QEC only in the case of single-rank Pauli noise (specified by a single Hermitian Lindbladian L) [34]. Without loss of generality we can set L = σ z (the Pauli-Z matrix). Since S = span{1, σ z }, at most two parameters may be estimated in a qubit system with the HS (as dim(S ⊥ ) = 2). However, it turns out that when the multi-parameter HNLS condition is met there is no benefit in performing the more sophisticated JNT-QEC compared to SEP-QEC, which is shown analytically in Appx. B.

Two qubits in a magnetic field
In order to appreciate the superiority of JNT-QEC over SEP-QEC, let us consider a two-qubit model which is a multi-parameter generalization of the one from [37]. Consider two localized qubits, coupled to a magnetic field, which is constant in both time and space, apart from some small fluctuations in the z direction. These fluctuations are assumed to be uncorrelated in time, but maximally anticorrelated in space (for the two qubits they have always opposite signs). Such a system may be effectively described by Eq. (1) with H = 1 that the minimal cost when each parameter is estimated with the optimal individual parameter strategy is ∆ 2ω x,y,z = 1 4T 2 . At the same time, in accordance with the discussion presented in Sec. 5, the best precision achievable using SEP-QEC for the standard cost matrix W = 1, is C for the formal derivation of the above formulas. Below we present the result of numerical optimization of the JNT-QEC approach (found by the algorithm presented in Theorem 2 and reconstructed to its analytical form). We will use the standard Bell states notation: Entanglement with ancilla will be abbreviated in the subscript |ψ ⊗ |i A ≡ |ψ i . Using the numerical algorithm we have found out, that the optimal code space has the form where the input state is |ψ ω=0 = |c 0 . Note that the presence of the last term in |c 3 (entangled with |4 A ) is necessary to satisfy the QEC conditions. The value of ϕ can be found analytically and the minimal total cost of estimation ∆ 2 Wω is achieved for: while the corresponding optimal cost is: As we see a significant improvement has been achieved here compared to SEP-QEC.

SU (d) generators' estimation
Finally let us consider an example of estimating parameters associated with SU (d) generators which shows an asymptotic advantage (with the number of parameters) of the JNT-QEC over the SEP-QEC protocol. Let H S be d-dimensional Hilbert space, which for a more intuitive notation may be regarded as the one associated with a spin-j particle (where d = 2j + 1).
First, let us recall the noiseless case, where the Hamiltonian H = form an orthonormal basis of Hermitian operators on H S . We focus on the estimation around point ω = [0, . . . , 0] with the cost matrix W = 1. Since all G i are orthonormal the maximum spread between their maximum and minimal eigenvalues as well as any normalized combination of them is √ 2. Using the same line of reasoning as presented in Sec. 5.3 we conclude that The bound for ∆ 2 WωJNT may be derived analytically. To achieve that, we use the following chain of inequalities involving the QFI matrix: where the first one is the CR inequality and the rest are general algebraic properties of positive semidefinite matrices. What remains to be done is to derive a proper bound for the trace of the QFI matrix. For any input state |ψ in ∈ H S ⊗ H A we have: Taking into account the normalization of G i and noting that is the Casimir operator of the SU (d) algebra, and as such is proportional to the identity, we get that Therefore: After substituting the above to Eq. (48) we get The example of a state which saturates the above bound is [52,53]. For such a state the QFI matrix is given by F ij = δ ij 4T 2 d , so the second and the third inequalities in Eq. (48) become equalities. As Im( ψ ω |Λ i Λ j |ψ ω ) ∝ ψ ω |[G i ⊗ 1, G j ⊗ 1]|ψ ω = 0, the first one (the CR bound) is saturable as well. Note that the role of ancillae here is to make optimal measurements with respect to different parameters compatible. As a result, we see that optimal ∆ 2 WωJNT is approximately 2 √ P times smaller than ∆ 2 WωSEP . Let us now consider a noisy version with a single Lindblad operator J z = j k=−j k |k k|. From Theorem 1 we know that only parameters associated with generators G i / ∈ span R {1, J z , J 2 z } may be estimated with the HS. Therefore we consider the Hamiltonian H = , with the standard cost matrix W = 1 (such cost makes the problem independent on choosing peculiar set of {G i }). In Fig. 3, we present numerical results for such a problem, and we observe a significant advantage over the SEP-QEC protocol as well as strong indication of the asymptotic P 3/2 scaling identical to the noiseless case. Even though the optimal JNT-QEC code cannot be written down analytically in a concise way, in Appx. D we provide an analytical suboptimal construction achieving the P 3/2 scaling, supporting the numerical results. The scaling advantage we prove here is not trivial because there are no decoherence-free subspaces in the system.

Conclusions and outlook
We have generalized previous results on single-parameter error-corrected metrology to the multiparameter scheme, obtaining a necessary and sufficient condition (HNLS) for the achievability of the HS. In case of scenarios when HNLS is satisfied, we developed an efficient numerical algorithm (formulated as an SDP) to find the optimal QEC protocol, including the optimal input state, QEC codes and measurements. Our algorithm is applicable to arbitrary system dynamics as long as the HS is achievable (including noiseless cases), which contrasts previous works where special dynamics of quantum system is assumed [43,52,53] or the estimation is performed on fixed quantum states [66]. However, it still remains open in which situations the requirement of noiseless ancillae could be removed [37,38] and whether QEC is still helpful in the case when HNLS is violated, as in the single-parameter case [41,74]. We also remark that our way of formulating the Knill-Laflamme conditions [72] as a positive semidefinite constraint Eq. (25) is novel and may have applications beyond error-corrected quantum metrology.
Finally, we note that in this paper we have followed the frequentist estimation approach, and in principle more stringent HS bounds might be derived when following the Bayesian approach as was demonstrated recently in the single parameter case [75].

WωSEP
In the reasoning presented in the beginning of Sec. 5 we have assumed that, when estimating a single parameter, the minimal variance of the estimator ∆ 2ω i achievable in SEP-QEC is no bigger than the one achievable in JNT-QEC (even if we focus only on this particular parameter). This statement is not so obvious as both protocols impose different constraints on the code space: where none of them is weaker or stronger than the other one.
Therefore, here we show directly that having JNT-QEC with accuracy ∆ 2 WωJNT , one may always construct SEP-QEC giving ∆ 2 WωSEP = P · ∆ 2 WωJNT . Let H C be the code space used in JNT-QEC. Without loss of generality let us assume that W is diagonal (otherwise we apply a transformation of parameters that diagonalizes W ) and dim(H C ) ≥ 2P + 1 (otherwise we trivially extend it using an additional ancilla). From the Matsumoto bound (see Eq. (11)) the optimal cost is given by: where for ω = 0 we have ∂ j |ψ ω ω=0 = −iT G C j |ψ ω=0 . We define H Ci = span R {|c i 0 , |c i 1 } with: where |x i are the result of optimization Eq. (53). As |ψ ω=0 , |x i ∈ H C , obviously the noise acts trivially on H Ci . Moreover: Therefore SEP-QEC with initial state ρ in = 1 P P i=1 |c i 0 c i 0 | and H Ci defined above leads to the total cost ∆ 2 WωJNT . From this reasoning we see clearly that the largest possible advantage of JNT-QEC over SEP-QEC is to decreasing the total cost by a factor P .

B Optimality of SEP-QEC in the single-qubit model
Let us consider the most general case of a two-parameter quantum estimation problem under Markovian noise in a 2-dimensional Hilbert space when the HS achievability condition is satisfied. Without loss of generality we assume that the Lindblad operator is L = σ z , the Hamiltonian is H = ω x σ x + ω y σ y and the cost matrix is diagonal: (otherwise, one can always apply a proper transformation in the parameter space ω = ωA −1 which diagonalizes the cost matrix, without changing orthonormality of the generators). Below we show that for such a problem there is no advantage of JNT-QEC over SEP-QEC. First let us consider the optimal SEP-QEC. Each generator has eigenvalues |λ x/y± = ±1 and each parameter may be estimated with precision ∆ 2ω x/y = 1 4T 2 [34]. Moreover, as the cost is diagonal, there is no point in applying an additional transformation in the optimization procedure given in Eq. (32)-indeed, any matrix A satisfying ∀ i (AA T ) ii = 1 preserves the eigenvalues of the generators as well as Tr(AW A T ). Therefore the optimal SEP-QEC cost is simply given by . We will show below that this cannot be outperformed by the optimal JNT-QECs. For the sake of notation simplicity, when tensoring an operator acting on H S with the identity on H A , the part ⊗1 will be omitted and we will denote σ i ⊗ 1 simply as σ i (unless this leads to ambiguity). First, we note that the diagonal elements of the QFI matrix for the state |ψ ω are Moreover Therefore we may focus on an upper bound for i=x,y ψ ω |σ i Π H C σ i |ψ ω . Let {|c 0 , |c 1 , |c 2 } be an orthonormal basis of H C ⊆ H S ⊗ H A . These vectors can be written down as where |A i 0/1 are normalized states in H A and ϕ i ∈ [0, π 2 ] (a potential complex phase is incorporated in the definition of|A i 0/1 ). The QEC condition requires ∀ i,j c i |σ z |c j = λδ ij , which leads to the following two constraints: (i) ∀ i cos 2 (ϕ i ) − sin 2 (ϕ i ) = λ means all ϕ i are equal (therefore superscript i will be omitted); Note that there is no fixed relationship between sets {|A i 0 } i=2 i=0 and {|A i 1 } i=2 i=0 -in particular it may happen that span{|A i 0 } = span{|A i 1 }. Effective generators in the chosen basis are given as: We focus on estimation around point ω = [0, ..., 0] for which |ψ ω = |ψ in = |c 0 . Then Since for each k = 0/1, states {|A i k } 2 i=0 are mutually orthonormal, i=x,y where the first inequality is saturated if and only if both |A 0 0 ∈ span{|A i 1 } 2 i=0 and |A 0 1 ∈ span{|A i 0 } 2 i=0 . Using Eq. (57)  . This implies that the JNT-QEC stategy cannot outperform the best SEP-QEC strategy.
C The optimal SEP protocol for sensing all magnetic field components in presence of correlated dephasing noise in the two-qubit model Here we prove formally, that the optimal precision achievable in the second example from the main text is ∆ 2 WωSEP = 9 4T 2 . Let us briefly review the problem. We consider a two-atom system with Hamiltonian H = acts on the i th atom) and a single Lindblad operator As for all three generators the minimal and the maximal eigenvalues are respectively −1 and +1, we immediately see that ∆ 2ω x,y,z ≥ 1 4T 2 and from that ∆ 2 WωSEP ≥ P 2 4T 2 = 9 4T 2 (which cannot be improved by applying transformation A). Below we show that such a precision is indeed achievable.
First, ω z can be estimated using a decoherence-free subspace [76] span{|ψ z , 1 2 (σ (1) z +σ (2) z ) |ψ z }, where |ψ z = 1 √ 2 (|↑↑ + |↓↓ ), which leads to the precision ∆ 2 ω z = 1 4T 2 . In case of ω x the situation is slightly more complicated, as using the analogue approach the subspace span{|ψ x , 1 2 (σ (1) x + σ (2) x ) |ψ x } would not satisfy the QEC conditions. To get the desired estimation precision we need to find the state |ψ x which is optimal for measuring ω x (from the point of view of noiseless, single-parameter estimation) and for which satisfies the QEC condition: and, moreover, other generators act trivially inside H Cx : It is known that for single-parameter frequency estimation the optimal state corresponds to an equally weighted superposition of states with minimal and maximal eigenvalues of the generator. For 1 2 (σ (1) x + σ (2) x ) it will be: for any ϕ ∈ R. Note, that any superposition of |ψ ϕ x (with different ϕ) entangled with separated ancillae is still optimal for sensing ω x . Therefore, we can take |ψ x = 1 √ 2 (|ψ 0 x |1 A + |ψ π x |2 A ). We have: and then the code space Eq. (66) satisfies Eq. (67) and Eq. (68). It gives as ∆ 2 ω x = 1 4T 2 . Analogous reasoning could be provided for ω y . Therefore for ρ in = 1 3 i=x,y,z |ψ i ψ i | we have ∆ 2 Wω = 9 4T 2 in line with general considerations on the performance of the SEP-QEC codes as given in Sec. 5.

D Estimating the SU (d) generators
Below we present an example of a JNT-QEC protocol allowing one to achieve the total cost ∆ 2 Wω = Θ(P 3 The QFIs are which in our case simplifies to: As ψ R |[(G R kl ) H C R , (G R k l ) H C R ]|ψ R = 0, the CR bound is saturable and the total cost is k>l ∆ 2 ω R kl = 2j + 1 4 k>l 1 p 2 = 3j(2j + 1) 2 2 = Θ(P Imaginary off-diagonal generators. The reasoning is analogous to the previous case. Note that using different ancillary spaces for real and imaginary generators is needed. Even though ψ R |G R kl J z G R k l |ψ R ∝ δ kl,k l and ψ R |G I kl J z G I k l |ψ R ∝ δ kl,k l are satisfied automatically, ψ R |G R kl J z G I k l |ψ R ∝ δ kl,k l may not be true. Diagonal generators. The number of diagonal generators scales like Θ(j) (whereas for offdiagonal elements the scaling is Θ(j 2 )), implying that the estimation with respect to diagonal generators does not contribute significantly to the overall scaling. Therefore we could simply use the SEP-QEC approach. Following [36], any traceless generator may by written down as: We define states |c i+ , |c i− as purifications of these density matrices by using mutually orthogonal ancillary subspaces H Ai+ , H Ai− : and from that for code space span{|c i+ , |c i− } and input state 1 √ 2 (|c i+ + |c i− ) we have F ωi ≥ 1 2j+1 . For a single-parameter problem the CR bound is always saturable so using the SEP-QEC approach for the input state we have Results. Finally, combining all above, for the input state: (with properly applied QEC protocol) we get P i=1 ∆ 2 ω i = Θ(P