# A unified toolbox for multipartite entanglement certification

Ye-Chao Liu, Jannis Halbey, Sebastian Pokutta, Sébastien Designolle

► To cite this version:

Ye-Chao Liu, Jannis Halbey, Sebastian Pokutta, Sébastien Designolle. A unified toolbox for multipartite entanglement certification. 2025. hal-05192936

HAL Id: hal-05192936

https://hal.science/hal-05192936v1

Preprint submitted on 30 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Unified Toolbox for Multipartite Entanglement Certification

Ye-Chao Liu,¹ Jannis Halbey,¹ Sebastian Pokutta,¹ and Sébastien Designolle¹,²

¹Zuse Institute Berlin, Takustraße 7, 14195 Berlin, Germany

²Inria, ENS de Lyon, UCBL, LIP, 69342, Lyon Cedex 07, France

(Dated: 23rd July 2025)

We present a unified framework for multipartite entanglement characterization based on the conditional gradient (CG) method, incorporating both fast heuristic detection and rigorous witness construction with numerical error control. Our method enables entanglement certification in quantum systems of up to ten qubits and applies to arbitrary entanglement structures. We demonstrate its power by closing the gap between entanglement and separability bounds in white noise robustness benchmarks for a class of bound entangled states. Furthermore, the framework extends to entanglement robustness under general quantum noise channels, providing accurate thresholds in cases beyond the reach of previous algorithmic methods. These results position CG methods as a powerful tool for practical and scalable entanglement analysis in realistic experimental settings.

*Introduction.*—Entanglement [1] is a fundamental feature of quantum mechanics that distinguishes it from classical physics, and serves as a central resource in quantum technologies such as communication, computation, metrology, and teleportation. Accordingly, the detection and characterization of entanglement remain pivotal tasks in quantum information science [2]. For bipartite systems, significant progress has been made through various entanglement criteria. The Peres-Horodecki criterion, also known as the positive partial transpose (PPT) criterion, is a well-known example [3, 4]. Another one is the computable cross norm or realignment (CCNR) criterion, which is more experimentally accessible and can detect certain bound entangled states that the PPT criterion cannot [5, 6].

When it comes to multipartite systems, the challenges become significantly greater. Not only does the Hilbert space grow exponentially, but the multipartite entanglement also has much richer structures, ranging from fully separable to genuine multipartite entanglement (GME). While both PPT and CCNR criteria have been extended to detect GME [7], entanglement witnesses remain the most widely used tool for multipartite systems [8–11]. A entanglement witness $W$ is a necessary and sufficient criteria satisfying $\text{tr}(W\rho) \ge 0$ for all separable $\rho$, and $\text{tr}(W\rho) < 0$ for at least one entangled $\rho$. They can be analytically constructed for arbitrary pure states [12] but often fall short when dealing with general (mixed) entangled states, including bound entangled states.

To overcome these limitations, algorithmic methods have been developed to address arbitrary entangled states. The PPT criterion can be formulated as a semidefinite program (SDP), and its symmetric extensions, known as the Doherty-Parrilo-Spedalieri (DPS) hierarchy, can asymptotically approximate the separable space [13, 14]. On the other hand, by inner polytope approximation of the Hilbert space, there are also methods based on SDP that can certify separability [15]. State-of-the-art SDP-based methods have been demonstrated for systems up to six qubits, yet they remain inadequate for current experimental regimes involving more than ten qubits [16], where computational breakdowns often occur.

<figure><img src="image_2.png" /><figcaption>Figure 1. Schematic overview of our unified toolbox for multipartite entanglement certification based on conditional gradient (CG) methods. Starting from a target state $\rho$, the heuristic CG iteration evaluates the dual gap efficiency $r$ to determine the appropriate certification path. If $r \le r_0$, entanglement is rapidly detected; otherwise, further CG iterations are performed or a rigorous entanglement witness is constructed to certify entanglement while accounting for numerical error. Alternatively, geometric reconstruction from the same CG-based output can be used to certify separability.</figcaption></figure>

Beyond SDP-based methods, conditional gradient (CG) methods have emerged as a promising alternative for multipartite entanglement characterization. The idea was first explored by Kampermann et al. [17], and later developed into suitable algorithms within CG methods [18–21]. A key advantage of CG methods is their flexibility in handling arbitrary multipartite entanglement structures. However, existing CG-based techniques have primarily focused on separability certification and are constrained to small systems as well.

In this work, we introduce a CG-based framework capable of certifying entanglement for multipartite quantum systems. Our approach significantly extends the scalability of CG methods, enabling rapid entanglement detection in systems of up to ten qubits. Moreover, it yields rigorous entanglement witnesses suitable for experimental validation. Combined with separability certification via the same CG-based framework, our method constitutes a unified toolbox for multipartite entanglement analysis. We demonstrate the power of our method by providing optimal entanglement thresholds on the entanglement robustness for a class of bound entangled states, thereby
closing the gap between entanglement and separability thresholds. Last but not least, our toolbox can assess robustness under nonlinear noise models and realistic noise channels, making it broadly applicable in experimental scenarios. As illustrated in Fig. 1, our unified toolbox integrates fast detection, rigorous certification, and separability analysis into a coherent structure.

**Entanglement and Frank-Wolfe algorithm.**—We begin by recalling the definition of multipartite entanglement. A pure state is called k-separable if it can be written as $|\psi^{ks}\rangle = |\phi_1\rangle \otimes |\phi_2\rangle \otimes \cdots \otimes |\phi_t\rangle$, where each $|\phi_i\rangle$ may itself be entangled within a subset of the total system. For a mixed state, it is k-separable if it can be expressed as a convex combination of k-separable pure states: $\sigma^{ks} = \sum_i \lambda_i |\psi_i^{ks}\rangle\langle\psi_i^{ks}|$, where the coefficients $\lambda_i$ form a probability distribution. Usually, if n-partite quantum states are not n-separable, then they are called entangled states. And non-2-separable states are called GME.

Accordingly, the set of all k-separable states forms a convex subset $S_k$ of the state space, and any state outside $S_k$ is entangled. Characterizing multipartite entanglement then amounts to determining whether a given quantum state $\rho$ belongs to $S_k$. This naturally leads to the following optimization problem:

$$f(\rho) = \min_{\sigma \in S_k} f(\rho, \sigma), \qquad (1)$$

where $f(\rho, \sigma)$ denotes a distance measure between quantum states. If $f(\rho) = 0$, the state is k-separable; otherwise, it is entangled. A common choice for the distance measure is the distance induced by the Hilbert-Schmidt norm $\|\rho - \sigma\|$, also known as the Frobenius norm or $l_2$ norm. For algorithmic convenience, we consider half of its squared form $f(\rho, \sigma) = \frac{1}{2}\|\rho - \sigma\|^2$, which is convex and differentiable. Although not a strict distance, we will refer to $f(\rho)$ as a “distance” throughout, reflecting its geometric role in the optimization.

This is a constrained convex optimization problem, and the typical algorithmic tool is the Frank-Wolfe method, also known as the conditional gradient method. The basic procedure for entanglement characterization via the vanilla Frank-Wolfe algorithm is as follows:

### Vanilla Frank-Wolfe algorithm

**Input:** Distance function $f(x)$; initial point $\sigma_0 \in S_k$

1: **for** $t = 0, \dots, m$ **do**

2: &nbsp;&nbsp;&nbsp;&nbsp; $\psi_t \leftarrow \arg \min_{\psi \in S_k} \langle \psi | \sigma_t - \rho | \psi \rangle$

