Numerical construction of multipartite entanglement witnesses

Entanglement in multipartite systems is a key resource for quantum information and communication protocols, making its verification in complex systems a necessity. Because an exact calculation of arbitrary entanglement probes is impossible, we derive and implement a numerical method to construct multipartite witnesses to uncover entanglement in arbitrary systems. Our technique is based on a substantial generalization of the power iteration, an essential tool for computing eigenvalues, and it is a solver for the separability eigenvalue equations, enabling the general formulation of optimal entanglement witnesses. Beyond our rigorous derivation and direct implementation of this method, we also apply our approach to several examples of complexly quantum-correlated states and benchmark its general performance. Consequently, we provide an generally applicable numerical tool for the identification of multipartite entanglement.


I. INTRODUCTION
Quantum entanglement is one of the most fundamental concepts in physics.It was introduced in the pioneering works of Einstein et al. [1] and Schrödinger [2].The pure existence of this quantum phenomenon challenged previously established notions of correlations and paved the way towards a new interpretation of the nature of physics.Eventually, this led to new protocols used in quantum computing and communication, which utilize the resources of entangled quantum states [3].Examples of such classically infeasible tasks are quantum teleportation [4] and dense coding [5].Other early protocols concern quantum key distribution, known as BB84 [6] and E91 [7], and significantly improve communication security.Therefore, entanglement plays a key role in fundamental physics and technology-oriented applications.
A primary concern in the research of entanglement is the actual detection of this quantum correlation.Since a lot of protocols for quantum technologies rely on the presence of entanglement, the question whether or not an experimentally generated state is entangled has become a highly relevant topic.However, determining entanglement of general stateslikewise its counterpart, separability-is an NP-hard problem [8,9].
Another challenge specific to multipartite systems is the possibility that classical and quantum correlations can be differently distributed among the parties of an ensemble of systems.This leads to complex structures of multipartite entanglement; see, e.g., Refs.[10][11][12].Most notably, there are inequivalent forms of entanglement, which need to be distinguished.These are already present in systems of only three qubits, such as the prominent GHZ and W states [13].Beyond that, current experiments become more and more capable of producing large-scale entanglement [14][15][16].However, while entanglement is vital for characterizing such experiments, the tools to uncover highly quantum-correlated systems are rather limited, and the general verification remains an open problem.
Still, several criteria have been developed to successfully determine entanglement in bipartite and multipartite systems; see Refs.[17][18][19] for thorough lists of these entanglement tests.A prominent example is the partial transposition criterion [20], which has been generalized to general positive, but not completely positive maps [21].Furthermore, such maps are equivalent to entanglement witnesses [21][22][23][24].A crucial point of using witness operators is their nature of being observables, which can be directly implemented in experiments.Another main advantage is that no full quantum state reconstruction is required to apply such witnesses.Rather, a few measurements of the observable can be sufficient to experimentally uncover entanglement [25][26][27].
Consequently, witnesses have become a widely applied method for detecting entanglement.Their usefulness for quantum technologies has been shown to be promising by detecting entanglement of multipartite cluster states in theory and experiments; see, e.g., Refs.[28,29].Also, witnesses are not limited to specific systems; for example, they apply to trapped ions [30] as well as hybrid systems which correlate vastly different degrees of freedom [31].In addition, device-independent witnesses have been proposed for a robust verification of entanglement [32].For instance, such device-independent witnesses can be constructed via so-called matrix-product extensions [33].
An entanglement witness has a non-negative expectation value for separable states as it defines a hyperplane bisecting the set of states-one part containing at least all separable states and another part including exclusively entangled ones.In order to maximize the detectable range of entangled states, optimal witnesses have been introduced [34][35][36][37][38].A universally applicable approach is the method of separability eigenvalue equations (SEEs) which enables the construction of optimal witnesses in the bipartite and multipartite scenarios [39,40].The solution of the SEEs renders it possible, in principle, to formulate all entanglement witnesses.However, because of the general complexity of the separability problem, exact solutions are only known for specific scenarios.Still, this has already led to deeper insights into the complex forms of experimentally generated multipartite entanglement [41,42].
Once a witness-construction approach is realized, it can be applied to different physical systems and reveal more insight than the basic indication of entanglement.For example, entanglement in systems of indistinguishable particles can significantly differ from the case of distinguishable particles, but witnessing can be done in a similar manner [43][44][45].Furthermore, the quantification of entanglement can be based on witnesses as well [46][47][48][49].This also includes entanglement tests for the so-called Schmidt number in the bipartite systems [50][51][52], as well as its multipartite extension [12,53].
Since calculating witnesses is a hard problem and exact solutions are rare, a numerical approach is favorable.Numerical methods often use the convexity of the set of separable states.Prime examples are approaches based on semidefinite programming, used for the general, convex optimization of linear problems [54].The formulation of witnesses has the structure of exactly that kind of problem.Thus, semidefinite programming is a frequently applied method for probing entanglement [55][56][57][58][59][60].However, this approach addresses a general class of optimization tasks and is not specifically designed to address the properties of entangled systems.Consequently, such a general approach cannot present an optimal strategy to construct entanglement witnesses for arbitrary systems.Moreover, numerical standard approaches to solve the eigenvalue equations (EE), such as the well-known power iteration (PI) [61], do not apply to the construction of entanglement witnesses via the nonlinear SEEs.
In this contribution, we devise a numerical approach to construct multipartite entanglement witnesses by finding the maximal separability eigenvalue.Based on the properties of the SEEs, the analytical background is derived for our technique-termed the separability power iteration (SPI).As a special case, our approach includes the PI, which returns the maximal solution of EEs.We implement the SPI algorithm numerically.This is used to demonstrate that the directed design of our numerical approach is an efficient method compared to standard techniques applicable to arbitrary optimization problems.To outline possible applications, we use our algorithm, for example, to verify entanglement of weakly correlated, i.e., bound-entangled, states in the bipartite and multipartite scenarios.Therefore, an accessible algorithm is provided which renders it possible to construct entanglement probes for certifying multipartite quantum correlations.
We organize the paper as follows.Preliminary statements are made in Sec.II.Here, we introduce the framework used throughout the contribution and recollect information about entanglement.In Sec.III, the SPI algorithm to find the maximal separability eigenvalue of a positive operator is introduced.Proofs for the working behavior and the convergence of the algorithm are given.We analyze the performance of our algorithm in Sec.IV.In Sec.V, entanglement in a selection of bound-entangled states is analyzed.In Sec.VI, we discuss the connection between the SPI and experimental measurements as well as other entanglement criteria and show the broad applicability of our newly devised method to different problems.We conclude in Sec.VII, where we also summarize our re-sults.

