Certifying quantum separability with adaptive polytopes

The concept of entanglement and separability of quantum states is relevant for several fields in physics. Still, there is a lack of effective operational methods to characterise these features. We propose a method to certify quantum separability of two-and multi-particle quantum systems based on an adaptive polytope approximation. This leads to an algorithm which, for practical purposes, conclusively recognises two-particle separability for small and medium-size dimensions. For multiparticle systems, the approach allows to characterise full separability for up to five qubits or three qutrits; in addition, different classes of entanglement can be distinguished. Finally, our methods allow to identify systematically quantum states with interesting entanglement properties, such as maximally robust states which are separable for all bipartitions, but not fully separable. Copyright T.-A. Ohst et al . This work is licensed under the Creative Commons Attribution 4.0 International License.


Introduction
Entanglement is nowadays perceived as one of the hallmarks of quantum mechanics, which not only underlies the foundations of quantum mechanics and quantum technology, but also deeply influences our understanding of physics in various different areas, ranging from condensed matter physics [1,2] to gravity [3][4][5].Formally, entangled states are those that cannot be prepared by classical communication and local operations of the parties [6].Modelling the classical communication by a random variable λ, this implies that a bipartite entangled state ρ AB cannot be written as a convex combination of states factorised at the parties, where σ λ and τ λ are density operators in Alice's and Bob's spaces, respectively, and p λ are probability weights of λ.States of the form (1) are said to be separable.Significant effort has been devoted to methods to determine whether a state is entangled or not [7][8][9][10][11][12][13].As separable states form a convex set, depicted in Fig. 1a, a state can be 'witnessed' to be entangled if one can find a hyperplane that separates it from the set.On the contrary, proving a state to be separable is significantly complicated, requiring testing it against all possible entanglement witnesses or, equivalently, searching over all possible decompositions of the form (1). Accordingly, various methods have been proposed to demonstrate entanglement [7,8], but techniques to verify that a state is separable are scarce.In fact, it has been perceived as a 'notoriously difficult' problem in quantum information theory [14,15].On the other hand, certification of separability can be essential to optimise quantum information processing protocols, where entanglement between certain parties is difficult to establish and only local operations are available.This has been discussed, for example, in conference key agreement [16], superdense coding [17] and reconstruction of states from marginals [18,19].
Known approaches to the problem are based on sequences of semidefinite programs that are proven to certify separability for sufficiently high order [20,21].Later, iterative algorithms in which states are tried to be transferred into maximally mixed ones with infinitesimal non-separability-breaking transformations have been applied [22,23].Further works used a version of Gilbert's algorithm to the convex membership problem [24,25], the method of truncated moment sequences [26] and sets of inequalities in terms of Bloch representations [27].Recently, neural networks have been used for a parametrisation of separable states to tackle the problem of certifying separability [28] and even variational quantum algorithms for the problem exist [29].These present methods are highly computationally demanding, therefore applicable only for special states or low-dimensional systems.Filling this gap, we introduce a method of adaptive polytopes for certifying the separability of quantum states.We show that the resulting algorithm indicates a strong evidence of being nearly optimal independent of the structure of the states in all benchmarks.Being highly efficient, the algorithm is directly applicable to quantum systems of relatively high dimensions and systems with many particles, far beyond results given by the other known methods.In fact, the algorithm not only allows for the certification of a single targeted state, but also for families associated to it, and even for the investigation of different entanglement robustnesses.As an illustration, we apply the technique to explore the geometry of the boundary of the set of separable states for bipartite as well as multiparticle systems.