3: &nbsp;&nbsp;&nbsp;&nbsp; $\gamma_t \leftarrow \arg \min_{\gamma} f(\sigma_t + \gamma(\psi_t - \sigma_t))$

4: &nbsp;&nbsp;&nbsp;&nbsp; $\sigma_{t+1} \leftarrow \sigma_t + \gamma_t(\psi_t - \sigma_t)$

5: **end for**

**Output:** $\sigma_1, \dots, \sigma_m \in S_k$

Each iteration constructs a new state $\sigma_{t+1}$ such that $f(\rho, \sigma_{t+1}) \le f(\rho, \sigma_t)$. The product state $\psi_t$, obtained via a linear minimization oracle (LMO), corresponds to the extreme point in $S_k$ that is most aligned with the descent direction $\sigma_t - \rho$, i.e., it minimizes the directional derivative of the distance function along this direction.

This makes the update efficient and steers $\sigma_t$ progressively toward $\rho$. An essential advantage of Frank-Wolfe algorithms over other gradient-based approaches is to avoid costly projections onto $S_k$, which are as difficult as the problem (1) we are trying to solve, which is itself known to be NP-hard [22].

A related idea was previously introduced in the work of Kampermann et al. [17], and later studies have adopted variants based on Gilbert's algorithm [18–21]. Although Gilbert's algorithm [23] was used earlier in the quantum information literature, it is a special case of the more general Frank-Wolfe method [24, 25], which was introduced earlier and has since seen significant advances in algorithmic theory. In this work, we adopt the term conditional gradient (CG) method to unify terminology, encompassing both the vanilla Frank-Wolfe (or Gilbert's) algorithm and its modern improvements. Further details and algorithmic enhancements used in this work can be found in Appendix A and in recent review articles [26, 27]. Our code is available as a Julia package [28].

**Rapid entanglement detection.**—By avoiding explicit projections onto the separable space, the main computational cost of our method lies in evaluating the linear minimization oracle (LMO), which can be written as

$$\min_{|\phi_i\rangle \in \mathcal{H}_i, i=1, \dots, k} \langle \phi_1 | \otimes \cdots \otimes \langle \phi_k | (\sigma_t - \rho) | \phi_1 \rangle \otimes \cdots \otimes | \phi_k \rangle . \quad (2)$$

This optimization can be heuristically solved via alternating updates over subsystems. However, due to numerical error and the inherent limitations of finite iterations, the computed quantity $f(\rho, \sigma_t)$ generally does not converge to zero, even for separable states. Thus, $f(\rho, \sigma_t)$ cannot be directly used to certify entanglement or separability. In previous works [18], a geometric reconstruction approach was developed within the CG-based framework for separability certification.

Here, we take the complementary perspective to utilize CG methods for certifying entanglement. Our approach is based on the following proposition:

**Proposition 1.** A quantum state $\rho$ is entangled if and only if there exists a separable state $\sigma$ such that $\|\sigma - \tau\| < \|\rho - \tau\|$ for all $\tau \in S_k$.

*Proof.* Necessity.— Let $\sigma = \arg \min_{\sigma \in S_k} \|\rho - \sigma\|$. For any $\tau \in S_k$, we have $\text{tr}[(\sigma - \rho)(\sigma - \tau)] < 0$, and therefore we achieve $\|\rho - \tau\|^2 > \|\rho - \sigma\|^2$ following the parallelogram identity. Sufficiency.— If $\rho$ is separable, let $\tau = \rho \in S_k$. Then the assumption $\|\sigma - \rho\| < \|\rho - \rho\| = 0$ leads to a contradiction. □

We now introduce the dual gap $g_t = \text{tr}[(\sigma_t - \rho)(\sigma_t - \psi_t)]$, a standard convergence indicator in CG methods [27]. With the parallelogram identity, we can rewrite it as

$$g_t = \frac{1}{2} \|\rho - \sigma_t\|^2 + \frac{1}{2} \|\sigma_t - \psi_t\|^2 - \frac{1}{2} \|\rho - \psi_t\|^2 . \quad (3)$$

This motivates the following entanglement criterion:

**Corollary 2.** A quantum state $\rho$ is entangled if and only if the inequality $f(\rho, \sigma_t) > g_t$ holds.
<figure><img src="image_3.png" /><figcaption>Figure 2. The ratio $r_t = g_t/f(\rho, \sigma_t)$ for (a) three-qubit GHZ state mixed with varying levels of white noise, and (b) ten-qubit GHZ and Dicke states with 70% white noise. The dashed line denotes the threshold $r = 1$ for entanglement detection.</figcaption></figure>

Therefore, instead of simply checking whether $f(\rho, \sigma_t) > 0$, we establish a practical entanglement detection criterion by defining the dual gap efficiency

$$r_t := \frac{g_t}{f(\rho, \sigma_t)}, \qquad (4)$$

which compares the dual gap $g_t$ with the current distance $f(\rho, \sigma_t)$. The algorithm terminates and declares $\rho$ entangled once $r_t < 1$. As shown in Fig. 2(a), the three-qubit GHZ state with 80% white noise lies exactly on the boundary of $S_3$, resulting in a nearly constant $r_t$ throughout iterations. When the noise is reduced by 1%, introducing slight entanglement, the criterion is satisfied after a few thousand iterations. States with higher entanglement (i.e., less noise) are certified more rapidly.

This capability enables entanglement detection in large systems, particularly in practical scenarios. In experiments, entangled states often contain moderate noise, yet larger systems typically tolerate more. To illustrate this, Fig. 2(b) shows that 10-qubit GHZ and Dicke states mixed with 70% white noise can still be certified within dozens of iterations. The total runtime is under one minute using a 2.80 GHz CPU and 16 GB of memory. In the main text, we focus on detecting full separability or bipartite entanglement (i.e., $S_k$ with $k = n$ for n-qubit systems); for general entanglement structures, see Appendix B.

Our results reveal a clear connection between the amount of entanglement and the computational effort required. For a fixed system size, more entangled states require fewer iterations to satisfy the criterion. For instance, the three-qubit GHZ state with less noise converges faster; the Dicke state outperforms the GHZ state in the 10-qubit case, as it exhibits stronger robustness to white noise.

However, this property also implies that detecting weak entanglement near the boundary of $S_k$ may require prohibitively many iterations. Moreover, since the LMO problem (2) is solved heuristically, there is no guarantee of global optimality. This compromises the reliability of Corollary 2, which critically depends on the accuracy of $\psi_t$ and $g_t$. In practice, this can be mitigated by tightening the threshold, e.g., accepting $r_t < 1/5$, which suffices in most cases, although without theoretical guarantee.

In order to strictly avoid the above two possible problems, we introduce in the following an entanglement witness, which explicitly accounts for numerical uncertainty and therefore can certify entanglement rigorously. It can be seen as a subroutine within the same algorithmic framework.

**Rigorous entanglement witness.**—An entanglement witness for a generic quantum state can be constructed as [29–31]

$$W = \sigma - \rho_e + \mathrm{Tr}[\sigma(\rho_e - \sigma)] \cdot \mathbf{1}, \qquad (5)$$

where $\rho_e$ is the target entangled state to be certified, and $\sigma$ is the closest separable state in $S_k$ (normalization omitted). The main difficulty lies in finding $\sigma$, which corresponds exactly to the solution of Eq. (1).

Instead of the heuristic LMO procedure, we now consider a structured search over an $\epsilon$-net $S_k^\epsilon$ of the boundary of the separable set $S_k$. It is a finite set of pure separable states satisfying the condition that for every $\sigma$ in the boundary of $S_k$, there exists $\sigma^\epsilon \in S_k^\epsilon$ such that $\|\sigma - \sigma^\epsilon\| \le \epsilon$. A general construction of $S_k^\epsilon$ is given in Appendix C.

With this discretization, the original optimization problem can be replaced by a finite search, allowing for rigorous error control. In particular, the deviation from the true optimum is bounded as

$$0 \le \|\rho - \sigma^\epsilon\| - \|\rho - \sigma^*\| \le \epsilon, \qquad (6)$$

where $\sigma^*$ denotes the optimal solution to Eq. (1), and $\sigma^\epsilon$ is the solution obtained from $S_k^\epsilon$.

This enables us to construct a rigorous entanglement witness as follows:

**Proposition 3.** Given an entangled state $\rho$ and a separable state $\sigma$, let the product state $\phi$ in $S_k^\epsilon$ be the closest to the direction $\Lambda = \sigma - \rho$. Then the operator

$$W = \frac{\Lambda - (\beta - \epsilon) \cdot \mathbf{1}}{\|\Lambda\|} \qquad (7)$$

defines an entanglement witness for the quantum state $\rho$, where $\beta = \mathrm{tr}(\Lambda|\phi\rangle\langle\phi|)$ and $\epsilon = (1 - \eta)\|\Lambda\|$. Here, the factor $\eta$ depends on $\epsilon$ and the construction of $S_k^\epsilon$.

*Proof.* See the proof in Appendix D. $\square$

While the search over the $\epsilon$-net space can be computationally expensive, it is only required as a subroutine within the CG-based framework and not in every iteration. Its effectiveness depends on the quality of the reference separable state $\sigma$, which serves as an anchor for constructing the witness.

In practice, this structured search complements the fast heuristic method. The heuristic CG procedure is used initially to rapidly detect entanglement in most practical scenarios. If the ratio $r_t$ fails to converge to the detection threshold, or when experimental certification or precise boundary identification is needed, the rigorous entanglement witness can be invoked to ensure validity.
<figure><img src="image_4.png" /><figcaption>Figure 3. (a) Entanglement bounds $p_{ent,CG}$ (blue line) on the white noise robustness of Horodecki states $\rho^H(a)$, obtained via our CG method, compared with previous entanglement bounds based on a generalization of the Wootters formula $p_{ent,GW}$ (blue dots) and separability bounds based on SDP methods $p_{sep,SDP}$ (orange dots). (b) Comparison of the gap $\Delta p$ between entanglement and separability bounds: blue shows the range $[p_{ent,GW}, p_{sep,SDP}]$, green shows $[p_{ent,CG}, p_{sep,SDP}]$, and orange shows $[p_{ent,CG}, p_{sep,CG}]$.</figcaption></figure>

**Applications for closing gap.**—To demonstrate the power of our method, particularly its capacity for rigorous entanglement certification, we consider the problem of estimating the white noise robustness of entangled states. This task serves as a standard and widely used benchmark for evaluating entanglement and separability detection tools. It is to determine the maximal white noise level $p$ such that the mixed state $\rho(p) = (1-p)\rho + p \mathbb{1}/d$ remains separable; that is, $p_{sep} = \max\{p \ge 0 : \rho(p) \in S_k\}$.

As previously mentioned, the CG-based framework can also certify separability using a geometric reconstruction procedure [18]. In our framework, it can be directly achieved based on the same algorithmic results of Eq. (1) as the following Proposition 4. This allows us to construct a unified toolbox for closing the gap between entanglement and separability bounds in the white noise robustness problem.

**Proposition 4.** *Given a noisy (separable or entangled) quantum state $\rho(p) = (1-p)\rho + p \mathbb{1}/d$, and suppose $\|\rho(p) - \sigma\| = \delta$ for some $\sigma \in S_k$. Then the state $\rho(p+\frac{\epsilon}{1+\epsilon})$ must be separable, where $\epsilon \ge \delta/a$, and $a$ is the radius of the separable ball.*

*Proof.* See the proof in Appendix E. □

As a concrete example, we apply our method to the Horodecki states $\rho^H(a)$, a class of 3 × 3 bound entangled states introduced by P. Horodecki [4]. These states are entangled for all $0 < a < 1$, yet cannot be detected by the PPT criterion. The explicit form of $\rho^H(a)$ is provided in Appendix F 1.

As shown in Fig. 3(a), our CG-based approach yields significantly improved entanglement bounds $p_{ent,CG}$, outperforming the best known bounds from the generalized Wootters formula [32]. Notably, the new entanglement thresholds approach the previously established separability limits $p_{sep,SDP}$, nearly closing the gap. Fig. 3(b) further quantifies the improvement. Our method reduces the

<figure><img src="image_5.png" /><figcaption>Figure 4. Entanglement robustness thresholds for the Bell state $(|00\rangle + |11\rangle)/\sqrt{2}$ under various quantum noise channels: global depolarizing (white noise, GD), bit flip (BF), phase flip (PF), amplitude damping (AD), and phase damping (PD). Except for white noise (applied to both qubits), all channels act on one side of the state. Bars indicate theoretical thresholds for the noise strength (blue) and the corresponding fidelity with the ideal Bell state (orange). Error bars represent the certified interval between the entanglement and separability bounds obtained by our numerical algorithm.</figcaption></figure>

entanglement-separability gap by approximately an order of magnitude. Moreover, the entanglement and separability thresholds obtained from our CG-based approach yield a gap that is comparable to the tightest known gap between entanglement and separability bounds, defined by our results and SDP methods. We also note that in the parameter regime $0 < a < 0.02$, previous results exhibited inconsistencies between the entanglement and separability bounds. Our method not only narrows the gap but also resolves this contradiction, improving both numerical precision and physical consistency.

**Applications under arbitrary noise channels.**—From an algorithmic perspective, estimating the white noise robustness of a quantum state $\rho$ is computationally equivalent to certifying entanglement of its noisy counterpart $\rho(p)$. In SDP-based approaches, this typically requires only an additional variable in the optimization. In our CG-based framework, the parameter $p$ can be updated dynamically during iterations based on the ratio $r$.

Both approaches can generalize to arbitrary linear noise models of the form $(1-p)\rho + p\sigma$, $\sigma \in \mathcal{H}$. However, in realistic settings, quantum noise is often nonlinear and must be modeled as a quantum channel $\Phi$, defined via Kraus operators as $\Phi(\rho) = \sum_i K_i \rho K_i^\dagger$ with $\sum_i K_i^\dagger K_i \le \mathbb{1}$. For instance, the bit-flip channel is defined by $K_0 = \sqrt{1-p} \mathbb{1}$ and $K_1 = \sqrt{p} X$.

Different noise models yield different entanglement thresholds, which are of practical importance for experimental applications, e.g., determining fidelity thresholds for maintaining functionality in quantum protocols. The robustness of entanglement under general quantum channels lies beyond the reach of current SDP methods, but can be naturally addressed within our CG-based framework using the same convergence criterion based on $r_t$.

Fig. 4 presents the entanglement robustness thresholds for the Bell state $(|00\rangle + |11\rangle)/\sqrt{2}$ under several common
quantum noise channels relevant to quantum computation and teleportation. In each case, our numerical method yields certified entanglement and separability bounds with tight intervals, represented by error bars. The results show excellent agreement with known theoretical thresholds, demonstrating the accuracy and reliability of our method. Further details and comparisons with analytical values are provided in Appendix F 2.

**Conclusion.**—We have developed a unified framework for multipartite entanglement characterization based on the conditional gradient (CG) method, encompassing both heuristic and rigorous procedures. Our approach enables rapid entanglement detection for systems of up to ten qubits and constructs rigorous entanglement witnesses with explicit numerical error control. By applying this framework to the problem of white noise robustness for bound entangled states, we obtained the optimal entanglement thresholds, thereby closing the gap between entanglement and separability thresholds. We further demonstrated the applicability of our method in analyzing entanglement robustness under general quantum noise channels, providing accurate thresholds in regimes beyond the reach of SDP-based techniques. These results establish CG methods as a versatile and scalable tool for entanglement analysis in current experimental scenarios.

Beyond these results, our framework suggests several

promising directions for future work. Exploiting symmetries in quantum states and measurements may simplify optimization and improve interpretability, especially for highly symmetric states such as GHZ or Dicke states. The method can be extended to compute robustness measures by optimizing over both the separable space and the noisy space. Incorporating matrix product state (MPS) representations may scale the method to hundreds of qubits. Additionally, splitting methods can be integrated to address incomplete or partial information, broadening the applicability of our framework to a wider range of quantum characterization settings.

## ACKNOWLEDGMENTS

We are grateful to Márton Naszodi, H. Chau Nguyen, Ties-A. Ohst, Jiangwei Shang, Tamás Vértesi, and Liding Xu for helpful discussions. This work was supported by the DFG Cluster of Excellence MATH+ (EXC-2046/1, Project No. 390685689) funded by the Deutsche Forschungsgemeinschaft (DFG), and partly funded within the QuantERA II Programme that has received funding from the European Union’s Horizon 2020 research and innovation programme under Grant Agreement No 101017733 (VERIqTAS).

[1] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).