II. PRELIMINARIES
In this section, we revisit multipartite entanglement and its verification.In particular, we concentrate on the previously introduced method of SEEs and its relation to standard EEs, which is essential for the following investigations.Eventually, we summarize these methods in the context of the considered problem which is solved by our numerical approach, the SPI.

A. Multipartite entanglement
Say S is the set of all pure states that are separable in an N-partite system.This means that the elements of S take a tensor-product form, where |a j ∈ H j is an arbitrary state in the jth subsystem and a j |a j = 1 for j = 1, . . ., N. Furthermore, a mixed state σ is separable by definition [62] if it can be written as where P is a classical probability distribution over S .Conversely, a state ρ is defined to be entangled if it cannot be expressed in this way.The given form of separability is also called full separability of an N-partite system.To consider instances of partial entanglement, we can assume that each of the N parties is itself a composition of K j subsystems.This allows us to study arbitrary forms of partial separability-e.g., N-separability-in a system which, in total, consists of It is also worth mentioning that continuous-variable entanglement can always be detected in finite-dimensional subspaces [63].Hence, we can restrict ourselves to Hilbert spaces with a finite dimensionality, d j = dim H j < ∞.

B. Entanglement witnesses
Based on the convexity of the set of separable states [cf.Eq. ( 2)], so-called entanglement witnesses, Ŵ , have been introduced [21][22][23].They fulfill the property that for all separable states σ , the inequality tr( σ Ŵ ) ≥ 0 holds true.Consequently, entanglement is detected if this inequality is violated, tr( ρ Ŵ ) < 0. In particular, it has been shown that witness operators can be written in the form [24,39] where g max is the maximal expectation value of L for separable states.Therefore, the following approach is equivalent to the method of witnessing [39,40]: For any entangled state ρ, there is a Hermitian operator L such that the entanglement of ρ is certified by the criterion tr( L ρ) > g max . (4) The other way around, a state σ is separable if for all L the inequality tr( L σ ) ≤ g max holds true.Moreover, it has been shown that it is sufficient to consider (normalized) positivedefinite operators only; see, e.g., Ref. [39].We refer to operators satisfying as positive operators in this work.To determine the bound g max , applied in the entanglement criterion (4), we introduce the SEEs [39,40] (see also Appendix E).

C. Separability eigenvalue equations
There are two equivalent forms of the SEEs [40].For this work, the more important representation of the SEE reads L|a 1 , . . ., a N = g|a 1 , . . ., a N + |χ .
Here, the vector |χ is N orthogonal to |a 1 , . . ., a N Namely, we have a 1 , . . ., a j−1 , x, a j+1 , . . ., a N |χ = 0 for all j = 1, . . ., N and for all |x ∈ H j .The normalized vector |a 1 , . . ., a N is the separability eigenvector.The real value g is the separability eigenvalue, which can also be written as the expectation value of L with respect to the separability eigenvector, The disturbance to a standard EE, created by the N-orthogonal vector |χ , couples the individual subsystems represented by the states |a j .Thereby, it creates a highly nonlinear equation which, in general, cannot be solved straightforwardly.Furthermore, we can relate the separability eigenvalues to our necessary and sufficient entanglement criterion given in inequality (4).Namely, we have [40] g max = max{g : g solves Eq. ( 6)}.
Let us stress that the maximal separability eigenvalue is the solution to an optimization problem that maximizes the function a 1 , . . ., a N | L|a 1 , . . ., a N for normalized, pure, and separable states.Moreover, using relation (7), the value of g max is determined through the corresponding separability eigenvector.Finding this vector |a 1 , . . ., a N is the goal of our algorithm to be introduced.Furthermore, the SEE in Eq. ( 6) takes the form of a perturbed EE.In fact, for a single party (N = 1), the vector |χ necessarily vanishes, which means that Eq. ( 6) corresponds to the EE.This relation between the SEE and the EE is relevant for our algorithm.
Furthermore, let us also recall properties of the SEEs, which are of particular importance for this work.First, the separability eigenvectors of the operator µ L + ν 1, for real numbers ν and µ = 0, are identical to those of the operator L [40].This allows us to restrict ourselves to positive operators, as mentioned above.
The second property to be discussed here addresses the relations between the operators L = |ξ ξ | and L = tr N ( L), (9) where |ξ = 0 is an arbitrary vector in the N-partite system and tr N denotes the partial trace over the Nth subsystem.This also implies that L is positive semidefinite and acting on an (N − 1)-partite system.The theorem of cascaded structures [40] states that the nonzero separability eigenvalues of L and L are identical, which also implies that Moreover, the separability eigenvectors of L and L read |a 1 , . . ., a N and |a 1 , . . ., a N−1 , respectively, where the Nth component obeys This means that |a N is parallel to a vector that is obtained from |ξ by projecting its first N − 1 components onto a j |.We emphasize that the optimization of the expectation value of the operator L over |a 1 , . . ., a N ∈ S , i.e., | ξ |a 1 , . . ., a N | 2 → max, corresponds to a maximization using L , which is defined in one subsystem less than used for L. Also recall that the operator L is, in general, not a rank-1 operator anymore, and the cascaded structure is applicable only to rank-1 operators.