Polytope approximation
We start with rewriting the set of separable states in Eq. ( 1), denoted SEP, as the convex hull of product states, where S A ⊗ S B = {σ ⊗ τ : σ ∈ S A , τ ∈ S B }, with S A and S B being the sets of states at Alice's side and at Bob's side, which we will call the generalised Bloch spheres.Let P be a convex subset of Alice's operators with unit trace.Following (2), we consider conv(P ⊗ S B ).If we choose an inner polytope P in and an outer one P out to approximate Alice's generalised Bloch sphere, such that P in ⊆ S A ⊆ P out (see Fig. 1b for an illustration) then it follows that While the set SEP = conv(S A ⊗S B ) cannot be efficiently computed, we show that conv(P in/out ⊗S B ) can be.In fact, they can be formulated as a standard optimisation of a linear objective function and semidefinite constraints, known as a semidefinite program (SDP), for which efficient algorithms exist [30].Indeed, let the polytope P be described by a set of N vertices, P = {σ λ } N λ=1 , then ρ AB ∈ conv(P ⊗ S B ) means that there exist N positive operators on Bob's space {τ λ } N λ=1 such that ρ AB = N λ=1 σ λ ⊗ τλ where the τλ do not necessarily have trace 1.This can be understood as a minimisation of a constant function with semidefinite constraints, known as a feasibility SDP.More quantitatively, approximating the Bloch sphere by a polytope from inside and outside allows one to directly lower bound and upper bound various types of entanglement robustnesses, see Appendix B. As an illustration, we consider here the white noise mixing threshold to an entangled state such that it becomes separable, where ρ AB t = tρ AB + (1 − t)1/d AB and d AB is the total dimension of the system.By choosing a polytope P = {σ λ } N λ=1 to approximate S A , one obtains an SDP to approximate χ(ρ AB ) by In practice, if P approximates the generalised Bloch sphere from inside (outside), one obtains a lower (upper) bound of χ(ρ AB ), see also Fig. 1a.Outer polytopes give rise to a useful tool to detect entanglement if one of the parties is a qubit and the dual version of the SDP (5) can then be used to construct tailored entanglement witnesses [31,32], see Appendix A. This is because the simplicity of the geometry the Bloch sphere of the qubit allows for an efficient construction of the outer polytopes approximating it.This is no longer the case for systems of qudits.For that reason, in the following, we concentrate on lower bounds of χ(ρ AB ) by the inner polytope approximation for separability certification.Such inner polytope approximations can be used for an accurate and efficient certification of separability, even if the local dimensions are larger than two.

Adaptive polytopes
A key insight is that the structure of the procedure allows us to iteratively improve the choice of the inner polytope.To be precise, starting with a polytope P at Alice's side, one finds the set of unnormalised states {τ λ } N λ=1 at Bob's side in Eq. (5).Upon normalisation, these give an inner polytope Q to approximate the Bloch sphere at Bob's side, which can be used for the algorithm (5) with Alice and Bob being interchanged and after performing a system swap on the state ρ.Importantly, this newly obtained polytope forms an approximation that is at least as strong as the preceding polytope, which can be seen from the symmetrical roles of τλ and σ λ in the optimisation problem (5).Conceptually, it also makes no difference whether Alice's and Bob's system have the same or different dimensions.In the latter case, the size of the matrices τλ simply switches between every round of the algorithm.The polytope adaption algorithm may be summarised in the following way: 1. Initialise an arbitrary inner polytope at Alice's side P A = {σ λ } N λ=1 .
3. Exchange A and B, use P B as polytope approximation and return to step 2.
The algorithm stops upon convergence of the visibility χ.These steps are illustrated in Fig. 2.This iterative procedure gives a series of SDPs approximating SEP with increasing accuracy, which in practice converges rather quickly.Also, no significant difference in the performance is observed when using symmetric and random polytopes as initiation; in the following, the latter are used.It is also important to mention that the algorithm has a mild increase in complexity when number of iterations, polytope vertices and the Hilbert space dimension get higher.The number of scalar variables for a single SDP Figure 2: Schematic visualisation of polytope adaptions.The optimal feasible points computed in the SDP in Eq. ( 5) where Alice's generalized Bloch sphere is approximated by a polytope are used to construct a new polytope for the next approximation of Bob's generalized Bloch sphere.In this way, the algorithm of polytope adaptions alternates between polytope approximations of Alice's and Bob's system. in one iteration step scales as O(N d 2 B ) for optimisation on Bob's side and as O(N d 2 A ) on Alice's side where N is the number of polytope vertices and d B is the local dimension of Bob.Each iteration step has the same number of variables so that the runtime increases linearly in the number of iterations.Our implementation is written in the Julia programming language [33] and the semidefinite programs were solved using the Mosek solver [34].The code is available on a public repository [35].
Our obtained inner approximation to SEP turns out to saturate various upper bounds such as those given by the positive partial transposition (PPT) criterion [36] or symmetric extensions [37] in all test cases, uncovering the optimality of both.For example, we obtain the exact values of χ for the isotropic and Werner states for local dimensions up to 10, and 10 4 random states of dimension 5 × 5 distributed according to the Hilbert-Schmidt measure [38] with a numerical accuracy of 10 −4 .The computation time for states of a 5 × 5 system takes around 2 − 3 seconds per iteration when the polytope has 200 vertices and 3 iterations are on average sufficient for obtaining the correct value.These times refer to a computer with a CPU Intel(R) Core(TM) i5-3570 processor (4 cores) running at 3.40GHz using 16 GB of RAM.

