Unified and computable approach to optimal strategies for multiparameter estimation

Abstract

Precise estimation of physical parameters underpins both scientific discovery and technological development. A central goal of quantum metrology and sensing is to exploit quantum resources like entanglement to devise optimal strategies for estimating physical parameters as precisely as possible. While substantial progress has been made in single-parameter quantum metrology, the multiparameter scenario remains significantly more challenging due to the issue of parameter incompatibility. In this work, we present a unified and computable approach for the simultaneous estimation of multiple parameters that attains the ultimate precision permitted by quantum mechanics. The core of our approach is to integrate the quantum tester formalism into the recently proposed tight Cramér-Rao type bound. This formulation enables us to figure out the highest achievable precision via upper and lower bounds that are computable via semidefinite programs. More importantly, within this formulation, diverse quantum resources, including entanglement, coherence, quantum control, and indefinite causal order, are treated on equal footing and systematically optimized for the purpose of achieving the ultimate precision in multiparameter estimation. As a result, our approach is applicable to various metrological strategies both in the presence and absence of noise. To demonstrate its utility, we revisit three-dimensional magnetic-field estimation, uncovering the strengths and limitations of existing analytical results and further establishing a strict hierarchy among different types of strategies.

I Introduction

Much of quantitative science deals with measuring unknown parameters of a physical process. With the advent of quantum technologies, it has been found that quantum resources such as entanglement and coherence allow for pushing the measurement precision limit beyond what is achievable with classical means Giovannetti et al. 2004; Giovannetti et al. 2006. A central goal in quantum metrology and sensing is therefore to devise optimal strategies for exploiting quantum resources to estimate parameters as precisely as possible, which has motivated vibrant research activities over the past two decades Giovannetti et al. 2011; Pezzè et al. 2018.

The problem of devising optimal strategies has been well addressed in the single-parameter setting Liu et al. 2021; Liu et al. 2024. This is largely because the ultimate precision achievable is fully characterized by the quantum Cramér–Rao bound, which can be saturated by suitable measurements under broad conditions Braunstein and Caves 1994; Zhang et al. 2015; Lu et al. 2015; Zhang and Gong 2020; Zhong et al. 2020; Xu et al. 2020; Li et al. 2021; Zhang and Tong 2022; Zhang and Tong 2024; Zhou and Zhang 2025; Zhang and Tong 2025. This theoretical achievement has led to a variety of experimentally relevant protocols realizing quantum-enhanced measurements in systems ranging from optical interferometry to atomic clocks and solid-state sensors Jiao et al. 2023; Pei et al. 2025; Montenegro et al. 2025. The situation, however, changes dramatically in the multiparameter regime. When several parameters are estimated simultaneously, the optimal measurements associated with different parameters may be mutually incompatible, preventing the simultaneous saturation of the single-parameter quantum Cramér–Rao bounds. This issue, known as parameter incompatibility, renders the characterization of ultimate precision substantially more subtle and has been recognized as a central obstacle in multiparameter quantum metrology Lu and Wang 2021; Albarelli and Demkowicz-Dobrzański 2022. As such, identifying both the fundamental precision limits and the corresponding optimal strategies remains an open challenge in multiparameter scenarios.

The past few years have witnessed considerable efforts devoted to addressing this challenge, driven by the growing relevance of multiparameter estimation in applications like quantum imaging Tsang et al. 2016; Lupo and Pirandola 2016; Chrostowski et al. 2017; Chrostowski et al. 2017; Řehaček et al. 2017, magnetic field sensing Baumgratz and Datta 2016; Hou et al. 2020; Zhuang et al. 2024, and Hamiltonian parameter estimation Yuan 2016. It has become increasingly clear that quantum-enhanced measurements may be enabled by a variety of quantum resources, such as entanglement Giovannetti et al. 2004, coherence Giovannetti et al. 2006, quantum control Yuan 2016, and even indefinite causal order Zhao et al. 2020; Liu et al. 2023. While each of these ingredients has been demonstrated to provide advantages in specific scenarios, their roles in multiparameter estimation are often analyzed in isolation, making a systematic and unified treatment elusive. Furthermore, the presence of noise, which is ubiquitous in practical applications, substantially complicates the identification of truly optimal strategies for multiparameter estimation. Therefore, a unified framework capable of systematically optimizing over all available quantum resources is still lacking in multiparameter scenarios up to now.

In this work, we present a unified and computable approach for the simultaneous estimation of multiple parameters that attains the ultimate precision permitted by quantum mechanics. Our approach is grounded in the tight Cramér-Rao type bound (TCRB) recently introduced by Hayashi and Ouyang Hayashi and Ouyang 2023, which fully characterizes the highest achievable precision in estimating multiple parameters from a given output state. The new contribution of our work lies in integrating the quantum tester formalism Chiribella et al. 2008; Chiribella et al. 2009; Araújo et al. 2015 into the TCRB. Specifically, we formulate the problem of identifying the highest achievable precision as a conic programming problem, where the optimization is performed over all admissible quantum testers within a specified type of strategies. This formulation enables us to figure out the highest achievable precision via upper and lower bounds that are computable via semidefinite programs (SDPs). More importantly, within this formulation, diverse quantum resources, including entanglement, coherence, quantum control, and indefinite causal order, are treated on equal footing and systematically optimized for the purpose of achieving the ultimate precision in multiparameter estimation. As a result, our approach is applicable to various metrological strategies both in the presence and absence of noise.

We clarify that, although the derivations presented below closely follow those in Ref. Hayashi and Ouyang 2023, the results obtained here are fundamentally interesting and bear significant physical implications. Indeed, the TCRB was formulated to find optimal measurements given a fixed output state, without considering how this state is generated Hayashi and Ouyang 2023; Not. Consequently, the TCRB alone cannot determine the optimal strategies, as it does not account for how the output states are generated. In contrast, by incorporating quantum testers into the TCRB, our approach provides a complete characterization of the entire estimation protocol, encompassing state preparation, parameter encoding, intermediate operations, and final measurements. As such, all quantum resources are fully taken into account in our approach. To demonstrate the physical relevance of our approach, we apply it to the three-dimensional magnetic-field estimation Baumgratz and Datta 2016; Hou et al. 2020, a paradigmatic and fundamentally important task in quantum metrology. We first analyze the noiseless setting and show that existing analytical results exhibit both notable strengths and intrinsic limitations. We then turn to the noisy regime and establish a strict hierarchy among different types of strategies. Our approach thus represents a valuable tool capable of providing deeper physical insight into multiparameter estimation.