D. Preliminary discussion
In Fig. 1, we outline the previously discussed entanglement detection method using three different operators, labeled as L, L , and L .The tangent hyperplanes separate the set of separable states from states that are verified to be entangled.The touching points of the tangent represent the separability eigenvectors to the maximal separability eigenvalue.In general, the more operators are used, the better the hyperplanes can approximate the bounds of the set of separable states and the more entangled states can be identified.Note that one can construct a dense set of operators for such an approximation with arbitrarily high precision; see, e.g., Ref. [39].
Both the construction of multipartite entanglement witnesses and the approximation of the set of separable states depend on the solution of the SEEs.Specifically, we need to find the maximal separability eigenvalue, which is determined through its corresponding separability eigenvector.However, the SEEs present a sophisticated mathematical problem, which has at least the complexity of the standard eigenvalue problem [64].In fact, independently of our specific approach, the separability problem has been shown to be an NP-hard problem [8,9].
Furthermore, the SEE in Eq. ( 6) shares a number of properties with the EE, L|z = g|z .For the latter EE, there exists an algorithm to compute the eigenvector to the maximal eigenvalue of any positive operator L, the PI [61].In this algorithm, a vector |z is mapped onto a new normalized vector, |z = L|z / z| L2 |z 1/2 .An s-step iteration, |z , |z , |z , . .., |z (s) , yields a vector that approaches, for s → ∞, an eigenvector to the maximal eigenvalue of L for any initial vector that is not already an eigenvector to L.
In the following, we aim to generalize the PI to be applicable to the SEE.For this reason, we introduce an algorithm for a numerical implementation, which yields the desired solution of the SEEs-a separability eigenvector to the maximal separability eigenvalue.The resulting SPI algorithm is applicable to all positive operators L and enables the construction of witnesses to probe multipartite entanglement.

III. THE SPI ALGORITHM
In this section, we present the SPI algorithm-step by step.The flowchart of this algorithm to construct entanglement criteria is shown in Fig. 2. Our approach yields the separability eigenvector |a 1 , . . ., a N to the desired, maximal separability eigenvalue for a positive operator L [Eq. ( 7)].Before we study the individual, essential parts of the SPI in a rigorous mathematical framework, let us first get a general overview of how our algorithm operates by applying it to an example.

A. Proof of concept
For demonstrating the function of our algorithm, we consider the bipartite (N = 2) and positive operator FIG. 1. Visualization of three entanglement criteria.Entanglement is verified in the shaded half-spaces tr( L ρ) > g max (bottom, yellow area), tr( L ρ) > g max (right, cyan area), and tr( L ρ) > g max (left, magenta area), where the values g max , g max , and g max are the maximal separability eigenvalues of the operators L, L , and L respectively.The boundaries define hyperplanes tangent to the set of separable states (gray area).The bullet points correspond to the separability eigenvector to the maximal separability eigenvalue for each operator.
where V is the swap operator implies that the maximal separability eigenvalue is g max = 2, and it is attained for |a 1 ⊥ |a 2 [39].This exact result serves as our reference to assess the success of our algorithm for this example.Moreover, since the maximal standard eigenvalue is three, it follows from g max < 3 that this operator can be used to detect entanglement [52].In fact, the swap operator is related to the prominent partial transposition criterion to verify entanglement [20,21,39].
Our algorithm in Fig. 2 is initialized at point 1 with the operator (12) and the number of subsystems being N = 2.At 2 , let us begin with states |a 1 , a 2 , which are neither parallel nor orthogonal, to exclude the trivial cases.Namely, we have 0 < |γ| 2 < 1, where Say that in step 3 , we do not have convergence yet; i.e., we follow the branch labeled "false" and compute the vector in step 4 , Since N = 1 (step 5 ), we proceed to 6 and compute the operator, which is a single-subsystem operator.This step is referred to as forward iteration in Fig. 2. The idea behind this step is a result of the theorem of cascaded structures, which finds the maximal separable projection onto the state |Ψ ; see Sec.II C.This also allows us to apply the SPI to L + 1.As L is positive semidefinite by construction, the addition of 1 assures the positivity of L + 1 without modifying the separability eigenvectors.
Calling the SPI with N → N −1 = 1 in 7 and thereby going back to step 1 and going through steps 2 to 5 , we follow the branch for which N = 1 is true.This gives an iteration of steps 10 , 11 , 3 , 4 , and 5 , indicated through the dashed box in Fig. 2, which describes the PI.The PI is employed for solving the standard EE numerically by returning the eigenvector to the maximal eigenvalue of a positive operator with an arbitrarily high precision.So we can assume that the convergence 3 is true after some iterations of the PI.For the given operator L + 1 (thus, also for L ), the eigenvector to the maximal eigenvalue reads using γ [cf.Eq. ( 13)], the abbreviation and the normalization ν = [2Γ(Γ − 3 + 4|γ| 2 )] 1/2 .Thus, the PI basically returns the vector |a 1 in step 12 , which is used to continue with the case N = 2, where we exit 2. Flowchart of the SPI algorithm.Branches in the algorithm, either "while" loops or "if" conditions, are represented by magenta diamonds.Yellow rectangles represent assignments and function calls.Entry and exit points of the algorithm are shown as cyan ellipses.Box 1 refers to the input, which includes an operator and the number N of parties.A vector is generated in 2 to serve as our starting vector.If the convergence criterion 3 , studied in Sec.III B 4, is not met, we generate another N-partite vector in 4 .For N = 1 in 5 , the algorithm corresponds to the power iteration (PI) which finds the standard eigenvector to the maximal eigenvalue (dashed box).For N > 1, our extension consists of three essential parts (gray areas), which are the forward iteration 6 and the backward iteration 8 as well as a recursion calling the SPI for N − 1 parties in 7 .Eventually, when 3 is satisfied, the desired separability eigenvector is the output of our SPI algorithm and returned in 12 .
step 7 to perform the backward iteration step 8 .This gives a vector in the second subsystem.For convenience, this vector is renamed (step 9 ) and normalized (step 10 ); see Fig. 2. Again, the backward iteration is a result of the theorem of cascaded structures, which relates the separability eigenvectors of L for N − 1 subsystems with those of the initial operator L for N, cf.Sec.II C.This yields the state of the second subsystem, Thus, we obtain a new separable state |a 1 , a 2 in 11 , where the tensor-product state is formed.
What did we achieve with the construction of this new state?To answer this question, let us recall that the desired separability eigenvector of the operator under study has perpendicular components for the subsystems.Thus, in analogy to Eq. ( 13), we may compute the scalar product of the states of the subsystems, which yields γ = γ/Γ and This means the states |a 1 and |a 2 are closer to orthogonal than the initial states |a 1 and |a 2 .Equivalently, we can say that the expectation value of L increases, for which holds.Therefore, we get a convergent sequence of separability eigenvectors which, in the limit of infinite iterations, yields the desired exact maximal separability eigenvalue, a → 2 − 0 = g max for s → ∞.In conclusion of this example resulting in an entanglement test based on the swap operator [cf.Eqs. ( 12) and ( 4)], our SPI is constructed to deliver the separability eigenvector to the maximal separability eigenvalue.Applying properties of the theorem of the cascaded structure of SEEs, we identify the following essential steps: forward and backward iteration.The forward iteration allows the reduction of the number of subsystems by one in each recursion depth until the recursion depth reaches a maximum when the operator is a single-partite operator.Then the SEE reduces to the EE, and the PI is used to get the eigenvector to the maximal eigenvalue.The eigenvector is further used in the next step, the backward iteration, to obtain the remaining subsystem components of the separable product vector.After completing multiple instances of such a cycle, we obtain an arbitrarily precise approximation to our sought-after separability eigenvector.Now, we may consider the general case beyond the specific example, which was used to demonstrate the general operation of our generally applicable algorithm in Fig. 2.This gives the mathematically rigorous formulation of the SPI for arbitrary positive operators L and arbitrary numbers N of subsystems, which necessarily requires a rather technical treatment because of the complexity of the underlying separability problem.After this, we perform a benchmarking of our algorithm and apply it to various examples, which provides a more intuitive assessment of our method.