PPT-entangled states
We consider the two classic one-parameter families of quantum states of dimensions of 3 × 3 and 2×4 known as Horodecki states ρ H 3×3 (a) and ρ H 2×4 (b), where 0 ≤ a, b ≤ 1 [39]; for the explicit density operators see Appendix E. Both are entangled despite being PPT for 0 < a, b < 1 [39].
For the Horodecki state of dimension 3 × 3, many entanglement criteria have been used to obtain upper bounds for χ, and certain lower bounds are also known.As seen in Fig 3a, our lower bound outperforms the best known lower bound [25] and approaches the best known upper bound given in Ref. [40].For further comparison, we also implemented the SDP hierarchy of symmetric extensions to the fourth level for approximations of SEP from the outside [37].
In the case of ρ H 2×4 (b), the symmetric extension hierarchy is implemented up to level 5 on the qubit [37].In Fig. 3b, we observe that the lower bound given by our inner adaptive polytope algorithm nearly coincides with this upper bound up to the numerical accuracy, convincingly demonstrating its optimality.Moreover, as Alice's Bloch sphere is a 3-dimensional unit sphere, the outer polytope approximation can also be easily constructed, albeit without adaptation.We choose a fixed outer polytope approximation with 1002 vertices, which al-

(b) (a)
Figure 3: Plots of lower bound of the separability threshold for different state parameters of the 3 × 3 (a) and 2 × 4 (b) families of PPT-entangled states (red curves) obtained by polytope adaptions using 500 vertices.The dashed curves correspond to bounds given by different levels of symmetric extension (SE).The green and dotted curve in (a) labelled by "known lower" and "known upper" show the best previously known lower and upper bounds for the state family, see [25,40].
ready admits an entanglement detection generally better than symmetric extension of level 3, see Fig. 3b.
Importantly, the algorithm is applicable to generic states without special assumptions.This allows us to study and construct PPT entangled states beyond fine-tune families.Specifically, starting with a Horodecki PPT entangled state at a = 0.25 with visibility of χ[ρ H 2×4 (0.25)] = 0.9715, we design a see-saw algorithm to search for PPT states with higher entanglement robustness.In practice, we first compute the dual SDP of the polytope approximation (5) to obtain a witness for the corresponding approximation of SEP.In the second step, we maximise the violation of this witness over all PPT states.This new PPT entangled state is then chosen to be the input for the first SDP again.The see-saw algorithm stabilises at a state with significant lower visibility of 0.9461, which could be a candidate for an experimental realisation of robust bound entanglement.The method of polytope adaptions may therefore extend existing methods for the construction of bound entangled states [41].