This paper is organized as follows. In section II, we recall the quantum tester formalism. In section III, we present the unified formulation with the aid of quantum testers. In section IV and section V, we derive upper and lower bounds, respectively. We then revisit analytical results for the three-dimensional magnetic field estimation in section VI and establish a strict hierarchy among different types of estimation strategies in section VII. We conclude this work in section VIII.

II Preliminaries

Figure 1: Schematic of quantum testers with $N=2$. The dashed boxes in the left panels are concrete realizations of estimation strategies: (a) parallel strategies, (b) sequential strategies, (c) causal superposition strategies, and (d) general indefinite-causal-order strategies. Here $\rho$ is the input state, $\mathcal{E}_{\bm{\theta}}$ represents the parameter-encoding channel, and $U_{i}$ denotes the unitary control operations. In each of these strategies, a final measurement $\{\Pi_{x}\}_{x}$ is performed on the output state. The black and red curves in (c) represent two different causal orders. The right panels display the associated quantum testers.

To present our results clearly, we need to recapitulate the quantum tester formalism Chiribella et al. 2008; Chiribella et al. 2009; Araújo et al. 2015. Suppose that there are $N$ identical quantum channels $\mathcal{E}_{\bm{\theta}}$ at our disposal. Here, $\bm{\theta}=(\theta_{1},\theta_{2},\cdots,\theta_{\mathfrak{p}})^{T}$ collectively denotes the unknown parameters, where $\mathfrak{p}$ is the number of parameters in question. Let $\mathcal{H}_{I_{k}}$ and $\mathcal{H}_{O_{k}}$ represent the input and output Hilbert spaces of the $k$th copy of the channel, respectively. Denote by $\mathcal{L}(\mathcal{H})$ the set of linear operators over a Hilbert space $\mathcal{H}$. The $k$th copy of the channel $\mathcal{E}_{\bm{\theta}}$, as a linear map from $\mathcal{L}(\mathcal{H}_{I_{k}})$ to $\mathcal{L}(\mathcal{H}_{O_{k}})$, can be associated with a positive semidefinite operator in $\mathcal{L}(\mathcal{H}_{I_{k}}\otimes\mathcal{H}_{O_{k}})$

$$E_{\bm{\theta}}\coloneqq\mathrm{id}\otimes\mathcal{E}_{\bm{\theta}}\left(|I\rangle\hskip-2.7pt\rangle\langle\hskip-2.7pt\langle I|\right),$$ (1)

known as the Choi-Jamiołkowski (CJ) operator Jamiołkowski 1972; Choi 1975. Here, $\mathrm{id}$ denotes the identity map and $|I\rangle\hskip-2.7pt\rangle=\sum_{j}|j\rangle|j\rangle$. The CJ operator corresponding to the $N$ identical channels is

$$C_{\bm{\theta}}\coloneqq E_{\bm{\theta}}^{\otimes N}\in\mathcal{L}\left(\mathcal{H}_{IO}\right),$$ (2)

where $\mathcal{H}_{IO}\coloneqq\mathcal{H}_{I_{1}}\otimes\mathcal{H}_{O_{1}}\otimes\cdots\otimes\mathcal{H}_{I_{N}}\otimes\mathcal{H}_{O_{N}}$. A quantum tester is a generalization of a quantum measurement and may account for various ingredients of an estimation strategy, including the input state, quantum controls, and final measurement (see Fig. 1). Mathematically, a quantum tester is described by a collection of positive semidefinite operators $\{X_{x}\}_{x}$ in $\mathcal{L}(\mathcal{H}_{IO})$ such that

$$p_{\bm{\theta}}\left(x\right)\coloneqq\mathrm{tr}\left(C_{\bm{\theta}}X_{x}^{T}\right)$$ (3)

are probabilities, that is, $\sum_{x}{p_{\bm{\theta}}(x)}=1$ and $p_{\bm{\theta}}(x)\geq 0$ for all $x$. Here, the superscript $T$ denotes the matrix transpose.

III Unified formulation with quantum testers

To seek for the unified formulation, we exploit the quantum tester formalism to describe various types of strategies, including: ($i$) parallel strategies Giovannetti et al. 2006, where these $N$ channels are applied simultaneously on a multipartite entangled state [see Fig. 1(a)]; ($ii$) sequential strategies Giovannetti et al. 2006, involving successive queries of the channels, possibly interspersed with unitary control operations [see Fig. 1(b)]; ($iii$) causal superposition strategies Zhao et al. 2020, where the channels are probed in a superposition of different causal orders [see Fig. 1(c)]; ($iv$) general indefinite-causal-order strategies Liu et al. 2023, which encompass the most general causal relations among the channels and include causal superposition strategies as special cases [see Fig. 1(d)]. Each type of estimation strategies corresponds to a specific set of quantum testers. Hereafter, we use index $k$ to specify the type of strategies in question, that is, $k=i,ii,iii,iv$ refer to the mentioned four types of strategies ($i$)-($iv$), respectively. We use $X_{x}^{(k)}$ to denote the quantum tester associated with the strategies of type $k$ and define

$$X^{\left(k\right)}\coloneqq\sum_{x}{X_{x}^{\left(k\right)}}.$$ (4)

It is known Araújo et al. 2015; Zhou et al. 2024 that, for strategies of type $k=i,ii,iv$, the associated set of quantum testers can be characterized as

$$\mathcal{X}^{\left(k\right)}\!\coloneqq\!\left\{X^{\left(k\right)}|X^{(k)}\!\geq\!0,\Lambda^{\left(k\right)}\left(X^{\left(k\right)}\right)\!=\!X^{\left(k\right)},\mathrm{tr}\left(X^{\left(k\right)}\right)\!=\!d_{O}\right\}.$$ (5)

Here, $d_{O}=\mathrm{dim}\left(\mathcal{H}_{O_{1}}\otimes\cdots\otimes\mathcal{H}_{O_{N}}\right)$ and $\Lambda^{(k)}$ denotes a linear map from $\mathcal{L}(\mathcal{H}_{IO})$ to itself, whose explicit expression can be found in Appendix A. The case $k=iii$ can be treated in a similar way but with some additional care (see Appendix B).

To estimate $\bm{\theta}$, we adopt a locally unbiased estimator Helstrom 1976; Holevo 2011, denoted by $\hat{\bm{\theta}}(x)$, which satisfies the two conditions

$$\sum_{x}{p_{\bm{\theta}}(x)\hat{\bm{\theta}}(x)}=\bm{\theta},$$ (6)

$$\sum_{x}{\partial_{j}p_{\bm{\theta}}(x)\hat{\theta}_{i}(x)}=\delta_{ij},$$ (7)