[2] O. Gühne and G. Tóth, Entanglement detection, Phys. Rep. 474, 1 (2009).

[3] A. Peres, Separability criterion for density matrices, Phys. Rev. Lett. 77, 1413 (1996).

[4] P. Horodecki, Separability criterion and inseparable mixed states with positive partial transposition, Physics Letters A 232, 333 (1997).

[5] K. Chen and L.-A. Wu, A matrix realignment method for recognizing entanglement, Quantum Info. Comput. 3, 193 (2003).

[6] O. Rudolph, Further results on the cross norm criterion for separability, Quantum Inf. Proc. 4, 219 (2005).

[7] L. Novo, T. Moroder, and O. Gühne, Genuine multipartite entanglement of permutationally invariant states, Phys. Rev. A 88, 012305 (2013).

[8] M. Horodecki, P. Horodecki, and R. Horodecki, Separability of mixed states: necessary and sufficient conditions, Phys. Lett. A 223, 1 (1996).

[9] B. M. Terhal, Bell inequalities and the separability criterion, Physics Letters A 271, 319 (2000).

[10] M. Lewenstein, B. Kraus, J. I. Cirac, and P. Horodecki, Optimization of entanglement witnesses, Phys. Rev. A 62, 052310 (2000).

[11] D. Bruß, J. I. Cirac, P. Horodecki, F. Hulpke, B. Kraus, M. Lewenstein, and A. Sanpera, Reflections upon separability and distillability, J. Mod. Opt. 49, 1399 (2002).