B. Analytic framework
Based on the theorem on cascaded structures for the SEEs, the SPI iterates over the number of parties from N to one.For N = 1, the SPI and PI are identical, resembling the underlying fact that the SEE and EE are the same in this case too.Beyond the PI, the SPI algorithm includes two main steps, denoted as forward and backward iteration.Clearly, the major goal of our maximization algorithm for a positive operator L is to get a new separable state |a 1 , . . ., a N from the preceding state |a 1 , . . ., a N , which increases the expectation value, a 1 , . . ., a N | L|a 1 , . . ., a N > a 1 , . . ., a N | L|a 1 , . . ., a N .Here, let us discuss the details, the proofs of some of the required theorems are provided in the corresponding appendixes.

Initial considerations
Let us make some more general observations, which we then apply to the separability problem under study.A positive operator L induces a scalar product, for arbitrary |x and |y .Therefore, the Cauchy-Schwarz inequality holds true, | x|y L| 2 ≤ x|x L y|y L, where the equality is equivalent to |x |y .Also, we have x|x L > 0 for all |x = 0. To apply these features, we have to restrict ourselves to positive operators L. Note that in our following proofs, we rely on the properties of the scalar product; for example, a positive-semidefinite operator L would be insufficient [65].Say T is a closed and bounded subset of a finitedimensional Hilbert space.For a |z ∈ T , one can define an iterated state as where we use the function "arg max," which returns the argument for which the maximum is reached.In other words, Since |z is also an element of the set T over which we maximize, we can conclude that z|z L ≤ | z |z L|.Applying the Cauchy-Schwarz inequality, we get Considering the second and fourth terms, as well as the first and third terms, we find the increasing sequence From the definition of |z and the properties of the Cauchy-Schwarz inequality, we can also conclude that the equality holds true if and only if |z |z .Therefore, we can state that the iteration |z , |z , |z , etc. produces a sequence of increasing expectation values, However, the elusive arg max function (22) has to be computed for this purpose.
In fact, this can be done for separable states, T = S .We may use the abbreviation |Ψ = L|a 1 , . . ., a N .To maximize the projections of this state onto separable ones, we can apply the theorem of the cascaded structure.This means that the maximal projection of this state onto separable states is obtained by |a 1 , . . ., a N for the maximal separability eigenvalue of |Ψ Ψ|.In Sec.II C and in the flowchart of the SPI in Fig. 2, we describe how this is achieved: We reduce the number of parties N and solve the SEE for L = tr N |Ψ Ψ| (forward iteration, 6 ) to get |a 1 , . . ., a N−1 (step 7 ), which then determines the remaining component |a N from a 1 , . . ., a N−1 , • |Ψ (backward iteration, 8 ).
In summary, the cascaded structure describes how to compute the desired arg max function for separable states.This describes the underlying principle of the SPI, which allows us to compute the bounds g max for the necessary and sufficient entanglement criteria (4).