Three-qubit systems
In multiparticle systems, one can distinguish between different types of entanglement and separability.Specifically, a state ρ ABC is called fully separable (FSEP) if it can be decomposed as convex combination of product states, Otherwise, it is entangled.The state ρ ABC is separable for the bipartition A|BC if it can be written as where τ BC λ may be entangled, and similarly for the remaining bipartitions AB|C and AC|B.States that are separable for any bipartition are called fully biseparable (FBSEP).Finally, the Figure 4: Two dimensional cross-sections of three-qubit states.The section represented in (a) corresponds to the plane that includes the GHZ-, W-and maximally mixed state.Note that the boundary for genuinely multiparticle entanglement is is already determined in Ref. [47].In (b) and (c), the planes under consideration are obtained by the maximally mixed state, some random state and the GHZ and W state and in (d), two random states are taken into account.set of biseparable (BSEP) quantum states are those which can be decomposed as a convex combination of states that are separable for these bipartions, where the superscripts indicate the membership to the corresponding separability class, e.g., ρ AB|C 1 ∈ SEP(AB|C).Quantum states which are not biseparable are called genuinely multiparticle entangled.
To apply the adaptive polytope method for FSEP, we introduce a polytope P A on the system A and demand that the vertices of P A are paired to positive-semidefinite operators on BC that are PPT.This approach is valid as the set of PPT-states coincides with the separable states for two qubits.
Checking the membership to BSEP amounts to a choice of suitable polytopes P A , P B and P C , one for each subsystem.Then, the feasible set in the SDP is given by the r.h.s. of Eq. ( 8) where each summand is replaced by the corresponding bipartite polytope approximation.
Similarly in the case of FBSEP, every subsystem is once approximated by a polytope but in opposition to BSEP it is demanded that the state ρ ABC is tested for membership to every bipartite separability class by a polytope approximation, see Appendix C.
As a demonstration, the inner approximation with a polytope of 300 vertices gives the white noise thresholds of 0.199 for full separability and 0.42857 for biseparability of the GHZ state in accordance with the known exact values [42,43].In the case of the W state, we obtain a full separability threshold of 0.178 and 0.479 for biseparability which are both known to be exact [40,44].Moreover, we are able to study certain two-dimensional cross-sections of the set of three qubit states, see Fig. 4. Remarkably, the introduced method delivers a suitable approximation for the set of fully separable states which is generally hard to achieve in contrast to the bipartite scenario.We use this advantage to study the set of fully biseparable states which are still entangled.Although examples for such states exist [45,46], the robustness of this phenomenon was unexplored.

Robust fully biseparable states
We apply our method to approach the question on which fully biseparable state is maximally robust against full separability.To answer this question, first random entangled states are initialised and their visibility to FSEP is computed with the dual SDP from the polytope approximation.The solution provides an approximate entanglement witness, whose inner product

State
Lower bound Upper bound GHZ ( ).An asterisk ( * ) indicates the bound known to be tight.A sharp (♯) indicates the value obtained by the PPT criterion.The values always correspond to full separability.The computation time for the extreme cases (5 qubits, 3 qutrits) is a few minutes with the same hardware specification as mentioned above.with all fully bi-separable states is minimised in the second step.Then, the FSEP-visibility of the resulting state is computed and a corresponding new approximate entanglement witness is determined.These steps are repeated until convergence, see Appendix F. Although this iteration is not guaranteed to terminate in a global optimum, we obtain a stable minimal separability threshold χ ≈ 0.57.Further analysis shows that the corresponding obtained state admits a surprisingly simple form, for angles θ ∈ [0, 2π) that are solutions to (sin θ cos θ ) 2 = 1/6, see Appendix F for a geometric analysis.The four vectors in expression (9) are with where ⊕ denotes addition modulo 2 and α = α ⊕ 1, see Appendix F for more details.