[12] M. Bourennane, M. Eibl, C. Kurtsiefer, S. Gaertner, H. Weinfurter, O. Gühne, P. Hyllus, D. Bruß, M. Lewenstein, and A. Sanpera, Experimental detection of multipartite entanglement using witness operators, Phys. Rev. Lett. 92, 087902 (2004).

[13] A. C. Doherty, P. A. Parrilo, and F. M. Spedalieri, Complete family of separability criteria, Phys. Rev. A 69, 022308 (2004).

[14] M. Navascués, M. Owari, and M. B. Plenio, Complete criterion for separability detection, Phys. Rev. Lett. 103, 160404 (2009).

[15] T.-A. Ohst, X.-D. Yu, O. Gühne, and H. C. Nguyen, Certifying quantum separability with adaptive polytopes, SciPost Phys. 16, 063 (2024).

[16] C. Song, K. Xu, W. Liu, C.-p. Yang, S.-B. Zheng, H. Deng, Q. Xie, K. Huang, Q. Guo, L. Zhang, P. Zhang, D. Xu, D. Zheng, X. Zhu, H. Wang, Y.-A. Chen, C.-Y. Lu, S. Han, and J.-W. Pan, 10-qubit entanglement and parallel logic operations with a superconducting circuit, Phys. Rev. Lett. 119, 180511 (2017).

[17] H. Kampermann, O. Gühne, C. Wilmott, and D. Bruß, Algorithm for characterizing stochastic local operations and classical communication classes of multiparticle entanglement, Phys. Rev. A 86, 032307 (2012).

[18] J. Shang and O. Gühne, Convex optimization over classes of multiparticle entanglement, Phys. Rev. Lett. 120, 050506 (2018).

[19] M. Wieśniak, P. Pandya, O. Sakarya, and B. Woloncewicz, Distance between bound entangled states from unextendible product bases and separable states, Quantum Rep. 2, 49 (2020).

[20] P. Pandya, O. Sakarya, and M. Wieśniak, Hilbert-Schmidt distance and entanglement witnessing, Phys. Rev. A 102,
012409 (2020).

[21] Y. Hu, Y.-C. Liu, and J. Shang, Algorithm for evaluating distance-based entanglement measures, Chinese Phys. B 32, 080307 (2023).

[22] L. Gurvits, Classical deterministic complexity of Edmonds' problem and quantum entanglement, in Proc. STOC'03, STOC '03 (Association for Computing Machinery, New York, NY, USA, 2003) pp. 10–19.

[23] E. G. Gilbert, An iterative procedure for computing the minimum of a quadratic form on a convex set, J. SIAM Control 4, 61 (1966).

[24] M. Frank and P. Wolfe, An algorithm for quadratic programming, Nav. Res. Logist. Q. 3, 95 (1956).

[25] E. S. Levitin and B. T. Polyak, Constrained minimization methods, USSR Comput. Math. & Math. Phys. 6, 1 (1966).

[26] I. M. Bomze, F. Rinaldi, and D. Zeffiro, Frank-Wolfe and friends: a journey into projection-free first-order optimization methods, 4OR 19, 313 (2021).

[27] G. Braun, A. Carderera, C. W. Combettes, H. Hassani, A. Karbasi, A. Mokhtari, and S. Pokutta, Conditional gradient methods, arXiv:2211.14103 (2022).

[28] See our code as a Julia package at https://github.com/ZIB-IOL/EntanglementDetection.jl.

[29] A. O. Pittenger and M. H. Rubin, Convexity and the separability problem of quantum mechanical density matrices, Linear Algebra Its Appl. 346, 47 (2001).

[30] R. A. Bertlmann, H. Narnhofer, and W. Thirring, Geometric picture of entanglement and Bell inequalities, Phys. Rev. A 66, 032319 (2002).

[31] R. A. Bertlmann, K. Durstberger, B. C. Hiesmayr, and P. Krammer, Optimal entanglement witnesses for qubits and qutrits, Phys. Rev. A 72, 052331 (2005).

[32] Z.-H. Chen, Z.-H. Ma, O. Gühne, and S. Severini, Estimating entanglement monotones with a generalization of the Wootters formula, Phys. Rev. Lett. 109, 200503 (2012).

[33] J. Guélat and P. Marcotte, Some comments on Wolfe's 'away step', Math. Program. 35, 110 (1986).

[34] C. A. Holloway, An extension of the Frank and Wolfe method of feasible directions, Math. Program. 6, 14–27 (1974).

[35] K. K. Tsuji, K. Tanaka, and S. Pokutta, Pairwise conditional gradients without swap steps and sparser kernel herding, in International Conference on Machine Learning

(PMLR, 2022) pp. 21864–21883.