for $i,j=1,\cdots,\mathfrak{p}$. Here, $\partial_{j}$ stands for the partial derivative with respect to $\theta_{j}$, $\hat{\theta}_{i}(x)$ denotes the $i$-th component of $\hat{\bm{\theta}}(x)$, and $\delta_{ij}$ is the Kronecker delta. The performance of the estimator is characterized by its covariance matrix

$$\Sigma=\sum_{x}{p_{\bm{\theta}}\left(x\right)\left[\hat{\bm{\theta}}(x)-\bm{\theta}\right]\left[\hat{\bm{\theta}}(x)-\bm{\theta}\right]^{T}}.$$ (8)

To quantify the estimation error, we need to find a scalar function of the covariance matrix. A widely used choice is the weighted trace of the covariance matrix, $\tr(W\Sigma)$, where $W$ is a real positive semidefinite matrix that reflects the relative importance of different parameters Helstrom 1976; Holevo 2011. So, we need to solve the following optimization problem:

$$\min\quad \mathrm{tr}\left(W\Sigma\right)$$ (9a)

$$\mathrm{s}.\mathrm{t}.\quad \sum_{x}{p_{\bm{\theta}}(x)\hat{\bm{\theta}}(x)}=\bm{\theta},$$ (9b)

$$\sum_{x}{\partial_{j}p_{\bm{\theta}}(x)\hat{\theta}_{i}(x)}=\delta_{ij},\quad i,j=1,\cdots,\mathfrak{p}.$$ (9c)

Note that $p_{\bm{\theta}}(x)$ is determined by the quantum tester $\{X_{x}^{(k)}\}_{x}$ via eq. 3. The optimization in eq. 9 is thus over all admissible quantum testers of type $k$ and all locally unbiased estimators. Note also that we can remove the first type of constraints in eq. 9 without affecting the final result Demkowicz-Dobrzański et al. 2020. That is, eq. 9 can be simplified as

$$\min\quad \mathrm{tr}\left(W\Sigma\right)$$ (10a)

$$\mathrm{s}.\mathrm{t}.\quad \sum_{x}{\partial_{j}p_{\bm{\theta}}(x)\hat{\theta}_{i}(x)}=\delta_{ij},\quad i,j=1,\cdots,\mathfrak{p}.$$ (10b)

The optimal value of eq. 10 characterizes the highest achievable precision in estimating $\bm{\theta}$ via the strategies of type $k$.

To rewrite eq. 10 in a convenient form, we introduce an argumented weight matrix $\tilde{W}$ defined over the space $\mathbb{C}^{\mathfrak{p}+1}$

$$\tilde{W}\coloneqq 0\oplus W,$$ (11)

i.e., the direct sum of $0$ and $W$. We define the vectors in $\mathbb{C}^{\mathfrak{p}+1}$

$$|h_{x}\rangle\coloneqq\left(\begin{array}[]{c}1\\ \hat{\bm{\theta}}(x)-\bm{\theta}\\ \end{array}\right).$$ (12)

Using Eqs. (11) and (12) and noting that $p_{\bm{\theta}}(x)=\mathrm{tr}\left(C_{\bm{\theta}}X_{x}^{\left(k\right)T}\right)$, we can rewrite the objective function in eq. 10 as

$$\mathrm{tr}(W\Sigma)=\mathrm{tr}\left[\tilde{W}\otimes C_{\bm{\theta}}\sum_{x}{|h_{x}\rangle\langle h_{x}|\otimes X_{x}^{(k)T}}\right].$$ (13)

Further, letting $A_{i}=\left(|0\rangle\langle i|+|i\rangle\langle 0|\right)\big/2$ with $\left\{|i\rangle\right\}_{i=0}^{\mathfrak{p}}$ denoting the standard orthonormal basis of $\mathbb{C}^{\mathfrak{p}+1}$, we can rewrite the left-hand side of the constraints in eq. 10 as

$$\sum_{x}{\partial_{j}p_{\bm{\theta}}(x)\hat{\theta}_{i}(x)}=\mathrm{tr}\left[A_{i}\otimes\partial_{j}C_{\bm{\theta}}\sum_{x}{|h_{x}\rangle\langle h_{x}|\otimes X_{x}^{\left(k\right)T}}\right].$$ (14)

Therefore, eq. 10 can be equivalently formulated as

$$\underset{X_{x}^{\left(k\right)},\ket{h_x}}{\min}\,\, \mathrm{tr}\left[\tilde{W}\otimes C_{\bm{\theta}}\sum_{x}{|h_{x}\rangle\langle h_{x}|\otimes X_{x}^{\left(k\right)}}\right]$$ (15a)

$$\mathrm{s}.\mathrm{t}.\,\, \sum_{x}{X_{x}^{\left(k\right)}}\in\mathcal{X}^{\left(k\right)},$$ (15b)

$$X_{x}^{\left(k\right)}\geq 0,\quad\forall x,$$ (15c)

$$\mathrm{tr}\left[A_{i}\otimes\partial_{j}C_{\bm{\theta}}\sum_{x}{|h_{x}\rangle\langle h_{x}|\otimes X_{x}^{\left(k\right)}}\right]=\delta_{ij},$$

$$\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad i,j=1,\cdots,\mathfrak{p},$$ (15d)

where $X_{x}^{(k)T}$ is replaced by $X_{x}^{(k)}$ as $\mathcal{X}^{(k)}$ is closed under transposition.

We observe that eq. 15 involves products of variables $\ket{h_x}\bra{h_x}$ and $X_{x}^{(k)}$, which render the optimization challenging. To address this issue, we introduce the separable cone

$$\mathcal{S}\coloneqq\mathrm{conv}\left\{P_{1}\otimes P_{2}|P_{1}\in\mathrm{Pos}\left(\mathbb{C}^{\mathfrak{p}+1}\right),P_{2}\in\mathrm{Pos}\left(\mathcal{H}_{IO}\right)\right\},$$ (16)

where $\mathrm{Pos(\mathcal{H})}$ denotes the set of positive semidefinite operators on $\mathcal{H}$, and conv represents the convex hull Boyd and Vandenberghe 2004. Identifying $\sum_{x}{|h_{x}\rangle\langle h_{x}|\otimes X_{x}^{\left(k\right)}}$ as an element $Y^{(k)}\in\mathcal{S}$, we can show that eq. 15 is equivalent to the optimization

$$\underset{Y^{(k)}\in\mathcal{S}}{\min}\,\, \mathrm{tr}\left[\tilde{W}\otimes C_{\bm{\theta}}Y^{(k)}\right]$$ (17a)