The SPI
To apply the general relations above, let us begin with the forward iteration step.By the following Theorem 1, it is guaranteed that finding the separability eigenvector corresponds to determining the maximal separability eigenvalue for the (N − 1)-partite case.More specifically, it enables us to reduce the number of subsystems for the SEE by one.This theorem is a direct consequence of the SEE in Eq. ( 6) and its properties.In the SPI algorithm, the theorem is applied in the forward iteration step 6 .To find the full N-partite separability eigenvector, a reverse step has to be taken.Theorem 2 states how the N subsystem separability eigenvector can be generated from the N − 1 subsystem separability eigenvector.The proof of this theorem directly follows from the cascaded structure, cf.Eq. (11).It relates the separability eigenvector for N parties to those of a lower number of parties, N − 1. Thereby, if a solution to the (N − 1)-partite SEEs for L is known, we directly find the Nth component of the full solution |a 1 , . . ., a N .In the flowchart in Fig. 2, we see the application in the backward iteration step 8 .
The combination of Theorems 1 and 2 is fundamental for the SPI to work.In fact, one might visualize the working principle of the algorithm as a nested cascading structure.The forward iteration is recursively applied until we reach the case N = 1.In that case, the standard PI is performed.After that, the backward iteration finalizes the individual recursion layers of the SPI until we obtain the new N-partite separable vector.Then, we can start a new cycle of forward iterations, the PI, and backward iterations until the convergence is reached.The algorithm will terminate successfully and return the complete separability eigenvector corresponding to the maximal separability eigenvalue g max for detecting entanglement in terms of inequality (4).
To verify the statement that the algorithm converges to the maximal separability eigenvalue, a few observations have to be shown first.Let us take a closer look at the sequence of product vectors created by the SPI.In every step, we find an element of all product states, |a 1 , . . ., a N ∈ S , which projects maximally onto the action of operator L onto the previously generated product state |a 1 , . . ., a N .This iteration is done until we reach convergence.This generates a monotonously growing sequence of expectation values of L, which is stated in the following theorem.This theorem is a special case of the general considerations made in Sec.III B 1. In addition, the global phase of the separable state |a 1 , . . ., a N can be chosen freely, which we conveniently select such that we have positive projections onto |Ψ , i.e., a 1 , . . ., a N | L|a 1 , . . ., a N = | a 1 , . . ., a N | L|a 1 , . . ., a N |.Let us stress that the arg max function, i.e., finding the maximal projection onto |Ψ = L|a 1 , . . ., a N , is obtained from the cascaded structure; see also Theorems 1 and 2.
Because of Theorem 3, the SPI produces a sequence of increasing expectation values of L. This observation is an important aspect for the proof of convergence of the SPI, which is shown in two parts.Both theorems rely on the sequence (g (s) ) s of expectation values generated by the SPI in each step s, where Here, in analogy to the example in Sec.III A, the vector |a is the approximation to the separability eigenvector for the maximal separability eigenvalue after s iterations of the SPI.First, we consider the local convergence of the algorithm.
Theorem 4 (Local convergence).For any starting vector, the sequence (g (s) ) s of expectation values generated by the SPI converges, i.e., the limit exists and is bounded as 0 ≤ ḡ ≤ g max .See Appendix C for the proof.
For an arbitrary starting vector, a sequence of expectation values of L for separable states is generated.The generated sequence converges independent of the choice of starting vector.Combining the statements from Theorems 3 and 4, we conclude that there is a monotone growth of expectation values towards a maximum.This maximum does not necessarily need to be the maximal separability eigenvalue g max as shown in Theorem 4. We therefore require an additional observation to prove global convergence of the SPI, which is stated in Theorem 5.
Theorem 5 (Global convergence).Let Σ be a set of separable starting vectors and (g (s) Φ ) s be sequences of expectation values generated by the SPI for a starting vector |Φ ∈ Σ.Furthermore, say { ḡΦ ∈ R : |Φ ∈ Σ and ḡΦ = lim s→∞ g (s) Φ } defines the set of optimal expectation values (limits of the converged sequences) for each starting vector.The maximal separability eigenvalue for the operator L is g max = max Φ∈Σ { ḡΦ }, which is the maximum of the limit to the series of expectation values for each starting vector.See Appendix D for the proof.
The set Σ of different starting vectors |Φ that we consider is covered in Sec.III B 3.Even in the worst-case scenario, it is far smaller than the set S of all separable states.

Starting vectors
An important aspect for the implementation of the algorithm is the choice of a starting vector, cf.Theorem 5.Because a proper choice can significantly decrease the runtime of the algorithm, let us provide more details on this aspect.
Assume we start with a separability eigenvector corresponding to any-except for the largest-separability eigenvalue.Then, the algorithm converges immediately, and the resulting separability eigenvalue will not be maximal.It is worth mentioning that such a behavior is already well known for the PI.It is straightforward to check whether an initial vector is a (separability) eigenvector.Similar results might happen for starting vectors that are too close to any (separability) eigenvector.
To circumvent such problems, the SPI can be run multiple times with different starting vectors, chosen as an operator basis.Namely, the set {|a 1,k , . . ., a N,k a 1,k , . . ., a N,k |} k of starting vectors spans all operators of the underlying Hilbert space.This allows us to cover all parts of the operator space and, of course, also resolves the related problem for the PI.
This choice is valid as the set of separability eigenvectors can be used to find a decomposition of any state, similarly to the spectral decomposition found by regular eigenvectors.In fact, any positive-semidefinite operator can be decomposed in terms of projectors of separability eigenvectors; see Refs.[66,67] for proofs of the bipartite and multipartite cases, respectively.Specifically, the decomposition of at least one vector of the operator basis needs to contain the sought-after separability eigenvector.This warrants the choice of using the operator basis as starting vectors.
Another efficient ad hoc ansatz that we used for the implementation of the SPI is described as follows: First, a preliminary run of the SPI is done to find a product vector projecting maximally onto the vector L|v , where |v is the eigenvector to the maximal standard eigenvalue of L. This choice is inspired by the fact that the wanted vector |a 1 , . . ., a N maximizes the expectation value of L with respect to product vectors.The eigenvector |v maximizes the expectation value of L without the restriction to separable states.Second, the product vector that lies maximally parallel to |v serves as our starting vector.Finding such a maximal projection is in fact exactly what we get when running the SPI for a positive operator 1 + |v v|.Finally, the resulting product vector serves as the initial vector for the SPI algorithm applied to L.
Our numerical results and comparison with other methods confirm the assumption that the constructed starting vector is sufficient, as the described procedure returns the same values.Still, a rigorous proof of this observation requires further investigations.Until then, the choice of an operator basis of starting vectors is preferable in the general case.

Convergence criterion
The flowchart in Fig. 2 requires a check for convergence in step 3 .Theorems 4 and 5 guarantee, in theory, the convergence of the SPI.In a practical implementation of the algorithm, however, the computer needs to know when convergence is reached in a numerical sense.
We apply a convergence criterion that is based on the SEE.In the sth cycle of the algorithm, we obtain the vector |χ (s) = ( L − g (s) 1)|a N ; see Eq. ( 6).By definition (see Sec. II), |a N is a separability eigenvector if and only if |χ (s) is N orthogonal.Likewise, convergence is reached if and only if |χ (s) is N orthogonal.Theorem 4 guarantees that we approach this scenario-meaning that lim s→∞ a In fact, this N-orthogonality requirement can be used to quantify the closeness to the solution.For this purpose, we can evaluate if the following inequality is satisfied: (s) |<ε, for a sufficiently small ε and all x ∈ B j , where B j is a basis of H j .The machine precision of the representation of numbers on a computer bounds the value of ε.When the inequality is satisfied, the possible numerical convergence is achieved and the current iteration |a N is the desired approximation to the separability eigenvector.