More parties and higher dimensions
The idea of the adaptive algorithm can be further extended to certify full separability for systems of more parties and higher dimensions.Consider an n-partite state ρ 1•••n .One can initiate a polytope approximation for the products states of n − 1 parties, One then can check whether there exist positive semidefinite matrices {τ n λ } N λ=1 for the last party such that This is again an SDP.The output of this SDP gives a polytope for the last party.One then iterates over the parties to systematically improve the approximation.In addition, for a pair of systems with total dimension lower than 6, one can use the PPT criterion to simplify the procedure.As an example, we study the full separability of up to five qubits as well as three qutrits.Comparison of the results with known states in the literature gives excellent agreement as indicated in Table 1.Extension to many-body quantum systems with more parties but with special structure or symmetry is interesting for future study.

Conclusion
We introduced a method to tackle the problem of certifying the separability of quantum states which is optimal under benchmarks.The resulting algorithm allows for a highly accurate description of the set of separable states from inside.The two main illustrative applications of the algorithm are the precise entanglement characterisation of states that cannot be found by the PPT criterion as well as a precise approximation of fully separable states in a multiparticle system of up to five qubits or three qutrits.Methodologically, we suggested a new approach to multi-linear optimisation problems that clearly demands for future applications.Concretely, related ideas can be useful in the study of quantum networks using local operations and shared randomness and in the characterisation of entanglement-breaking quantum channels.It is also directly applicable to the analysis of the role of memory in quantum processes building on the results in Ref. [50].Future applications that might reach far beyond entanglement theory are also within the realm of possibility.

A Semidefinite programming characterisation of conv(P ⊗ S B )
In this appendix, we demonstrate how one can characterise conv(P ⊗ S B ) for a polytope P by means of semidefinite programming.More specifically, let ρ AB be a bipartite state.We would like to determine if ρ AB ∈ conv(P ⊗ S B ). Quantitatively, we want to compute Let the polytope P be given by the vertices {σ λ } N λ=1 .The problem becomes For the sake of convenience, we formulate the dual problem for computing a related and equally useful quantity χ P (ρ AB ) = 1/[1 + r P (ρ AB )] with The correspondence of ( 14) and ( 15) can be seen in two steps.In the first step, we derive an expression of the set of non-normalised separating hyperplanes for conv(P ⊗ S B ), which is the dual of the membership problem to conv(P ⊗ S B ).The membership problem is and its dual reads min Tr(ρ AB W AB ) w.r.t.
This will then correspond to an outer approximation of the set of entanglement witnesses in the separability problem provided that P is an inner polytope.
In the second step, we use the representation of the entanglement robustness measure in terms of normalised entanglement witnesses derived in Ref. [31].There, it is shown that the robustness with respect to white noise r(ρ AB ) = min s : can be equivalently expressed as The set N on which the expression is minimised in (19) contains here entanglement witnesses with a special normalisation, namely In our case where we replace the set of states on A by a polytope, we have to replace N in (19) by N P , which is given by to obtain r P .We can conclude that χ P (ρ AB ) = 1/[1 + r P (ρ AB )] with which is in accordance to (15).

B Estimation of other robustness measures
In our analysis, we focus on determining the visibility, or equivalently the robustness of entanglement with respect to white noise defined in Eq. ( 18).This is also often called random robustness [51].In this appendix, we discuss other robustness measures of entanglement [51].
First, robustnesses resulting from different choices of separable states other than the maximally mixed state as noise mixing component in the random robustness (18) can be computed in a completely similar way.More generally, we can consider the so-called (absolute) robustness of entanglement with respect to the whole separable set, defined as [51] R(ρ AB ) = min Upon approximating Alice's set of states by a polytope P with vertices This can be computed by the following semidefinite program: In addition, the so-called generalised robustness, introduced in Ref. [52] as can efficiently be computed with polytope approximations.The difference to the absolute robustness is that also entangled states are taken into account as possible sources of noise.

C Polytope approximations in multiparticle systems
Turning to the multiparticle scenario, we start with considering the three-qubit system.We use α ∈ {A, B, C} to denote one of the three parties, and ᾱ to denote its other complementary parties.