$$\mathrm{s}.\mathrm{t}.\,\, \mathrm{tr}\left[A_{i}\otimes\partial_{j}C_{\bm{\theta}}Y^{(k)}\right]=\delta_{ij},\,\,i,j=1,\cdots,\mathfrak{p},$$ (17b)

$$\mathrm{tr}_{1}\left[|0\rangle\langle 0|\otimes\mathbb{I}_{IO}Y^{(k)}\right]\in\mathcal{X}^{\left(k\right)},$$ (17c)

where $\tr_{1}$ is the partial trace over the first subsystem and $\mathbb{I}_{IO}$ is the identity operator on $\mathcal{H}_{IO}$. The equivalence between eq. 15 and eq. 17 can be established as follows. Note that any feasible variable $\sum_{x}{|h_{x}\rangle\langle h_{x}|\otimes X_{x}^{\left(k\right)}}$ in eq. 15 is also a feasible variable in eq. 17. This implies that the optimal value of eq. 17 is no larger than that of eq. 15. On the other hand, given any feasible variable $Y^{(k)}$ in eq. 17, we can decompose it in the form

$$Y^{\left(k\right)}=\sum_{x}{\ket{\phi_{x}}\bra{\phi_{x}}\otimes P_{x}},$$ (18)

where $\ket{\phi_{x}}$ is a vector in $\mathbb{C}^{\mathfrak{p}+1}$ and $P_{x}\in\mathrm{Pos}(\mathcal{H}_{IO})$. Since $A_{i},|0\rangle\langle 0|$, and $\tilde{W}$ are all real symmetric matrices, only the real part of $\ket{\phi_{x}}\bra{\phi_{x}}$ contributes to the constraints and the objective function. Therefore, without loss of generality, we can assume that $\ket{\phi_{x}}$ is real. Besides, note that the components $\ket{\phi_{x}}\bra{\phi_{x}}\otimes P_{x}$ in the sum decomposition (18) do not contribute to the constraints and the objective function in Eq. (17) when $\langle 0|\phi_{x}\rangle=0$. We can further assume $\langle 0|\phi_{x}\rangle\neq 0$ without loss of generality. Actually, we can assume that $\langle 0|\phi_{x}\rangle=1$ without loss of generality, since $\ket{\phi_{x}}$ can be rescaled arbitrarily by absorbing the scaling factor into $P_{x}$. Now we can identify $\ket{\phi_{x}}$ with $|h_{x}\rangle$ defined in eq. 12 and $P_{x}$ with $X_{x}^{\left(k\right)}$; that is, a feasible variable $Y^{(k)}$ in eq. 17 corresponds to a feasible variable $\sum_{x}{|h_{x}\rangle\langle h_{x}|\otimes X_{x}^{\left(k\right)}}$ in eq. 15. It follows that the optimal value of eq. 15 is no larger than that of eq. 17. This completes the proof of the equivalence.

We highlight that Eq. (17) is the unified formulation we seek for. Notably, once eq. 17 is solved, we can determine the highest achievable precision in the strategies of type $k$, and moreover, we may construct an optimal strategy that attains this precision. Specifically, associated with an optimal solution of the form (18), the optimal estimator is $\hat{\theta}_{i}(x)=\theta_{i}+\langle i|\phi_{x}\rangle$ for $i=1,\cdots,\mathfrak{p}$. Meanwhile, the optimal quantum tester can be identified as $\{{P}_{x}\}_{x}$. It is worth noting that the physical realizations of quantum testers have been extensively studied. In particular, a general scheme for implementing quantum testers corresponding to parallel, sequential, and causal superposition strategies was developed in Refs. Chiribella et al. 2009; Bavaresco et al. 2021. Unfortunately, a concrete implementation for general indefinite-causal-order strategies remains an open problem Kurdziałek et al. 2023; Bavaresco et al. 2021. Below, we show how to effectively solve eq. 17 by constructing some upper and lower bounds that are computable via SDPs.

IV Semidefinite programs for upper bounds

To construct upper bounds for eq. 17, we randomly generate $m$ unit vectors $\left\{|w_{x}\rangle\right\}_{x=1}^{m}$ in $\mathbb{C}^{\mathfrak{p}+1}$. Actually, these vectors can be chosen to be real. Then, letting $Y^{(k)}=\sum_{x}{|w_{x}\rangle\langle w_{x}|\otimes X_{x}^{\left(k\right)}}$, where $X_{x}^{\left(k\right)}$ is a positive semidefinite matrix, we can convert eq. 17 into the form

$$\underset{X_{x}^{\left(k\right)}}{\min}\,\, \sum_{x}{\langle w_{x}|\tilde{W}|w_{x}\rangle\mathrm{tr}\left(C_{\bm{\theta}}X_{x}^{\left(k\right)}\right)}$$ (19a)

$$\mathrm{s}.\mathrm{t}.\,\, \sum_{x}{\left|\langle w_{x}|0\rangle\right|^{2}X_{x}^{\left(k\right)}}\in\mathcal{X}^{\left(k\right)},$$ (19b)

$$X_{x}^{\left(k\right)}\geq 0,\quad\forall x,$$ (19c)

$$\sum_{x}{\langle w_{x}|A_{i}|w_{x}\rangle\mathrm{tr}\left(\partial_{j}C_{\bm{\theta}}X_{x}^{\left(k\right)}\right)}=\delta_{ij},$$

$$\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad i,j=1,\cdots,\mathfrak{p},$$ (19d)

where the optimization variables are $\left\{X_{x}^{\left(k\right)}\right\}_{x=1}^{m}$. Apparently, eq. 19 provides an upper bound for eq. 17 and different sets of $\left\{|w_{x}\rangle\right\}_{x=1}^{m}$ correspond to different upper bounds. It follows from eq. 5 that Eq. (19b) can be explicitly rewritten as the following two linear constraints:

$$\mathrm{tr}\left[\sum_{x}{\left|\langle w_{x}|0\rangle\right|^{2}X_{x}^{\left(k\right)}}\right]=d_{O},$$ (20)

and

$$\left(\mathrm{id}-\Lambda^{(k)}\right)\left[\sum_{x}{\left|\langle w_{x}|0\rangle\right|^{2}X_{x}^{\left(k\right)}}\right]=0.$$ (21)