[36] G. Braun, S. Pokutta, and D. Zink, Lazifying conditional gradient algorithms, in Proceedings of the 34th International Conference on Machine Learning (2017) pp. 566–575.

[37] J. Halbey, S. Rakotomandimby, M. Besançon, S. Designolle, and S. Pokutta, Efficient quadratic corrections for Frank-Wolfe algorithms, arXiv:2506.02635 (2025).

[38] P. Wolfe, Finding the nearest point in a polytope, Math. Program. 11, 128 (1976).

[39] O. Gühne and M. Seevinck, Separability criteria for genuine multiparticle entanglement, New J. Phys. 12, 053002 (2010).

[40] X.-Y. Chen, L.-Z. Jiang, and Z.-A. Xu, Necessary and sufficient criterion for *k*-separability of *n*-qubit noisy GHZ states, Int. J. Quantum Inf. 16, 1850037 (2018).

[41] T. Gao and Y. Hong, Detection of genuinely entangled and nonseparable *n*-partite quantum states, Phys. Rev. A 82, 062113 (2010).

[42] B. Jungnitsch, T. Moroder, and O. Gühne, Taming multiparticle entanglement, Phys. Rev. Lett. 106, 190502 (2011).

[43] L.-L. Ge, M. Xu, T. Hu, L.-Z. Jiang, and X.-Y. Chen, Tripartite separability of four-qubit W and Dicke mixed state in noise environment, Eur. Phys. J. Plus 136, 1 (2021).

[44] N. Ananth and M. Senthilvelan, On the non-k-separability of Dicke class of states and n-qudit w states, Int. J. Theor. Phys. 55, 1854 (2016).

[45] G. Aubrun and S. J. Szarek, Alice and Bob meet Banach: the interface of asymptotic geometric analysis and quantum information theory, Mathematical Surveys and Monographs, Volume 223 (American Mathematical Society, 2017).

[46] G. Pisier, The Volume of Convex Bodies and Banach Space Geometry, Cambridge Tracts in Mathematics (Cambridge University Press, 1989).

[47] H. Edelsbrunner and D. R. Grayson, Edgewise subdivision of a simplex, Discrete Comput. Geom. 24, 707 (2000).

[48] L. Gurvits and H. Barnum, Separable balls around the maximally mixed multipartite quantum states, Phys. Rev. A 68, 042312 (2003).
Appendix A: Frank-Wolfe method

Conditional Gradient (CG) methods, also known as Frank-Wolfe (FW) methods [24, 25], are a class of projection-free first-order algorithms for constrained convex optimization problems of the form

$$ \min_{x \in \mathcal{X}} f(x), $$

where f is a smooth convex function and $\mathcal{X}$ is a compact convex set. Instead of performing projections onto $\mathcal{X}$, which can be computationally expensive, CG methods access the feasible region via a Linear Minimization Oracle (LMO). At each iteration, the algorithm solves a linear subproblem to compute the so-called FW vertex

$$ v_t = \operatorname*{argmin}_{v \in \mathcal{X}} \langle \nabla f(x_t), v \rangle. $$

The next iterate is then a convex combination of $v_t$ and $x_t$, which ensures feasibility. While FW methods are attractive for their simplicity and projection-free structure, their vanilla form requires an LMO call at every iteration and achieves only a sublinear convergence rate of $O(1/t)$ in the general convex case.

To improve convergence, active-set variants of CG methods maintain an explicit convex combination of the current iterate over a set of previously used vertices, the so-called active set $A_t$. Notable examples include Away-step Frank-Wolfe [33], Fully-Corrective Frank-Wolfe [34], and Blended Pairwise Conditional Gradients (BPCG) [35]. Additionally to the standard FW step, these methods perform local updates changing the weights of the vertices. This enables dropping previously selected vertices from the active set, leading to sparser solutions and faster convergence. In the case of polytopes and strongly convex objectives, these methods can achieve linear convergence.

BPCG is a particularly efficient instance of the Frank-Wolfe method. In each iteration, it computes the local FW vertex $s_t$ and the away vertex $a_t$,

$$ s_t = \operatorname*{argmin}_{v \in A_t} \langle \nabla f(x_t), v \rangle, $$

$$ a_t = \operatorname*{argmax}_{v \in A_t} \langle \nabla f(x_t), v \rangle. $$

Unlike the global FW vertex, computing these vertices does not require a call to the LMO but can be done efficiently by enumerating the active set.

BPCG then compares the standard FW direction $d_t^{FW} = v_t - x_t$ with the local pairwise direction $d_t^{LPW} = s_t - a_t$ based on their respective gap value. The Frank-Wolfe gap $g_t^{FW} = \langle \nabla f(x_t), x_t - v_t \rangle$ and the local pairwise gap $g_t^{LPW} = \langle \nabla f(x_t), a_t - s_t \rangle$ measure the alignment of each direction with the negative gradient and thus their potential to decrease the objective. BPCG proceeds in the direction corresponding to the larger gap.

Furthermore, BPCG uses a lazification technique inspired by Braun et al. [36] storing an estimate of the Frank-Wolfe gap $\Phi_t$. By this, one avoids LMO calls in iterations where the local update is used. This is especially useful in the context of complicated feasible regions, where the LMO is expensive, which is precisely our case in this work.

To accelerate our method even further, we equip BPCG with quadratic correction steps (QC) recently proposed in Halbey et al. [37]. In cooperation with our group they conducted preliminary experiments on our problem and included the results in their work. The QC step is tailored for quadratic objectives, $f(x) = \frac{1}{2}x^T A x + b^T x$, like the considered projection problem. While the local pairwise step updates only two weights in each iteration, QC aims to optimize the weights of the whole active set in a single step. Minimizing a given objective over the convex hull of the current active set is in general as hard as the original problem. However, minimizing a quadratic objective over the affine hull of given vertices can be done by solving a linear system,

$$ \sum_{u \in A_t} \lambda[u] \langle Au, v - w \rangle = 0 \quad \forall v \in A_t \setminus \{w\}, \quad \sum_{u \in A_t} \lambda[u] = 1, $$

where $\lambda[u]$ denotes the barycentric coordinate of the vertex $u$ and $w$ is a fixed vertex in $A_t$. The affine minimizer only lies in the convex hull of the active set (and is therefore feasible) if and only if the barycentric coordinates $\lambda[u]$ are all non-negative.

To handle this issue we follow the QC-MNP approach inspired by the Minimum Norm Point (MNP) algorithm by Wolfe [38]. One considers the line segment between the current iterate and the affine minimizer. Using a simple ration test, one can find the point on the line segment and in the convex hull that is closest to the affine minimizer. In practice, we perform a QC-step every time a new vertex is added to the active set, otherwise we continue with cheap pairwise updates in BPCG.

The framework of the algorithm we actually used for entanglement characterization is as follows:
Lazified BPCG with QC-MNP

**Input:** Distance function $f(x)$; initial point $x_0 \in S_k$

1: $\Phi_0 \leftarrow \max_{v \in \mathcal{X}} \langle \nabla f(x_0), x_0 - v \rangle / 2$

2: $A_0 \leftarrow \{x_0\}$

3: $b = \text{false}$

4: **for** $t = 0, \dots, T$ **do**

5: $a_t \leftarrow \operatorname*{argmax}_{v \in A_t} \langle \nabla f(x_t), v \rangle$

6: $s_t \leftarrow \operatorname*{argmin}_{v \in A_t} \langle \nabla f(x_t), v \rangle$

7: **if** $\langle \nabla f(x_t), a_t - s_t \rangle \ge \Phi_t$ **then**

8: **if** $b$ **then**