C.1 Full separability
Mathematically, the set of fully separable states is defined by FSEP = conv(S A ⊗ S B ⊗ S C ).To construct the inner (outer) polytope approximation in this case, it is useful to observe that conv(S A ⊗ S B ⊗ S C ) = conv(S A ⊗ PPT BC ), where PPT BC denotes the set of two-qubit states that have a positive partial transpose, which is the same as SEP BC .By approximating S A by a polytope P = {σ A λ } N λ=1 from inside (outside), certification of a state ρ ABC to be fully separable amounts to finding λ ≥ 0 and This is an SDP in which we test the membership to the convex set conv(P ⊗ PPT BC ).The primal and dual versions of the corresponding estimation of visibilities χ FSEP P (ρ ABC ) are given as follows:

C.2 Full biseparability
Motivated by separability of bipartite states, one can consider three-qubit states that are separable with respect to any bipartition.Formally, this set can be expressed as FBSEP = ∩ α SEP(α| ᾱ) (fully biseparable).Being able to systematically approximate SEP(α| ᾱ) by the polytope approximation, we can also directly describe FBSEP.
In analogy to the already discussed cases of separability, FBSEP may be extended for general convex sets of operators P α = {σ α λ } N α λ=1 to FBSEP(P A , P B , P C ) defined by FBSEP(P A , P B , The membership of ρ ABC to FBSEP(P A , P B , P C ) is then equivalent to the existence of positivesemidefinite operators {τ ᾱ λ } N α λ=1 where τᾱ λ is an operator acting on ᾱ = {A, B, C} \ {α} such that The corresponding semidefinite programm to compute the visibility of a state ρ ABC to FBSEP is given as
Members of this set are called biseparable.Recalling that SEP(α| ᾱ) = conv(S α ⊗S ᾱ), let us now approximate the Bloch sphere S α by polytopes P α from inside (outside) for α = A, B, C. Then one has conv[∪ α conv(P α ⊗ S ᾱ)] as an inner (outer) approximation for BSEP.One can also easily see that these inner (outer) approximations of BSEP can be described by SDPs.Fixing three inner (outer) polytope approximations of N α vertices for the Bloch sphere, λ=1 , each respectively acting on the system ᾱ with α = A, B, C, such that which is an SDP.Explicitly, one can estimate the visibility of a state ρ ABC in the following way.

D Implementation of the polytope adaption technique
In this appendix, three different polytope adaption schemes for different situations of separability certification are explained.It will be discussed how the method can be implemented for the computation of the random robustness and absolute robustness in bipartite systems and the computation of the random robustness to full separability in multiparticle systems.

D.1 Bipartite random robustness
In the case of bipartite systems and mixing with the maximally mixed state, both system parts are iteratively replaced by polytopes that are obtained by the feasible optimum of the preceding semidefinite program given by ( 14).The method of polytope adaptions may thus be summarised by the following simple algorithm: 1. Initialise an arbitrary inner polytope P A 2. Compute χ with respect to P A , extract the corresponding τB i and construct a polytope P B with vertices {τ i /Tr(τ i )} 3. Exchange A and B, i.e. update ρ AB → Π AB ρ AB Π AB (Π AB : permutes systems A and B) and P B → P A , and go back to step 2 if convergence is not achieved.
Figure 5 illustrates how the number of needed iterations depends on the number of initialised vertices.

D.2 Bipartite absolute robustness
In the case in which one is interested in the robustness with respect to the whole separable set (see Section B), the method of polytope adaptions has to be slightly modified.Both of the lists of positive operators that one obtains as output from the first SDP (25)

D.3 Full separability in multiparticle systems
One situation in which polytope adaptions can also be applied is the treatment of full separability of a multiparticle quantum state.As described in the main text, we applied the method for states consisting of up to five qubits or up to three qutrits.We use a cyclic optimisation over all local subsystems.That means, in each iteration step all but one chosen local subsystem are approximated by a polytope.The feasible solution of this SDP in one iteration step then determines a new polytope for this chosen local subsystem.In the next iteration, another local system is chosen.In the case of systems involving qubits and qutrits, it is also possible optimise over two subsystems in one iteration by making use of the PPT criterion when it is necessary and sufficient for separability.

E Robust PPT-entangled states
The method of polytope approximations may be used to construct robust PPT-entangled states.
To do this, one has to initialise some PPT entangled state ρ AB in .A possible choice of such initial states are members of the families of Horodecki states ρ H 3|3 (a) and ρ H 2|4 (b) given by a 0 0 0 a 0 0 0 a 0 a 0 0 0 0 0 0 0 0 0 a 0 0 0 0 0 0 0 0 0 a 0 0 0 0 0 a 0 0 0 a 0 0 0 a 0 0 0 0 0 a 0 0 0 A reasonable choice of such a state is for example the 2 × 4 dimensional Horodecki state at its most entangled position, i.e. at b = 0.25.The strategy to construct robust PPT-entangled states involves then the following steps: 1. Set t = 0.
2. Perform a polytope adaption for ρ AB in , set t = χ(ρ AB in ) and save the resulting polytope P res on Alice's side.If the value of t did not change, exit.Minimas of the curve correspond to points that are found by the optimisation procedure.(b): A sketch of the algorithm to find local maximally robust fully biseparable states.Approximate witnesses to the set FSEP are constructed followed by a minimisation of the inner product of this operator among all fully biseparable states.
3. Perform the dual SDP (15) with respect to P res and save the resulting approximate entanglement witness Y AB .
4. Minimise the overlap Tr(Y AB ρ AB ) over all states ρ AB that have a positive partial transpose to obtain ρ AB out .
5. Update ρ AB in = ρ AB out and go back to step 2.
The resulting robust PPT-entangled state ρ AB res that we obtain has a visibility of χ(ρ AB res ) = 0.9461.In Fig. 6, the cross section spanned by the initial and the resulting states is plotted.It can be seen that the algorithm ended up in an unrelated PPT-entangled state lying in a different corner of the set of states.

F Robust fully biseparable states
To find and optimise states that are separable with respect to any bipartition but still not fully separable we invoke a see-saw algorithm which is very similar to the one used for robust PPTentangled states.We followed precisely the following steps to collect states which are locally maximally robust: 1. Set t = 0 and initialise a random three-qubit state ρ ABC in .
2. Perform a polytope adaption for ρ ABC in , set t = χ(ρ ABC in ) and save the resulting polytope P res on Alice's side.If the value of t did not change, exit.A state which results from this algorithm may be further transformed with local unitary operators for the minimisation of matrix entries which helps in the analysis of that state.The states that we obtain with this procedure take the form spanned by the vectors for certain angles θ .Notice that the index i here corresponds to (α, β) in the main text.Hence, each member in the family of states ρ(θ ) is proportional to a projection on a fourdimensional subspace that is spanned by A−BC product vectors.For them, we have one of the Bell states on the BC−system and some superposition parametrised by a coherence angle θ on the A−system, respectively.The state family shares some noticeable similarity with the post measurement state in the teleportation protocol using Bell-states [53].The only difference between them is the appearance of a factor i in |γ 2 (θ )〉 and γ 4 (θ ) .Physically, this difference is expressed by an extra phase flip operation on the subsystem A conditioned on measuring the Bell-states ψ + or ψ − .A comparison with the state without the respective phase flips reveals that the robustness is drastically increased.To see that ρ(θ ) is fully biseparable for all θ , it is useful to note that ρ(θ ) is "X-shaped" in the computational basis, i.e. we have with a = cos 2 (θ ), b = sin 2 (θ ) and c = sin(θ ) cos(θ ).By the necessary and sufficient criteria on bipartite separability of three qubit X-shaped states given in [54], it follows that ρ(θ ) is fully biseparable according to Eqns. ( 1) -(3) in Ref. [54] if and only if cos 2 (θ ) sin 2 (θ ) ≥ | sin(θ ) cos(θ )| which is trivially fulfilled by equality for all θ .On the other hand, it is still an interesting question to ask for which angles the state ρ(θ ) is most robust against full separability when mixing with white noise.There are eight points between 0 and 2π for which the visibility is minimal, see Fig. 7, and they are given by θ 1−4 =≈ 0.48+nπ/2 and θ 5−8 ≈ 1.09+nπ/2 (n = 0, 1, 2, 3).A reason for this exact numerical values may be seen when considering the reduced states of |γ i (θ )〉 on the A subsystem.Their Bloch vectors take the form so that their mutual differences are either given by ||r 1 − r 2 || 2 = 4(cos 4 (θ ) + sin 4 (θ )) or ||r 1 − r 3 || 2 = 16(cos 2 (θ ) sin 2 (θ )).A straightforward calculation reveals that all vectors have the same distance (and therefore span a regular tetrahedron inside the Bloch sphere, see Fig. 8, exactly when the equation cos 2 (θ ) sin 2 (θ ) = 1/6 is fulfilled, e.g. for θ = arccos 1 2 + 1 12 .Interestingly those are exactly the angles θ 1−8 for which the robustness of ρ(θ ) is maximised.The local maxima at e.g.θ ≈ 0.777 correspond to the situation in which the Bloch vectors lie in a two-dimensional plane, see Fig. 8.

SEPFigure 1 :
Figure 1: (a) Schematic sketch of the convex set of separable states.The visibility χ(ρ AB ) of a quantum state ρ AB corresponds to the fraction of separable states in the convex hull with the maximally mixed state, Eq. (4).Inner and outer polytopes approximating Alice's set of states (generalised Bloch sphere) give rise to lower and upper bounds of χ, respectively.(b) Sketch of the polytope approximation.Alice's generalised Bloch sphere (on the left) is replaced with a polytope while Bob's generalised Bloch sphere (right) is left unchanged.The parameter λ indicates the random variable correlating Alice's and Bob's states.Upon approximating Alice's states by a polytope, it takes values as vertices of the polytope.

Figure 5 :
Figure 5: Comparison of the number of needed iterations for different numbers of polytope vertices using random states in 2 × 2 -7 × 7 dimensional systems.The algorithm stops if the difference of calculated values within one round of iteration is below 10 −4 .The convergence of the visibilities is confirmed by comparing with upper bounds given by the PPT criterion.

Figure 6 :
Figure6: Cross-section of the 2-dimensional plane in the state space spanned by the most robust 2 × 4 Horodecki state ρ in and the obtained stronger entangled PPT state ρ out .Remarkably, the stronger PPT entangled state is located in a different region of the global state space in the sense that their convex hull is mainly included in SEP.

Figure 7 :
Figure 7: (a): Plot of the maximal mixing parameter of ρ(θ ) for different angles θ .Minimas of the curve correspond to points that are found by the optimisation procedure.(b): A sketch of the algorithm to find local maximally robust fully biseparable states.Approximate witnesses to the set FSEP are constructed followed by a minimisation of the inner product of this operator among all fully biseparable states.

Table 1 :
Comparison of lower bounds for χ computed with the method of adaptive polytopes and known upper bounds from the literature.The three qutrit W, GHZ and four qubit Cluster states are defined by the vectors W3 have to be used in the next iteration step.The list { ηB λ } must be invoked to construct a separable state and the list {τ B λ } for a new polytope approximation in every step of the iteration.1. Initialise an arbitrary inner polytope P A given by the vertices {σ λ } 2. Compute χP A with respect to P A by the SDP (25), extract the corresponding τB λ and ηB AB extracted from the previous SDP.4. Evaluate the corresponding σ λ , update P A with the vertices {σ λ /Tr(σ λ )} and go back to step 2 if convergence is not reached.
4. Minimize the overlap Tr(Y ABC ρ ABC ) over all states ρ ABC that are fully biseparable to obtain ρ ABC out .If the value is non-negative, go back to step 1 (thresholds for fully separable and fully biseparable coincide).5. Update ρ ABC in = ρ ABC out and go back to step 2.