Therefore, the optimization problem in eq. 19 is a SDP and can be solved efficiently using numerical methods Grant and Boyd 2008; CVX Research 2012. Intuitively speaking, in the course of randomly generating more vectors $\left\{|w_{x}\rangle\right\}_{x=1}^{m}$, the set of operators $\left\{\sum_{x}{|w_{x}\rangle\langle w_{x}|\otimes X_{x}^{\left(k\right)}}|X_{x}^{\left(k\right)}\geq 0\right\}$ becomes larger, and consequently, this set could approximate the separable cone $\mathcal{S}$ increasingly well. It is therefore expected that the optimal value of eq. 19 converges to that of eq. 17 as $m$ increases with a high probability. We would like to mention that similar ideas have been employed to approximate a desired set of operators via random sampling Zhang et al. 2018.

V Semidefinite programs for lower bounds

To construct lower bounds for eq. 17, we resort to the technique of symmetric extension in entanglement theory Doherty et al. 2002; Doherty et al. 2004; Tavakoli et al. 2024. The main idea of this technique is to relax the set of separable states to a larger set that can be characterized by semidefinite constraints. Recall that, if $\rho$ is a separable state acting on the Hilbert space $\mathcal{H}_{1}\otimes\mathcal{H}_{2}$, then $\rho$ has an extension $\tilde{\rho}$ that acts on the extended Hilbert space $\mathcal{H}^{\otimes n}_{1}\otimes\mathcal{H}_{2}$ with the following two properties. First, $\tilde{\rho}$ is symmetric under interchanges of any two copies of subsystem $\mathcal{H}_{1}$. Second, tracing out any $n-1$ copies of subsystem $\mathcal{H}_{1}$ yields the original state $\rho$. These two properties can be easily verified by noting that $\tilde{\rho}=\sum_{i}p_{i}\rho_{i}^{\otimes n}\otimes\sigma_{i}$ if $\rho=\sum_{i}p_{i}\rho_{i}\otimes\sigma_{i}$, where $\{p_{i}\}$ is a probability distribution, and $\rho_{i}$ and $\sigma_{i}$ are two quantum states on $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, respectively. Moreover, if $\rho$ is entangled, there always exists one positive integer $n$ such that the above extension $\tilde{\rho}$ cannot be constructed. As the separate cone $\mathcal{S}$ shares the same structure as the set of separable states, we can employ the symmetric extension technique to obtain a series of converging lower bounds for eq. 17. Specifically, let $Y^{(k)}$ be a separable operator acting on the Hilbert space $\mathbb{C}^{\mathfrak{p}+1}\otimes\mathcal{H}_{IO}$. According to the symmetric extension technique, $Y^{(k)}$ has an extension $Y_{n}^{(k)}$ that acts on the extended Hilbert space $(\mathbb{C}^{\mathfrak{p}+1})^{\otimes n}\otimes\mathcal{H}_{IO}$ so that eq. 17 can be relaxed as follows:

$$\underset{Y^{(k)}_{n}}{\min}\quad \mathrm{tr}\left[\tilde{W}\otimes C_{\bm{\theta}}Y^{(k)}\right]$$ (22a)

$$\mathrm{s}.\mathrm{t}.\quad \mathrm{tr}_{1}\left[|0\rangle\langle 0|\otimes\mathbb{I}_{IO}Y^{(k)}\right]\in\mathcal{X}^{\left(k\right)},$$ (22b)

$$\mathrm{tr}\left[A_{i}\otimes\partial_{j}C_{\bm{\theta}}Y^{(k)}\right]=\delta_{ij},\,\,i,j=1,\cdots,\mathfrak{p},$$ (22c)

$$\big(U_{\pi}\otimes\mathbb{I}_{IO}\big)Y_{n}^{\left(k\right)}\big(U_{\pi}^{\dagger}\otimes\mathbb{I}_{IO}\big)=Y_{n}^{\left(k\right)},\quad\forall\pi\in\mathfrak{S}_{n},$$ (22d)

$$\mathrm{tr}_{1,\cdots,n-1}\left(Y_{n}^{\left(k\right)}\right)=Y^{\left(k\right)},$$ (22e)

$$Y_{n}^{\left(k\right)}\geq 0.$$ (22f)

Here the optimization variable is $Y_{n}^{(k)}$ and $Y^{(k)}$ is an intermediate variable. $\mathfrak{S}_{n}$ denotes the symmetric group of degree $n$ and $U_{\pi}$ is the unitary operator that permutes the $n$ copies of subsystem $\mathbb{C}^{\mathfrak{p}+1}$ according to the permutation $\pi$. Owing to the permutation symmetry constraint (22d), the partial trace in eq. 22e can be taken over any $n-1$ copies of subsystem $\mathbb{C}^{\mathfrak{p}+1}$. Here we set the partial trace to be over the first $n-1$ copies. According to the symmetric extension technique Doherty et al. 2002; Doherty et al. 2004; Tavakoli et al. 2024, the optimal value of eq. 22 is no larger than that of eq. 17 and it coincides with the latter for sufficiently large $n$. We remark that, to make the computation in Eq. (22) more tractable, one can further impose a stronger symmetry by restricting to the symmetric subspace of $(\mathbb{C}^{\mathfrak{p}+1})^{\otimes n}$ and imposing the positive partial transpose constraints Doherty et al. 2004; Tavakoli et al. 2024.

VI Application 1: Revisiting analytical results

To illustrate the physical significance of our approach, we apply it to the three-dimensional magnetic field estimation problem. Consider a spin-$1/2$ particle subjected to a magnetic field $\bm{B}=(B_{1},B_{2},B_{3})^{T}$. The dynamics of the system is governed by the Hamiltonian $H=\sum_{i=1}^{3}{\mu B_{i}\sigma_{i}/2}$, where $\sigma_{i}$, $i=1,2,3$, denote the Pauli matrices. The coefficient $\mu$ represents the proportionality coefficient relating the spin operator $\sigma_{i}/2$ to the magnetic moment. From a parameter estimation perspective, the above problem is mathematically equivalent to estimating $\bm{\theta}=\left(\theta_{1},\theta_{2},\theta_{3}\right)^{T}$ encoded in the Hamiltonian $H=\theta_{1}\sigma_{1}+\theta_{2}\sigma_{2}+\theta_{3}\sigma_{3}$, with $\theta_{i}=\mu B_{i}/2$. The associated signal channel is given by the unitary evolution

$$U_{\bm{\theta}}=e^{-iHt}$$ (23)

where $t$ is the evolution time. We here consider the noiseless case, leaving the discussion on the noisy scenario to Sec. VII. The purpose of this section is to employ our approach to reveal the strengths and limitations of some analytical results, which exist for parallel and sequential strategies.