9: $x_{t+1}, A_{t+1} \leftarrow \text{QC-MNP}(A_t)$

10: $\Phi_{t+1} \leftarrow \Phi_t$

11: $b \leftarrow \text{false}$

12: **else**

▷ away vertex

▷ local FW vertex

13: $d_t \leftarrow a_t - s_t$

14: $\gamma_t^* \leftarrow \lambda[a_t]$

15: $\gamma_t \leftarrow \operatorname*{argmin}_{\gamma \in [0, \gamma_t^*]} f(x_t - \gamma d_t)$

16: $x_{t+1} \leftarrow x_t - \gamma_t d_t$

17: $\Phi_{t+1} \leftarrow \Phi_t$

18: **if** $\gamma_t < \gamma_t^*$ **then**

19: $A_{t+1} \leftarrow A_t$

20: **else**

21: $A_{t+1} \leftarrow A_t \setminus \{a_t\}$

22: **end if**

23: **end if**

24: **else**

25: $v_t \leftarrow \operatorname*{argmin}_{v \in \mathcal{X}} \langle \nabla f(x_t), v \rangle$

26: $b \leftarrow (v_t \in A_t)$

27: **if** $\langle \nabla f(x_t), x_t - v_t \rangle \ge \Phi_t / 2$ **then**

28: $d_t \leftarrow x_t - v_t$

29: $\gamma_t \leftarrow \operatorname*{argmin}_{\gamma \in [0,1]} f(x_t - \gamma d_t)$

30: $x_{t+1} \leftarrow x_t - \gamma_t d_t$

31: $\Phi_{t+1} \leftarrow \Phi_t$

32: $A_{t+1} \leftarrow A_t \cup \{v_t\}$

33: **else**

34: $x_{t+1} \leftarrow x_t$

35: $\Phi_{t+1} \leftarrow \Phi_t / 2$

36: $A_{t+1} \leftarrow A_t$

37: **end if**

38: **end if**

39: **end for**

**Output:** $x_1, \dots, x_T \in \mathcal{X}$

▷ local pairwise update

▷ descent step

▷ drop step

▷ global FW vertex

28: $d_t \leftarrow x_t - v_t$

29: $\gamma_t \leftarrow \operatorname*{argmin}_{\gamma \in [0,1]} f(x_t - \gamma d_t)$

30: $x_{t+1} \leftarrow x_t - \gamma_t d_t$

31: $\Phi_{t+1} \leftarrow \Phi_t$

32: $A_{t+1} \leftarrow A_t \cup \{v_t\}$

33: **else**

34: $x_{t+1} \leftarrow x_t$

35: $\Phi_{t+1} \leftarrow \Phi_t / 2$

36: $A_{t+1} \leftarrow A_t$

37: **end if**

38: **end if**

39: **end for**

▷ FW step

▷ gap step

Appendix B: Examples for detecting different entanglement structures

Here, we recall the definition of multipartite entanglement. A multipartite pure state is called *k*-separable if it can be written as $|\psi^{ks}\rangle = |\phi_1\rangle \otimes |\phi_2\rangle \otimes \dots \otimes |\phi_t\rangle$, where each $|\phi_i\rangle$ may itself be entangled within a subset of the total system. A mixed state is *k*-separable if it can be expressed as a convex combination of *k*-separable pure states: $\sigma^{ks} = \sum_i \lambda_i |\psi_i^{ks}\rangle \langle \psi_i^{ks}|$, where the coefficients $\lambda_i$ form a probability distribution. Usually, if *n*-partite quantum states are not *n*-separable, then they are called entangled states. And non-2-separable states are called GME.

In the main text, we focus on the detection of multipartite entanglement of *n*-qubit states against full separability. In this section, we illustrate that CG methods can naturally handle rich structures of multipartite entanglement. We consider the white noise robustness problem for $S_k$ with $k = 2$ and 3 for several common quantum states. Note, when $k = 2$, it's the case of genuine multipartite entanglement (GME), where detection of GME is an important problem in the field.
<table><caption>Table I. Entanglement robustness thresholds for several multipartite quantum states against white noise, with respect to k-separability for $k = 2$ (GME) and $k = 3$. Each entry shows the entanglement bound $p_{ent}$ and the gap $\Delta p$ to separable bound, obtained by our CG-based framework. Known entanglement and separable thresholds from the literature are given in parentheses for comparison. One number in parentheses represents the exact thresholds.</caption><thead><tr><th rowspan="2">Quantum states</th><th colspan="2">Entanglement bound ± gap</th></tr><tr><th>k = 2 (GME)</th><th>k = 3</th></tr></thead><tbody><tr><td>GHZ(5 Qubits)</td><td>0.508 + 0.016 (<sup>16</sup>&frasl;<sub>31</sub> ≈ 0.516 [39])</td><td>0.760 + 0.006 (<sup>16</sup>&frasl;<sub>21</sub> ≈ 0.762 [40])</td></tr><tr><td>GHZ(4 Qubits)</td><td>0.528 + 0.009 (<sup>8</sup>&frasl;<sub>15</sub> ≈ 0.533 [39])</td><td>0.797 + 0.005 (<sup>4</sup>&frasl;<sub>5</sub> = 0.8 [40])</td></tr><tr><td>W(5 Qubits)</td><td>0.515 + 0.070 (0.478 [41], -)</td><td>0.739 + 0.013 (-, -)</td></tr><tr><td>W(4 Qubits)</td><td>0.528 + 0.003 (0.526 [42], 0.545 [18])</td><td>0.722 + 0.033 (<sup>216-16√6</sup>&frasl;<sub>235</sub> ≈ 0.752 [43])</td></tr><tr><td>Dicke(5 Qubits, 2 ex.)</td><td>0.518 + 0.047 (0.516 [44], -)</td><td>0.717 + 0.055 (0.615, -)</td></tr><tr><td>Dicke(4 Qubits, 2 ex.)</td><td>0.540 + 0.003 (0.539 [42], -)</td><td>0.769 + 0.032 (0.788, -)</td></tr></tbody></table>

## Appendix C: Construction on an $\varepsilon$-net of the set of pure separable states

With [45] and in particular Lemma 9.2 therein, the construction of an $\varepsilon$-net of the set of pure separable states can be reduced to the construction of an $\varepsilon$-net of a corresponding hypersphere. In this appendix we then concentrate on this case, giving an explicit construction, in contrast with the arguments used in [45].

### 1. Definitions and notations

Let $H$ be a $d$-dimensional normed space with unit ball $B$ and unit sphere $S$. Let $K$ be a nonempty centrosymmetric compact convex set in $H$ and $X = \{x_i\}_{i=1}^m$ be a subset of $\text{ext}(K)$ such that $0 \in \text{conv}(X) = P$.

**Definition 5.** The homothetic distance, also called the geometric distance or the shrinking factor, of $X$ with respect to $K$ is the maximum $\eta$ such that $\eta K \subset P$.

**Definition 6.** $X$ is an inner $\varepsilon$-approximation of $K$ if there exists $\varepsilon > 0$ such that $\frac{1}{1+\varepsilon}K \subset P$.

**Definition 7.** $X$ is an $\varepsilon$-net of $\text{ext}(K)$ if there exists $\varepsilon > 0$ such that $\text{ext}(K) \subset X + \varepsilon B$.

In this appendix we focus on the case where $H$ is real and $K = B$ is the unit ball.

**Observation 8.** An inner $\varepsilon$-approximation of $B$ has homothetic distance $\eta \ge \frac{1}{1+\varepsilon}$ with respect to $B$.

**Observation 9.** An $\varepsilon$-net of $S$ has homothetic distance $\eta \ge \sqrt{1 - \frac{\varepsilon^2}{2}}$ with respect to $B$.