IV. BENCHMARK
We now want to find out how the SPI performs as opposed to other methods that allow the construction of arbitrary entanglement witnesses.A simple brute-force approach to obtain the entanglement criterion (4) for L is to find all separable pure states and calculate the expectation value of L.Then, g max is the maximum of these values.
As the search space for these vectors is over a continuum, one could use a generally applicable global optimization algorithm, such as genetic algorithms [68].This presents a stateof-the-art method to solve optimization problems.It is rather fast and inspired by evolutionary processes in biology.Thus, we implemented such a genetic algorithm to evaluate the performance of the SPI.A genetic algorithm requires a fitness function to be minimized, which will be f (v) = − v| L|v , the negative of the expectation value of L. An intermediate step ensures that the argument vector |v is indeed a product vector.
During the runtime, the genetic algorithm will minimize f (v) and converge towards a vector |v 0 with f The resulting minimization will give the maximal separability eigenvalue g max = − f (v 0 ), or at least a close approximation.
To show the advantages of the proposed algorithm, SPI, as opposed to this simple maximization strategy, we compare the two approaches for the following, different scenarios: We consider a bipartite system (N = 2) and vary the dimensions, and increase the number of parties N. As discussed previously, we choose the convergence criterion in Sec.III B 4. The starting vector is chosen as the maximal separable projection on the (standard) eigenvector corresponding to the maximal eigenvalue of L.
To exclude any bias, the chosen operators are randomly generated by first defining a random operator M acting on the D-dimensional space, where D = d 1 • . . .• d N .Then, we construct a positive and normalized operator L for which we want to find the maximal separability eigenvalue as The SPI and brute-force approaches have been tested for 100 randomly selected operators.To make the runtimes comparable, the same set of random operators was used for both approaches.
Figure 3 shows the average runtime for the SPI compared to the brute-force approach.These results come from running both algorithms on a desktop computer.The runtime of the SPI is, on average, at least two orders of magnitude lower for the considered sample size of 100 randomly generated test operators.In bipartite systems (Fig. 3, top panel), we see a smaller scaling behavior of the SPI, whereas the scaling is about the same for an increasing number of qubits (bottom plot).Moreover, focusing on the numbers of subsystems (Fig. 3, bottom panel), we see that the SPI finds the maximal separability eigenvalue for a state acting on a 13-fold Hilbert space (dimensionality D = 2 13 = 8 192) in roughly the same time as the other approach manages to find in the ninefold case (dimensionality D = 2 9 = 512).It is also worth mentioning that all curves of the presented study in Fig. 3 can be roughly approximated by exponential functions of the overall dimensionality (D = d 2 [top] and D = 2 N [bottom]), representing the expected exponentially increasing runtime of the separability problem.The dip (in favor of the SPI) at N = 9 in the bottom plot cannot be explained at this point and requires further investigations.
Our benchmark indicates the superior potential of the SPI algorithm to numerically construct entanglement tests.Specifically, it outperforms the competing approach for highdimensional scenarios, which includes the dimensionality of the individual parties as well as the number of parties itself.This enables a comparably efficient tool for the identification of entanglement in complex physical systems.Keep in mind that the runtimes shown in Fig. 3 are from running the SPI on a desktop computer; computation clusters might improve the performance even further by a large margin.
Beyond the genetic algorithm, there exist more specialized algorithms, treating Eq. ( 7) as a maximization of a multivariate polynomial.Such approaches are also NP-hard problems, meaning they can not be solved in polynomial time by a nondeterministic Turing machine, and only lower bounds of the global maximum can be found in polynomial time [69].We apply one state-of-the-art realization of such an algorithm to find the maximum of a polynomial [70], using semidefinite programming, instead of the problem of finding the maximal separability eigenvalue of an operator.Semidefinite programming is a frequently applied technique used for entanglement tests, cf., e.g., Refs.[54][55][56][57][58][59][60].Already in a 3 × 3 case, the algorithm in Ref. [70] failed to be conclusive and, in fact, returned a lower value than our SPI.For use as an optimal witness, the true maximal separability eigenvalue is crucial; thus, the result of the competing algorithm could lead to a false indication of entanglement.In all other tested cases in which the algorithm was conclusive, our SPI was superior in terms of speed and accuracy.

V. EXAMPLES
As a proof of principle, let us apply our algorithm to detect entanglement of states of special relevance.Specifically, we study the two-qutrit Horodecki state [71] and the fourqubit Smolin state [72].Both states have been classified as bound-entangled states.In the case of the Horodecki state, this arises from the dimensions of the state, which is acting on a 3 × 3-dimensional Hilbert space.The Smolin state acts on a 2 × 2 × 2 × 2-dimensional Hilbert space, and the boundentangled nature arises from the fact that the state is separable with respect to all bipartitions consisting of two subsystems each, yet still entangled in all other partitions.By applying the SPI algorithm, we aim at confirming the weak entanglement properties of those bound entangled states for which the well-known partial transposition test [20,73] fails to be conclusive.
For our entanglement analysis based on the criterion (4), a positive-definite, Hermitian operator L is required.For simplicity, the test operator will be chosen as Lβ = ρβ .We calculate the maximal separability eigenvalues g β for every Lβ where 0 ≤ β ≤ 5.In this entanglement test, a state is verified to be entangled if which corresponds to the criterion based on the entanglement witness Ŵβ = g β 1 − Lβ .
The results are shown in Fig. 4. The entanglement criterion (33) is satisfied in the magenta colored areas.The blank area corresponds to parameters 2 ≤ α ≤ 3 for which no entanglement could be detected, which agrees with the prediction in Ref. [71].In all other cases (cyan area), there exists at least one other value β for which Lβ verifies entanglement.Thus, we correctly and straightforwardly certify entanglement of all  Horodecki states, which are positive under partial transposition, using our SPI approach.Beyond the bipartite case, let us apply our method to the multipartite scenario for which the partial transposition criterion does not apply in principle.For this reason, we study the four-partite Smolin state [72], where σx , σy , and σz denote the Pauli spin matrices.We restrict ourselves to a test operator of the simple form L = Ŝ.
The SPI algorithm was run for all partitions.The results are listed in Table I.For applying entanglement criterion Eq. ( 4), we additionally compute tr( L Ŝ) = 1/4.Thus, in agreement with the results in Ref. [72], entanglement could be verified for all partitions, except for the bipartition which consists of two subsystems each, i.e., {1, 2}:{3, 4}.
In this section, we demonstrated the direct application of our SPI algorithm to construct entanglement probes, for example, to identify bound instances of entanglement.We deliberately chose such weakly entangled states, which have been characterized previously to challenge our method and compare our numerical results with sophisticated exact analysis.In particular, entanglement was verified in bipartite qudit and multipartite qubit states.The entanglement of the states under study is a challenge for other directly applicable methods as the partial transposition criterion gives inconclusive results.