The highest achievable precision is already well understood in the noiseless case for sequential strategies. An analytical expression for this precision is derived in Ref. Yuan 2016 and the corresponding bound is shown to be always attainable by a physically implementable sequential strategy. In contrast, analyzing parallel strategies is much more challenging. To address this challenge, a line of research has focused on proposing some heuristic probe states and using the precision attained by these states to benchmark the highest achievable precision in parallel strategies. For example, the authors of Ref. Baumgratz and Datta 2016 have proposed a class of permutationally invariant heuristic probe states with a two-body reduced density matrix given by

$$\rho^{\left[2\right]}=\frac{1}{4}\mathbb{I}_{2}\otimes\mathbb{I}_{2}+\frac{1}{12}\sum_{k=1}^{3}{\sigma_{k}\otimes\sigma_{k}}$$ (24)

and a one-body reduced density matrix $\rho^{\left[1\right]}=\mathbb{I}_{2}/2$. It has been shown that the estimation error for these states is

$$\mathrm{tr}\left(\Sigma\right)=\frac{3}{4N\left(N+2\right)}\left[\frac{1}{t^{2}}+\frac{2\left\|\bm{\theta}\right\|^{2}}{\sin^{2}\left(\left\|\bm{\theta}\right\|t\right)}\right].$$ (25)

Here $\left\|\bm{\theta}\right\|$ denotes the Euclidean norm of the vector $\bm{\theta}$, and the weight matrix $W$ is chosen to be the identity matrix. Subsequently, Ref. Hou et al. 2020 has improved upon this result by introducing a new class of heuristic probe states that achieve higher estimation precision. Nevertheless, it remains unclear whether this improved precision is the highest.

Figure 2: Estimation error as a function of $t$ for $\theta_{1}=\theta_{2}=1/2$ and $\theta_{3}=\sqrt{2}/2$ with $N=2$. The blue curve depicts the estimation error given by eq. 25. The yellow diamonds represent the estimation error associated with the heuristic states proposed in Ref. Hou et al. 2020. The red circles are the upper bounds on the optimal estimation error, obtained from our approach by randomly generating $m=125$ real unit vectors $\left\{|w_{x}\rangle\right\}_{x=1}^{m}$.

With the aid of our approach, we can answer this question. As an example, we consider the parameter configuration $\theta_{1}=\theta_{2}=1/2$ and $\theta_{3}=\sqrt{2}/2$ and compare the estimation error achieved by these heuristic probe states with the upper bound obtained from our approach. Hereafter, we use $B_{+}^{(k)}$ and $B_{-}^{(k)}$ to denote the upper and lower bounds obtained from our approach for the strategies of type $k$. For instance, $B_{+}^{(i)}$ and $B_{-}^{(i)}$ denote the upper and lower bounds for parallel strategies, respectively. The results are shown in fig. 2, where we plot the estimation error as a function of $t$ for $N=2$. Throughout, all numerical calculations are carried out with MOSEK ApS 2019. As can be seen from fig. 2, the estimation error achieved by all the heuristic probe states is larger than the upper bound $B_{+}^{(i)}$, obtained from our approach by randomly generating $m=125$ real unit vectors $\left\{|w_{x}\rangle\right\}_{x=1}^{m}$. This observation indicates that these heuristic probe states are not optimal, thereby revealing the limitations of these analytical results.

Figure 3: Illustration of the coincidence between the analytical bound in eq. 26 and the upper/lower bounds obtained from our approach. The yellow diamonds represent the precision associated with the heuristic states proposed in Ref. Hou et al. 2020. The blue line corresponds to the analytical bound in eq. 26. The red circles and green crosses denote the upper and lower bounds obtained from our approach with $m=125$ and $n=2$, respectively. Here, we set $\theta_{1}=\theta_{2}=1/2$ and $\theta_{3}=\sqrt{2}/2$ with $N=2$.

We note that another notable analytical result is the analytical bound derived in Ref. Hou et al. 2020. It states that the estimation error for parallel strategies is lower bounded by

$$\mathrm{tr}\left(\Sigma\right)\geq\frac{\left(1+2\frac{\left\|\bm{\theta}\right\|t}{\left|\sin\left(\left\|\bm{\theta}\right\|t\right)\right|}\right)^{2}}{4N\left(N+2\right)t^{2}}.$$ (26)

This bound is asymptotically attainable in the limit of $N\rightarrow\infty$, by employing the heuristic probe states introduced in Ref. Hou et al. 2020. However, for small $N$, these states fail to saturate the bound in general, as illustrated in Fig. 3. As a result, it remains unclear whether the analytical lower bound in eq. 26 is tight for finite and small values of $N$. We numerically demonstrate that this bound may indeed be tight even for small $N$. An explicit example is shown in fig. 3, where we compare the analytical lower bound in eq. 26 with the upper and lower bounds obtained from our approach for the parameter configuration $\theta_{1}=\theta_{2}=1/2$ and $\theta_{3}=\sqrt{2}/2$ with $N=2$. As depicted in fig. 3, the upper and lower bounds, $B_{+}^{(i)}$ and $B_{-}^{(i)}$, coincide with the analytical lower bound in eq. 26, indicating that the bound (26) is indeed tight for this specific case. To further support our observation, we have conducted extensive numerical calculations for various parameter configurations. Figure 4 shows the differences between the analytical lower bound (26) and the upper bound $B_{+}^{(i)}$ obtained from our method for different parameter configurations satisfying $\left\|\bm{\theta}\right\|=1$. As can be seen from fig. 4, these differences are nearly zero, irrespective of the specific choice of $\bm{\theta}$. We have also performed numerical calculations for other norms of $\bm{\theta}$ and observed similar behaviors. We therefore present only the results for $\left\|\bm{\theta}\right\|=1$ here. Based on these numerical results, we conjecture that the analytical lower bound in eq. 26 remains tight even for small $N$. The rigorous proof of this statement should be an interesting topic for future work.

Figure 4: Differences between the analytical lower bound in eq. 26 and the upper bound $B_{+}^{(i)}$ obtained from our method for different parameter configurations satisfying $\left\|\bm{\theta}\right\|=1$. Each point corresponds to a specific choice of $\theta_{1}$ and $\theta_{2}$, with $\theta_{3}$ determined through the relation $\theta_{3}=\sqrt{1-\theta_{1}^{2}-\theta_{2}^{2}}$. The color denotes the magnitude of the difference. The white region represents parameter configurations that do not satisfy the constraint $\left\|\bm{\theta}\right\|=1$. Here, we set $t=3,m=700$, and $N=2$.

VII Application 2: Establishing strict Hierarchy