### 2. Implicit construction

In [46, Lemma 4.10], the following result is stated, whose proof we recall for completeness.

**Proposition 10.** There exists an $\varepsilon$-net of $S$ with cardinality smaller than $(1 + \frac{2}{\varepsilon})^d$.

*Proof.* Let $X$ be a maximal subset of $S$ such that $\|x_i - x_j\| \ge \varepsilon$ for all $i \ne j$. Clearly, by maximality, $X$ is an $\varepsilon$-net of $S$. To derive an upper bound on its cardinality $m$, we note that the balls $x_i + \frac{\varepsilon}{2}B$ are disjoint and included into $(1 + \frac{\varepsilon}{2})B$. Therefore,

$$ \sum_{i=1}^{m} \mathrm{vol} \left( x_i + \frac{\varepsilon}{2} B \right) \le \mathrm{vol} \left[ \left( 1 + \frac{\varepsilon}{2} \right) B \right], \quad \text{hence} \quad m \left( \frac{\varepsilon}{2} \right)^d \mathrm{vol}(B) \le \left( 1 + \frac{\varepsilon}{2} \right)^d \mathrm{vol}(B), \quad (C1) $$

and this implies $m \le (1 + \frac{2}{\varepsilon})^d$ as announced. $\square$

To the best of our knowledge, no better scaling is currently known, but the construction is implicit, so that its usage in practice is somehow limited. In the next subsection, we present an alternative construction that is fully explicit, and whose scaling is even better.
### 3. Explicit construction

We first recall without proof the main result from [47] on the edgewise subdivision of a simplex.

**Proposition 11.** For every integer $n \ge 1$ and every $d$-simplex $\sigma$ there is a subdivision of $\sigma$ into $n^d$ $d$-simplices of equal volume.

**Theorem 12.** Let $X \subset S$ be an inner $\epsilon$-approximation of $B$ and consider a simplicial complex triangulating the facets of $P = \text{conv}(X)$ and such that the circumcenter of each simplex lies inside this simplex. For an integer $n \ge 2$, construct the set $X_n$ by applying Proposition 11 to all of these simplices and projecting the resulting vertices onto $S$. Then $X_n$ is an inner $\epsilon_n$-approximation of $B$ for $\epsilon_n$ such that

$$ (1 + \epsilon_n)^2 = 1 + \frac{d-1}{2d} D_n^2 (1 + \epsilon)^2, \qquad (C2) $$

where $D_n$ is the maximal diameter of the subsimplices before projection onto $S$.

*Proof.* Let $\Delta$ be one of the $(d-1)$-simplices obtained after subdivision of the facets, but before projection onto $S$. Denote $c$ the center of its circumscribed sphere (which by assumption is also the center of its minimum enclosing ball) and $a$ one point on this sphere. By Jung's theorem, we have that

$$ \|a - c\| \le \frac{d-1}{2d} \text{diam}(\Delta)^2 \le \frac{d-1}{2d} D_n^2. \qquad (C3) $$

Since $c$ lies on a facet of $P$, we also have $\|c\| \ge \frac{1}{1+\epsilon}$. Combining these two inequalities and by the Pythagorean theorem (see Fig. 5), we obtain

$$ \left\| \frac{c}{\|a\|} \right\|^2 = \frac{\|c\|^2}{\|c\|^2 + \|a-c\|^2} \ge \frac{\left(\frac{1}{1+\epsilon}\right)^2}{\left(\frac{1}{1+\epsilon}\right)^2 + \frac{d-1}{2d} D_n^2}. \qquad (C4) $$

Since the closest point to the origin of the projection of $\Delta$ onto $S$ is $\frac{c}{\|a\|}$, this bound is valid for $X_n$. $\square$

<figure><img src="image_6.png" /><figcaption>Figure 5. Setup for the proof of Theorem 12.</figcaption></figure>

**Proposition 13.** Consider a regular $(d-1)$-simplex with edge length $s$. For $d \le 3$ or $n = 1$, the diameters of its subsimplices are all $\frac{s}{n}$; in the other cases, they are at most $\frac{s\sqrt{2}}{n}$.

*Proof.* The diameter of a simplex is simply the largest distance between two of its vertices, which immediately gives $D_1 = s$. From the proof of the independence lemma in [47], we know that the vectors defining the subsimplices are, up to a factor $n$, the same as those defining the initial simplex, the various shapes coming from a permutation of the order of these vectors. For $d \le 3$, this permutation does not affect the edge length of the different shapes. For higher dimensions, the worst case gives a factor $\sqrt{2}$ by the Pythagorean theorem. $\square$
**Example 14.** Let $X$ be the set of vertices of the $d$-dimensional cross-polytope. Observe first that its homothetic distance with respect to the unit ball is $\eta = \frac{1}{\sqrt{d}}$, corresponding to $\epsilon = \sqrt{d} - 1$. Then choose an integer $n \ge 2$ and define $X_n$ according to the procedure in Theorem 12. All $2^d$ facets of $P$ are regular $(d-1)$-simplices with edge length $\sqrt{2}$. Each facet gives rise to $\binom{n+d}{d}$ vertices of $X_n$ (the $(n+1)$-th $(d-1)$-simplex number), so that the total cardinality of $X_n$ is lower than $2^d \binom{n+d}{d}$. With Proposition 13 and Theorem 12, we obtain that the homothetic distance $\eta_n$ of $X_n$ with respect to the unit ball is such that

$$ \eta_n \ge \frac{n}{\sqrt{n^2 + 2(d-1)}}, \qquad (C5) $$

so that $X_n$ is an $\epsilon$-net of $S$ for

$$ \epsilon = \sqrt{\frac{2}{1 + \frac{n^2}{2(d-1)}}} \sim \frac{2\sqrt{d-1}}{n}. \qquad (C6) $$

This gives the following scaling for the cardinality of this family of $\epsilon$-nets:

$$ \frac{2}{(d-1)!} \left( \frac{4\sqrt{d-1}}{\epsilon} \right)^{d-1}, \qquad (C7) $$

which improves on Proposition 10 as announced.

## Appendix D: Proof of Proposition 3

Proposition 3 is built on the construction on $\epsilon$-net of quantum space, which is discussed above. As the product state $\phi$ is the closest state to the direction $\Lambda$ in $S_k^\epsilon$, there is

$$ \mathrm{tr}(\Lambda|\phi_i\rangle\langle\phi_i|) \le \mathrm{tr}(\Lambda|\phi\rangle\langle\phi|) = \beta, \quad \forall \phi_i \in \mathrm{conv}(S_k^\epsilon). \qquad (D1) $$

According to Definition 5, we have

$$ \eta S_k \subset \mathrm{conv}(S_k^\epsilon) \subset S_k. \qquad (D2) $$

Thus, for every $\tau \in S_k$, we have

$$ \begin{aligned} \mathrm{tr}(\Lambda\tau) &= \mathrm{tr}(\Lambda\eta\tau) + (1-\eta)\mathrm{tr}(\Lambda\tau) \\ &\le \mathrm{tr}(\Lambda|\phi^*\rangle\langle\phi^*|) + (1-\eta)\mathrm{tr}(\Lambda\tau) \end{aligned} \qquad (D3) \\ &\le \mathrm{tr}(\Lambda|\phi^*\rangle\langle\phi^*|) - (1-\eta)\|\Lambda\|. $$

Note the final inequality applying the Cauchy-Schwarz inequality and $\mathrm{tr}(\tau) = 1$. Thus we obtain the witness $\mathcal{W}$ in Proposition 3 such that