VI. DISCUSSION
We introduced, implemented, and applied a method to numerically construct entanglement tests.In this section, let us discuss how this technique can be used in experiments, how it improves other entanglement probes, and how it can be generalized to detect other forms of entanglement.Finally, we discuss future research directions that become accessible with our approach and address the interdisciplinary importance of the introduced technique by relating it to a current problem in pure mathematics.

A. Experimental implementation
A major benefit of our approach is the direct applicability in experiments.Suppose that the set of observables { Mk : k = 1, . . ., m} describes a measurement scheme.In other words, the data yield the expectation values Mk = tr( Mk ρ).An example for such operators relates to a displaced photon-number correlation [74].In general, a family of positive operators L can be constructed from the considered measurements, by choosing real-valued coefficients µ k and adjusting ν to ensure positivity of L. The entanglement criterion (4) can be applied.On the one hand, the experimental expectation value is given by L = ν + ∑ m k=1 µ k Mk .On the other hand, we get the maximal expectation value for separable states, g max , from the application of our SPI to the family of operators L under study.Note that a variation over the coefficients µ k also enables an optimal entanglement verification based on the set of measured observables, similarly to the technique applied to Gaussian measurements in Refs.[41,42].

B. Relations to other entanglement criteria
As mentioned earlier, our entanglement criteria are identical to witnesses [Eq.(3)].Furthermore, based on the Choi-Jamiołkowski isomorphism [75,76], entanglement witnesses enable the formulation of positive, but not completely positive maps to probe entanglement [21,22].Thus, our numerical method can be used to construct previously unknown families of such maps.For instance, the test operators that verified the entanglement of the bound-entangled states (Sec.V) necessarily lead to maps that go beyond the partial transposition since the partial transposition cannot detect the entanglement of states considered in those examples.
In addition, in Ref. [33], an elegant approach was formulated that enables the construction of device-independent entanglement witnesses from device-dependent ones.This technique is based on a matrix-product extension that assigns to each subsystem an auxiliary Hilbert space, but requires the previous knowledge of a witness.Such desired initial witnesses can be provided by our algorithm and combined with the method from Ref. [33] to construct device-independent entanglement witnesses.

C. Outlook
Beyond the witnessing of multipartite entanglement, the SEE approach has been generalized.Thus, let us briefly discuss some future generalizations of our numerical method for the aim of exploring entanglement in a broader context.
The detection of K-entanglement, and thus of genuine entanglement, is possible by finding the maximum of all maximal separability eigenvalues for an operator with respect to partitions of the length K [42].It is therefore a straightforward extension to the SPI to find the optimal witness for Kentanglement with the introduced algorithm-the algorithm is run multiple times for different partitions and the maximum of the results is the required separability eigenvalue.
Furthermore, some physical problems require solutions of a generalized EE, L|Φ = λ P|Φ , where the right-hand side includes a contribution that is different from the identity, P = 1.Interestingly, the same holds true for the SEE.
One example is the verification of entanglement in systems of indistinguishable particles, which is based on a generalized SEE and where P represents the (anti)symmetrization operator for bosons (fermions) [45].Another example is the quantification of multipartite entanglement via generalized Schmidt-number witnesses [12].There, P takes the form of a spinor projection (details can be found in the supplement to Ref. [12]).A third example is the detection of multipartite entanglement in systems for which the number of subsystems is not fixed.For instance, the underlying generalized SEE applies to the construction of multiparticle-entanglement witnesses for fluctuating particle numbers [77].
Thus, a generalization of the SPI to account for such generalized SEEs, including P, will further enhance the range of applications.It is worth mentioning that the desired generalization is well known for the PI, which is likely to be applicable to the SPI in a similar manner.
Furthermore, the standard EE applied to the density operator leads to the spectral decomposition of the state.Similarly, the SEE can be used to expand the density operator in terms of separability eigenvectors and a quasiprobability distribution [66,67].The latter one includes negativities iff the state is entangled; see Ref. [78] for an application to uncover bound entanglement.However, this approach requires the computation of all separability eigenvectors.Therefore, similar to the subspace iteration for the PI, a generalization of the SPI to include all solutions, beyond the one that corresponds to the maximal separability eigenvalue, could lead to a broader applicability of entanglement quasiprobabilities.
As a final example let us consider the dynamics of quantum systems, which is described by the Schrödinger equation.To distinguish the entanglement-generating evolution from the separable dynamics, we recently introduced the separability Schrödinger equations [79], which relate to the SEE in the static case.Again, the SPI can be the starting point for the numerical implementation of this approach.
Thus, generalizations of the SPI have the potential to uncover multipartite entanglement in a much broader sense.Beyond the already-available construction of positive, but not completely positive maps and device-independent entanglement witnesses, our numerical approach builds the foundation for the future studies of entanglement.

D. Relations to mathematical problems
The question of positive polynomials is an interesting and, in the most general case, unsolved mathematical problem, which has been studied for a long time [80] and finds many applications [81].As already indicated in Sec.IV, any entanglement witness can be characterized by the non-negativity of a multivariate polynomial [69].All entanglement witnesses can be generated through the solution of the SEE.Therefore, the solution of the SEE enables the construction and characterization of positive multivariate polynomials; see also Appendix E. Consequently, the proposed SPI is an alternative approach to numerically solving the positivity problem of polynomials.
Another family of important problems in pure mathematics that could benefit from the SPI are partial differential equations, which are also closely related to many problems in physics.For instance, the applicability of the method of separation of variables corresponds to the question of whether or not solutions are factorizable, i.e., a tensor product.Since a separable eigenfunction is also a separability eigenfunction [39], i.e., eigenvector in the function space, the SPI can be applied to find factorizable solutions of the partial differential equation.
Moreover, nonlinear partial differential equations address questions such as finding the ground state to a nonlinear energy functional.If this functional is polynomial, a problem related to the previously mentioned characterization of multivariate polynomials can be formulated.Namely, the numerical approximation to the ground state can be obtained by the multipartite SPI as the maximum of the negative nonlinear energy functional, resulting in the minimal energy.