Let us now take noise into account. Assume that the signal channel $\mathcal{E}_{\bm{\theta}}$ is described by the unitary evolution in Eq. (23) followed by an amplitude damping noise. This scenario is relevant in various physical platforms, such as nitrogen-vacancy centers in diamond Doherty et al. 2013 and superconducting qubits Kjaergaard et al. 2020, where the spin-$1/2$ particles are often subject to the amplitude damping noise. Explicitly, the Kraus operators of $\mathcal{E}_{\bm{\theta}}$ are $K_{1}U_{\bm{\theta}}$ and $K_{2}U_{\bm{\theta}}$, with $K_{1}$ and $K_{2}$ denoting the Kraus operators for the amplitude damping noise

$$K_{1}=\left[\begin{matrix}1&0\\ 0&\sqrt{1-\gamma}\\ \end{matrix}\right],\quad K_{2}=\left[\begin{matrix}0&\sqrt{\gamma}\\ 0&0\\ \end{matrix}\right].$$ (27)

Here $\gamma$ is the noise strength. It is worth noting that an outstanding topic in quantum metrology is to explore the hierarchy problem among different types of strategies, which has been extensively studied in the single-parameter estimation Giovannetti et al. 2006; Demkowicz-Dobrzański and Maccone 2014; Yuan 2016; Liu et al. 2023; Kurdziałek et al. 2023 but remains largely unexplored in the multiparameter scenario. The purpose of this section is to employ our approach to investigate this problem in the noisy multiparameter estimation scenario.

Figure 5: Lower bound $B_{-}^{(i)}$ for parallel strategies and upper bound $B_{+}^{(ii)}$ for sequential strategies as functions of the noise strength $\gamma$. Here $B_{-}^{(i)}$ and $B_{+}^{(ii)}$ are obtained from our approach by setting $n=2$ and $m=1500$, respectively. The parameters are chosen as $\theta_{1}=\theta_{2}=1/2$ and $\theta_{3}=\sqrt{2}/2$ with $t=0.1$ and $N=2$.

Throughout this section, we set the parameter configuration as $\theta_{1}=\theta_{2}=1/2$ and $\theta_{3}=\sqrt{2}/2$ with $t=0.1$ and $N=2$. With the aid of our approach, we first examine the hierarchy relation between parallel and sequential strategies. The numerical results are shown in fig. 5, where we plot the lower bound $B_{-}^{(i)}$ for parallel strategies and the upper bound $B_{+}^{(ii)}$ for sequential strategies as functions of $\gamma$. Here $B_{-}^{(i)}$ and $B_{+}^{(ii)}$ are obtained from our approach by setting $n=2$ and $m=1500$, respectively. As can be seen from fig. 5, $B_{-}^{(i)}$ is strictly larger than $B_{+}^{(ii)}$, irrespective of the value of $\gamma$. This implies that the highest achievable precision in sequential strategies always exceeds that in parallel strategies, indicating a strict advantage of sequential strategies over parallel strategies.

Figure 6: Lower bound $B_{-}^{(ii)}$ for sequential strategies and upper bound $B_{+}^{(iii)}$ for causal superposition strategies as functions of the noise strength $\gamma$. Here $B_{-}^{(ii)}$ and $B_{+}^{(iii)}$ are obtained from our approach by setting $n=2$ and $m=1500$, respectively. The parameters considered here are chosen to be identical to those in Fig. 5. The inset enlarges the gap between the two bounds.

We next examine the hierarchy relation between sequential strategies and causal superposition strategies. Figure 6 shows the lower bound $B_{-}^{(ii)}$ for sequential strategies and the upper bound $B_{+}^{(iii)}$ for causal superposition strategies as functions of $\gamma$. Here $B_{-}^{(ii)}$ and $B_{+}^{(iii)}$ are obtained from our approach by setting $n=2$ and $m=1500$, respectively. Compared with the previous case (see Fig. 5), the gap between $B_{-}^{(ii)}$ and $B_{+}^{(iii)}$ is much smaller, as can be seen from fig. 6. Nevertheless, there is still a strict advantage of causal superposition strategies over sequential strategies, as can be seen in the inset of fig. 6.

We finally investigate the hierarchy relation between causal superposition strategies and general indefinite-causal-order strategies. The numerical results are presented in fig. 7, where we plot the difference $B_{-}^{(iii)}-B_{+}^{(iv)}$ as a function of $\gamma$. Here $B_{-}^{(iii)}$ and $B_{+}^{(iv)}$ are obtained from our approach by setting $n=2$ and $m=1500$, respectively. Again, we observe from fig. 7 that $B_{-}^{(iii)}$ is always larger than $B_{+}^{(iv)}$, indicating a strict advantage of general indefinite-causal-order strategies over causal superposition strategies. Now, combing all the results in figs. 5, 6 and 7, we conclude that there exists a strict hierarchy among the four types of strategies in the noisy scenario.

Figure 7: The difference $B_{-}^{(iii)}-B_{+}^{(iv)}$ as a function of the noise strength $\gamma$. Here $B_{-}^{(iii)}$ and $B_{+}^{(iv)}$ are obtained from our approach by setting $n=2$ and $m=1500$, respectively. The parameters considered here are chosen to be identical to those in Fig. 5.

VIII Conclusion

We have presented a unified and computable approach for the simultaneous estimation of multiple parameters that attains the ultimate precision permitted by quantum mechanics. The key innovation here is to combine the quantum tester formalism and the TCRB, which together provide a unified framework for dealing with diverse quantum resources, including entanglement, coherence, quantum control, and indefinite causal order. Our approach is therefore versatile and can be applied to various types of strategies both in the noiseless and noisy scenarios.

We have clearly demonstrated the physical significance of our approach by applying it to the three-dimensional magnetic-field estimation. In the noiseless case, our numerical results identify strategies that outperform previously proposed heuristic ones and provide strong evidence that the analytical lower bound in Ref. Hou et al. 2020 is saturable even for finite and small $N$. In the noisy scenario, our approach enables a systematic comparison among different classes of strategies and establishes a strict hierarchy among parallel, sequential, causal-superposition, and general indefinite-causal-order strategies. These results extend the notion of strategy hierarchy, well understood in single-parameter estimation, to the genuinely multiparameter regime.

We clarify that our approach can accommodate arbitrary weight matrices, general noise models, and various estimation strategies. Moreover, we remark that our approach may, in principle, be extended to Bayesian estimation problems with prior distributions Demkowicz-Dobrzański et al. 2020; Mukhopadhyay et al. 2025, providing a unified tool for both local and Bayesian multiparameter quantum metrology. We expect that the results presented in this work could facilitate a deeper understanding of ultimate achievable precision limits in multiparameter estimation and serve as a versatile benchmark for assessing and designing optimal quantum metrological strategies.