$$ \mathrm{tr}(\mathcal{W}\tau) \ge 0, \quad \forall \tau \in S_k. \qquad (D4) $$

## Appendix E: Proof of Proposition 4

*Proof.* In Ref. [18], a geometric reconstruction approach was developed within the CG-based framework for separability certification. In short, based on a quantum state $\rho$, we can construct a state $\rho_t = (1 - \epsilon)\rho - \epsilon \mathbb{1}/d$. For the constructed state $\rho_t$, the CG method can give us the separable state $\sigma_t$ close to $\rho_t$, as a numerical result. Extrapolating $\rho$ against $\rho_t$ to achieve $\rho_x$ such that the direction of $\rho_t - \sigma_t$ and $\rho_x - \mathbb{1}/d$ are parallel. Then if the distance $r_t = \|\rho_t - \sigma_t\|$ is smaller than the radius $a$ of the separate ball [48], the state $\rho_x$ is separable, therefore the state $\rho$ is separable as well. An illustration of the geometric reconstruction is shown in Fig. 6.

For the problem of estimating the white noise robustness of an entangled state $|\psi\rangle$, let the constructed state be the noisy state we consider, i.e.,

$$ \rho_t = \rho(p) = (1-p)|\psi\rangle\langle\psi| + p\mathbb{1}/d^2. \qquad (E1) $$
<figure><img src="image_7.png" /><figcaption>Figure 6. Geometric reconstruction for separability certification within the CG-based framework.</figcaption></figure>

Then as long as

$$r_t = \frac{1}{\epsilon} \| \rho_t - \sigma_t \| \le a, \qquad (E2)$$

the state

$$\rho = \frac{1}{1+\epsilon}\rho_t + \frac{\epsilon}{1+\epsilon}\frac{\mathbf{1}}{d} = \left(1 - \frac{\epsilon+p}{1+\epsilon}\right)|\psi\rangle\langle\psi| + \frac{\epsilon+p}{1+\epsilon}\frac{\mathbf{1}}{d} = \rho\left(\frac{p+\epsilon}{1+\epsilon}\right) \qquad (E3)$$

will be separable.

<figure><img src="image_8.png" /></figure>

# Appendix F: Details of applications

## 1. Horodecki state

The explicit density matrix of the family of Horodecki states [4] is

$$\rho^H(a) = \frac{1}{8a+1} \begin{pmatrix} a & \cdots & a & \cdots & a \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ a & \cdots & a & \cdots & a \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ a & \cdots & a & \cdots & a \end{pmatrix} + \begin{pmatrix} \frac{1+a}{2} & \cdots & \frac{\sqrt{1-a^2}}{2} & \cdots & \frac{1+a}{2} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \frac{1+a}{2} & \cdots & \frac{\sqrt{1-a^2}}{2} & \cdots & \frac{1+a}{2} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ \frac{1+a}{2} & \cdots & \frac{\sqrt{1-a^2}}{2} & \cdots & \frac{1+a}{2} \end{pmatrix}, \qquad (F1)$$

where the dots indicate zeroes.

## 2. Noise channels

Here we explicitly give the definitions of the noise channels we consider in the main text, including (a) global depolarizing (white noise), (b) bit flip, (c) phase flip, (d) amplitude damping, and (e) phase damping channels.

For one qubit, the depolarizing channel is also known as white noise, i.e.,

$$C_D(\rho) = (1-p)\rho + p\frac{\mathbf{1}}{2}. \qquad (F2)$$
For a bipartite state, the white noise is the global depolarizing channel. And under a balance local depolarizing channel, the noisy state will be

$$C_{\text{GD}}(\rho) = (1-p)^2\rho + p(1-p)\text{tr}_A(\rho) \otimes \frac{\mathbf{1}}{2} \\ + p(1-p)\frac{\mathbf{1}}{2} \otimes \text{tr}_B(\rho) + p^2\frac{\mathbf{1}}{2} \otimes \frac{\mathbf{1}}{2}. \qquad (\text{F3})$$

Hence, for a maximally entangled bipartite state, where $\text{tr}_A(\rho) = \text{tr}_B(\rho) = \mathbf{1}/d$, the global and local depolarizing channels are the same. We have known the robustness threshold $p = 1/(d+1)$ for the local dimension $d$.

For the bit flip channel and the other channel below, we describe them by Kraus operators, i.e.,

$$C(\rho) = \sum_{i} E_{i} \rho E_{i}^{\dagger}, \qquad (\text{F4})$$

where Kraus operators need to satisfy $\sum_i E_i E_i^\dagger \le \mathbf{1}$.

For the bit flip channel, the Kraus operators are

$$C_{\text{BF}}(\rho) = E_0\rho E_0^\dagger + E_1\rho E_1^\dagger, \\ E_0 = \sqrt{1-p}\mathbf{1}, \quad E_1 = \sqrt{p}X. \qquad (\text{F5})$$

The phase flip channel can be described as

$$C_{\text{PF}}(\rho) = E_0\rho E_0^\dagger + E_1\rho E_1^\dagger, \\ E_0 = \sqrt{1-p}\mathbf{1}, \quad E_1 = \sqrt{p}Z. \qquad (\text{F6})$$

For the amplitude damping channel, it is the description of energy dissipation, that is, effects due to loss of energy from a quantum system. For one qubit, it is defined as

$$C_{\text{AD}}(\rho) = E_0\rho E_0^\dagger + E_1\rho E_1^\dagger, \\ E_0 = |0\rangle\langle 0| + \sqrt{1-\gamma}|1\rangle\langle 1|, \quad E_1 = \sqrt{\gamma}|0\rangle\langle 1|, \qquad (\text{F7})$$

where $\gamma$ can be thought of as the probability of losing a photon.

The phase damping channel describes the loss of quantum information without loss of energy, which happens when a photon scatters randomly as it travels through a waveguide, or how electronic states in an atom are perturbed upon interacting with distant electrical charges. The channel for qubit systems is

$$C_{\text{PD}}(\rho) = E_0\rho E_0^\dagger + E_1\rho E_1^\dagger, \\ E_0 = |0\rangle\langle 0| + \sqrt{1-\gamma}|1\rangle\langle 1|, \quad E_1 = \sqrt{\gamma}|1\rangle\langle 1|, \qquad (\text{F8})$$

The numerical results for Fig. 4 are presented below.

<table><caption>Table II. Entanglement robustness thresholds for the Bell state ($|00\rangle + |11\rangle)/\sqrt{2}$ under various quantum noise channels: global depolarizing (white noise, GD), bit flip (BF), phase flip (PF), amplitude damping (AD), and phase damping (PD). Except for white noise (applied to both qubits), all channels act on one side of the state. Error bars represent the certified interval between the entanglement and separability bounds obtained by our numerical algorithm.</caption><thead><tr><th>Noise Channel</th><th>Noise Strength</th><th>Fidelity Threshold</th></tr></thead><tbody><tr><td>White Noise</td><td>0.6665 + 0.0005 (3/3)</td><td>50.01% - 0.04% (50%)</td></tr><tr><td>Bit Flip</td><td>0.4999 + 0.0002 (1/2)</td><td>50.01% - 0.02% (50%)</td></tr><tr><td>Phase Flip</td><td>0.4999 + 0.0002 (1/2)</td><td>50.01% - 0.02% (50%)</td></tr><tr><td>Amplitude Damping</td><td>0.9997 + 0.0003 (1)</td><td>25.87% - 0.87% (25%)</td></tr><tr><td>Phase Damping</td><td>0.9999 + 0.0001 (1)</td><td>50.50% - 0.50% (50%)</td></tr></tbody></table>