VII. CONCLUSION
In this paper, we introduce an algorithm, the SPI, to numerically construct arbitrary multipartite entanglement witnesses.This algorithm enables us to find the maximal separability eigenvalues, which directly results in measurable entanglement tests.Beyond the formulation of our method, we also provide the mathematical background for the SPI, which yields the maximal solution of the nonlinear separability eigenvalue problem addressing the complex entanglement problem in quantum physics.Furthermore, our framework is supplemented by performing a benchmark of our approach, applying it to uncover hard-to-detect forms of entanglement, and relating it to other methods in the theory of quantum entanglement and their experimental application.
Our algorithm shows two crucial steps-namely, forward and backward iteration-following directly from the cascaded structure of the separability eigenvalue equations.The forward iteration reduces the number of parties until we have a single-party problem, which is then used in the backward iteration to solve the multipartite problem.This property also allows us to prove the convergence of the SPI to reliably produce entanglement tests based on arbitrary observables.Interestingly, our algorithm includes the well-known power iteration, which is able to calculate the maximal (standard) eigenvalue, as a special case.
We show the efficiency of our approach in comparison with another method, which is mainly based on a genetic algorithm.The genetic algorithms presents a state-of-the-art approach to solve arbitrary optimization problems.The SPI is faster by two orders of magnitude, which is partly because of its directed design to specifically address the entanglement problem.For example, we analyze the runtime as a function of the dimension of a bipartite quantum system.In addition, we numerically solve the separability eigenvalue equations in a feasible time for operators up to a 13-party qubit Hilbert space, corresponding to 8 192 dimensions.
Furthermore, we apply the SPI to bound-entangled states whose entanglement detection is a cumbersome problem.For instance, the frequently applied partial transposition criterion fails to uncover the entanglement of the considered examples.Applying the SPI, we straightforwardly verify this weak form of entanglement, proving the advantage of our method.Moreover, we demonstrate with these examples that our algorithm renders it possible to uncover entanglement of all forms of partial entanglement in multipartite systems.It is also worth mentioning that entanglement of continuous-variable systems can be detected in finite subspaces, allowing us to apply our algorithm to these kind of states as well.
We outline the versatile nature of our method and its impact on future research by relating it to other open problems in quantum entanglement and beyond.For instance, the construction of entanglement witnesses, which is achieved by our SPI, is the basis for the formulation of positive, but not completely positive maps for entanglement detection and the construction of device-independent entanglement witnesses.Furthermore, we describe the construction of entanglement criteria based on measured quantities and outline several generalizations, which are-at their core-related to our method.
Thus, we devise a relatively simple, yet versatile approach to numerically construct entanglement tests in multipartite systems.The direct implementation of our method enables us to certify complex forms of quantum correlations based on measurable criteria.In addition, we derive the required mathematical background of our algorithm to ensure its operation and benchmark its performance.To the best of our knowledge, there exists no alternative method of entanglement verification that is applicable to complex systems that our method can manage.To summarize, we provide a full numerical framework for the detection of multipartite entanglement for theoretical studies and, more importantly, for application in current and future experiments using entanglement in quantum information and communication protocols.
As the Rayleigh quotient is invariant under the norm of the vector, we may assume a j |a j = 1.Consequently, the optimization of the Rayleigh quotient [cf.Eq. (E2)] yields the SEE in the first form as La 1 ,...,a j−1 ,a j+1 ,...,a N |a j = g|a j (E3) for j = 1, . . ., N. The SPI does not evaluate this first form; rather, it solves Eq. ( 6), the second form of the SEE, which has been shown to be equivalent to Eq. (E3) (a comprehensive proof can be found in the Supplement Material to Ref. [40]).
FIG.2.Flowchart of the SPI algorithm.Branches in the algorithm, either "while" loops or "if" conditions, are represented by magenta diamonds.Yellow rectangles represent assignments and function calls.Entry and exit points of the algorithm are shown as cyan ellipses.Box 1 refers to the input, which includes an operator and the number N of parties.A vector is generated in 2 to serve as our starting vector.If the convergence criterion3 , studied in Sec.III B 4, is not met, we generate another N-partite vector in4 .For N = 1 in 5 , the algorithm corresponds to the power iteration (PI) which finds the standard eigenvector to the maximal eigenvalue (dashed box).For N > 1, our extension consists of three essential parts (gray areas), which are the forward iteration6 and the backward iteration8 as well as a recursion calling the SPI for N − 1 parties in7 .Eventually, when 3 is satisfied, the desired separability eigenvector is the output of our SPI algorithm and returned in12 .

FIG. 3 .
FIG. 3. Benchmark results comparing the SPI (magenta line) to a brute-force approach, which is a genetic algorithm (cyan line).Top panel: The results for the runtime comparison between both approaches are shown when scaling the dimensions of each subsystem for a bipartite state, N = 2 and d 1 = d 2 = d.Bottom panel: The corresponding results are shown when scaling the number N of parties, which are qubit systems (d 1 = • • • = d N = 2).It can be seen that the SPI performs better than the competing approach by several orders of magnitude.

FIG. 4 .
FIG.4.Results of the entanglement test for the bound-entangled, two-qutrit Horodecki state.No state ρα with 2 ≤ α ≤ 3 (blank area) has been detected as entangled.In the magenta areas, the criterion Eq. (33) certifies entanglement for the given combination of α and β , which does not hold true for the cyan areas.
Now, we can perform the next cycle, which results in |γ| 2 > |γ | 2 > |γ | 2 .In fact, performing s steps of the SPI, we get vectors |a