DATA AVAILABILITY

The code used in this article is openly available from cod.

Appendix A Expressions for $\Lambda^{(k)}$

Here we present the explicit expressions for the maps $\Lambda^{(k)}$, $k=i,ii,iv$, introduced in eq. 5. In what follows, $\prescript{}{1-Q}{X^{(k)}}\coloneqq X^{(k)}-\prescript{}{Q}{X^{(k)}}$, where $\prescript{}{Q}{X^{(k)}}$ denotes the operation that traces out subsystem $Q$ and replaces it with the normalized identity operator, i.e., $\prescript{}{Q}{X^{(k)}}\coloneqq\mathrm{tr}_{Q}X^{(k)}\otimes\left(\mathbb{I}_{Q}/d_{Q}\right)$ where $d_{Q}=\mathrm{dim}(Q)$. The expression for $\Lambda^{(i)}$ is

$$\Lambda^{(i)}(X^{(i)})\coloneqq\prescript{}{O_{1}\cdots O_{N}}{X^{(i)}}.$$ (28)

The expression for $\Lambda^{(ii)}$ is

$$\Lambda^{(ii)}(X^{(ii)})\coloneqq\prescript{}{O_{N}}{X^{(ii)}}-\prescript{}{(1-O_{N-1})I_{N}O_{N}}{X^{(ii)}}-$$ (29)

$$\prescript{}{(1-O_{N-2})I_{N-1}O_{N-1}I_{N}O_{N}}{X^{(ii)}}-\cdots-\prescript{}{(1-O_{1})I_{2}O_{2}\cdots I_{N}O_{N}}{X^{(ii)}}.$$

The expression for $\Lambda^{(iv)}$ is

$$\Lambda^{(iv)}(X^{(iv)})\coloneqq\prescript{}{\left[1-\prod_{j=1}^{N}\left(1-O_{j}+I_{j}O_{j}\right)+\prod_{j=1}^{N}I_{j}O_{j}\right]}{X^{(iv)}}.$$ (31)

Appendix B Strategies of type $iii$

To simplify the notations, we focus on the case of $N=2$. All results presented here can be extended straightforwardly to arbitrary $N$. Denote $1\prec 2$ as the causal order such that the first channel is queried before the second channel, and similarly for $2\prec 1$. It has been shown Araújo et al. 2015; Bavaresco et al. 2021; Liu et al. 2023; Zhou et al. 2024 that the set $\mathcal{X}^{(iii)}$ can be characterized as

$$\mathcal{X}^{\left(iii\right)}\coloneqq\left\{X|X=pX^{\left(1\prec 2\right)}+\left(1-p\right)X^{\left(2\prec 1\right)},0\leq p\leq 1\right\},$$ (32)

where $X^{\left(1\prec 2\right)}\in\mathcal{X}^{(ii)}$ has causal order $1\prec 2$ and $X^{\left(2\prec 1\right)}\in\mathcal{X}^{(ii)}$ has causal order $2\prec 1$. With the characterization of $\mathcal{X}^{(iii)}$ given in eq. 32, the upper bounds for causal superposition strategies can be written explicitly as

$$\underset{\begin{array}[]{c}X_{x}^{\left(iii\right)}\\ \tilde{X}^{\left(1\prec 2\right)},\tilde{X}^{\left(2\prec 1\right)}\\ \end{array}}{\min}\,\, \sum_{x}{\langle w_{x}|\tilde{W}|w_{x}\rangle\mathrm{tr}\left(C_{\bm{\theta}}X_{x}^{\left(iii\right)}\right)}$$

$$\mathrm{s}.\mathrm{t}.\,\, \sum_{x}{\left|\langle w_{x}|0\rangle\right|^{2}X_{x}^{\left(iii\right)}}=\tilde{X}^{\left(1\prec 2\right)}+\tilde{X}^{\left(2\prec 1\right)},$$ (33c)

$$\Lambda^{\left(1\prec 2\right)}\left(\tilde{X}^{\left(1\prec 2\right)}\right)=\tilde{X}^{\left(1\prec 2\right)},$$ (33d)

$$\Lambda^{\left(2\prec 1\right)}\left(\tilde{X}^{\left(2\prec 1\right)}\right)=\tilde{X}^{\left(2\prec 1\right)},$$ (33e)

$$\mathrm{tr}\left[\tilde{X}^{\left(1\prec 2\right)}+\tilde{X}^{\left(2\prec 1\right)}\right]=d_{O},$$ (33f)

$$\tilde{X}^{\left(1\prec 2\right)}\geq 0,\quad\tilde{X}^{\left(2\prec 1\right)}\geq 0,$$ (33g)

$$X_{x}^{\left(iii\right)}\geq 0,\quad\forall x,$$ (33h)

$$\sum_{x}{\langle w_{x}|A_{i}|w_{x}\rangle\mathrm{tr}\left(\partial_{j}C_{\bm{\theta}}X_{x}^{\left(iii\right)}\right)}=\delta_{ij},$$

$$\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad i,j=1,\cdots,\mathfrak{p},$$ (33i)

where the variable $p$ in eq. 32 has been absorbed into $\tilde{X}^{\left(1\prec 2\right)}$ and $\tilde{X}^{\left(2\prec 1\right)}$, i.e., $\tilde{X}^{\left(1\prec 2\right)}=pX^{\left(1\prec 2\right)}$ and $\tilde{X}^{\left(2\prec 1\right)}=(1-p)X^{\left(2\prec 1\right)}$. Inserting eq. 32 into eq. 22, we can express the lower bounds for the causal superposition strategies as

$$\underset{\begin{array}[]{c}Y^{(iii)}_{n}\\ \tilde{X}^{\left(1\prec 2\right)},\tilde{X}^{\left(2\prec 1\right)}\\ \end{array}}{\min}\,\, \mathrm{tr}\left[\tilde{W}\otimes C_{\bm{\theta}}Y^{(iii)}\right]$$

$$\mathrm{s}.\mathrm{t}.\,\, \mathrm{tr}_{1}\left[|0\rangle\langle 0|\otimes\mathbb{I}_{IO}Y^{(iii)}\right]=\tilde{X}^{\left(1\prec 2\right)}+\tilde{X}^{\left(2\prec 1\right)},$$ (34c)

$$\text{Eqs.~(\ref{eq_refbegin}-\ref{eq_refend}) hold},$$ (34d)

Eqs. (33d-33g) hold. (34e)