{'fraction_non_alphanumeric': 0.028134228187919462, 'fraction_numerical': 0.015212527964205816, 'mean_word_length': 4.059856922937607, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 54, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extremal states of positive partial transpose in a system of three qubits

We have studied mixed states in the system of three qubits with the property that all their partial transposes are positive, these are called PPT states. We classify a PPT state by the ranks of the state itself and its three single partial transposes. In random numerical searches we find entangled PPT states with a large variety of rank combinations. For ranks equal to five or higher we find both extremal and nonextremal PPT states of nearly every rank combination, with the restriction that the square sum of the four ranks of an extremal PPT state can be at most 193. We have studied especially the rank four entangled PPT states, which are found to have rank four for every partial transpose. These states are all extremal, because of the previously known result that every PPT state of rank three or less is separable. We find two distinct classes of rank 4444 entangled PPT states, identified by a real valued quadratic expression invariant under local SL(2,C) transformations, mathematically equivalent to Lorentz transformations. This quadratic Lorentz invariant is nonzero for one class of states (type I in our terminology) and zero for the other class (type II). The previously known states based on unextendible product bases is a nongeneric subclass of the type I states. We present analytical constructions of states of both types, general enough to reproduce all the rank 4444 PPT states we have found numerically. We can not exclude the possibility that there exist nongeneric rank four PPT states that we do not find in our random numerical searches.


Introduction
As one of the most fascinating features of quantum mechanics entanglement has been intensively studied in the last decades [1]. The Bell inequalities [2] turned the philosophical discussion between Einstein and Bohr into a subject for experimental investigations [3,4,5,6]. The Bell inequalities are deduced from the hypothesis of local realism, and apply to statistical correlations in composite systems with two subsystems. One basic weakness of all such experiments is the inherent statistical uncertainties of the observed correlations.
The striking new feature of the experiments proposed by Greenberger, Horne, and Zeilinger [7] and Mermin [8], in composite systems with three or more subsystems, is that the correlations to be tested are absolute and no longer statistical. One single observation is sufficient to demonstrate that quantum mechanics violates local realism. Experiments of this kind have also been made [9].
A central problem in the study of entanglement is how to determine, theoretically or experimentally, whether a state in a composite system made up of two or more components is entangled or separable. The answer is simple for a pure state: it is entangled if it is not a tensor product. The separability problem for mixed states, on the other hand, has proven to be highly nontrivial [10] and does not yet have a solution which is satisfactory for practical use. The simplest and best known condition for the separability of a mixed state is the Peres criterion [11], which states that a separable state remains positive semidefinite under partial transpositions. This necessary condition is in general not sufficient, but applying it requires very little computational effort and the separability problem is therefore in essence reduced to determining whether mixed states with positive partial transposes (PPT states) are entangled or separable.
A partial transposition is a mathematical operation that transposes one or more factors in a tensor product, for example, The point is that it is a well defined operation even for a matrix which is not a tensor product (see Appendix A). A Hermitian matrix representing a mixed state of an n-partite system may appear in 2 n different versions related by partial transpositions, but these can be separated into two sets where every matrix in one set is the total transpose of a matrix from the other set. Since the total transpose T preserves eigenvalues, half of the 2 n partial transposes are superfluous when we evaluate whether or not a mixed state ρ is a PPT state.
In the present article we consider PPT states in a system of three qubits. In this system it is sufficient to watch the eigenvalues of ρ, ρ T 1 , ρ T 2 , and ρ T 3 , where T i is the partial transposition with respect to subsystem i. It is necessary to require that all these four are positive, since it is quite possible for three of them to be positive, while the fourth one is not. The four remaining partial transposes are obtained from these by a total transposition. For example, T 1 T 2 = T 3 T , because T = T 1 T 2 T 3 , partial transpositions commute, and T 2 i is the identity operation. We write D for the set of all unnormalized mixed states, and D 1 for the set of all mixed states normalized to unit trace. Similarly, we write S or S 1 for the set of all separable states, and P or P 1 for the set of all PPT states, unnormalized or normalized. All of these are convex sets.
An extremal point in a convex set is a point which is not a convex combination of other points in the set. The sets of normalized states, D 1 , S 1 , and P 1 are compact and hence completely described by their extremal points, in the sense that all nonextremal points may be written as convex combinations of the extremal points. The extremal points of D 1 are the pure states, and the extremal points of S 1 are the pure product states. The inclusions together with the fact that the pure product states are extremal points in both S 1 and D 1 , imply that the pure product states are also extremal points of P 1 . It is easy to prove that all pure nonproduct states are not in P 1 and therefore entangled. The pure product states are the only PPT states of rank one, since only pure states have rank one. Because S 1 and P 1 are not identical in the three qubit system, P 1 must have extremal points of rank higher than one giving rise to all the entangled PPT states. These extremal entangled PPT states are almost completely unknown, and that situation motivates our study presented here. Previous studies of multipartite entanglement have been mostly concerned with pure states. Dür, Vidal, and Cirac [12] classified pure states in the three qubit system into six equivalence classes, based on the type of entanglement possessed by a state. A pure state is of the form ψψ † with ψ ∈ 8 . Two (unnormalized) vectors ψ and φ are considered equivalent if with V i ∈ SL(2, ). One class contains the separable (unentangled) pure states where ψ is a 2 × 2 × 2 dimensional product vector. Three other classes contain the biseparable pure states where ψ is a product vector in one of the three splittings into one system of one qubit and one system of two qubits. All states in these three classes have only bipartite entanglement. The last two classes contain states with two inequivalent types of genuine tripartite entanglement. There is the W class, exemplified by the unnormalized state, in Dirac notation, Finally there is the GHZ class, exemplified by the Greenberger-Horne-Zeilinger state In matrix notation, with the W and GHZ states are The corresponding classification scheme for mixed states in the three qubit system was introduced by Acín et al. [13]. The scheme involves four convex sets with an onion structure where each set is defined by including increasingly larger sets of pure states. The innermost set is the set S of all separable states, consisting of all states that are convex combinations of pure product states. The second set B includes all the biseparable pure states. The third set W includes also pure states with W entanglement. The fourth set includes the last remaining class of pure states, the GHZ entangled states. Because it includes all pure states the fourth set is D, the set of all mixed states. The authors conjecture that all entangled PPT states are members of the third set W.
Karnas and Lewenstein [14] proved that all states of ranks two and three in the three qubit system are separable. They applied the range criterion of entanglement for the three qubit system as a special case of the 2 × 2 × N dimensional system. Bennett et al. [15,16] introduced the PPT states based on unextendible product bases (UPBs) as examples of entangled PPT states in the three qubit system. These three qubit UPB states have rank four. They are entangled because there is no product vector in the range of a UPB state, and they are extremal PPT states because no entangled states of lower rank exist.

Sorting states into equivalence classes
We want to classify the three qubit states into SL ⊗ SL ⊗ SL equivalence classes, defining two unnormalized density matrices ρ and σ to be equivalent if with V i ∈ SL (2, ). This definition is useful because equivalent density matrices have the same entanglement properties, although quantitative measures of entanglement will be different. Qualitative properties will be the same, such as tensor product structure of pure states, decomposition of mixed states as convex combinations of pure states, rank and positivity of states and all their partial transposes, and so on. The relation between the group SL(2, ) and the group of continuous Lorentz transformations is well known, and is reviewed here in Appendix B. From a density matrix in the three qubit system we define one quadratic and four quartic real Lorentz invariants, so called because they are invariant under SL⊗SL⊗SL transformations as in equation (8). They are also invariant under partial transpositions, because a partial transposition may be interpreted as a parity transformation, which is a discrete Lorentz transformation. This means, for example, that the ratio between one quartic Lorentz invariant and the square of the quadratic invariant will have the same value for all the states in one equivalence class and all their partial transposes. Taking the ratio between Lorentz invariants is necessary in order to cancel out any normalization factor in the density matrix. If two density matrices are not in the same equivalence class, their inequivalence will most likely be revealed when we calculate their invariants.

Outline of this article
We have investigated both extremal and nonextremal PPT states in the system of three qubits using both numerical and analytical methods. We classify the states according to the ranks m 0 , m 1 , m 2 , m 3 of ρ, ρ T 1 , ρ T 2 , and ρ T 3 .
In Section 2 we present the two main numerical algorithms that we have used in random searches for PPT states. We have searched systematically for extremal states of unspecified ranks, and also for states of specified ranks that are not necessarily extremal.
In Section 3 we present results from the searches. We find PPT states with a wide variety of rank combinations. The rank 1111 states are the pure states. The states of ranks 2222 and 3333 are all separable, in agreement with previously known results. An interesting class of extremal PPT states are those of rank 4444, they contain genuine tripartite entanglement, since they are separable in any bipartite splitting of the three qubit system. We have studied these states in detail and shown how to construct them analytically, as reported in Section 4.
All the other states we find have all ranks equal to five or higher. We find PPT states, extremal or nonextremal, of very nearly every rank combination. For an extremal PPT state the square sum of ranks can be at most 193. Up to this upper limit we also find extremal PPT states of very nearly every rank combination.
In Section 4 we present our understanding of the rank 4444 extremal PPT states. It turns out that the key to understanding them is the fact that they are biseparable in three different ways [17]. There are two distinct classes of such states, we call them simply type I and type II. The most obvious distinction is that the quadratic Lorentz invariant is nonzero for states of type I and zero for states of type II.
When we sort these states further into SL ⊗ SL ⊗ SL equivalence classes, we find that every equivalence class of type I states contains a density matrix which is real and symmetric under all partial transposistions. We define this matrix to be a standard form for all the states in the equivalence class. The UPB states are a subclass of the type I states, but they are not sufficiently generic that we find any of them in our random searches.
It is straightforward to construct the most general density matrix which is biseparable, real and symmetric under all partial transpositions. The empirical result of our numerical searches is that this construction reproduces all the rank 4444 PPT states of type I. The SL ⊗ SL ⊗ SL equivalence classes are parametrized by seven continuous real parameters.
A rank 4444 PPT state of type II can not in general be transformed to a real form. We have studied our numerical examples of such states and observed several special properties that they have. Based on these observations we present an analytical construction general enough to reproduce all the type II states found numerically. The SL ⊗ SL ⊗ SL equivalence classes in this case are parametrized by one continuous complex parameter.

Numerical methods
Most of the PPT states that we have studied numerically were found by means of two algorithms that search for extremal PPT states or for PPT states with a specified combination of ranks. In this section we present briefly these two main algorithms.

Random search for extremal PPT states
Extremal states in P 1 , the convex set of PPT states of unit trace, are found by an iterative process where we start with any PPT state, pick a random search direction restricted to the unique face where this state is an interior point, and follow this direction to the edge of the face, which is a face of lower dimension. The next search direction is restricted to this new face. The procedure is repeated until a face of dimension zero, which is an extremal point of P 1 , is located. An extremal point of P 1 defines an extremal ray, a one dimensional face, of P, the cone of unnormalized PPT states.
A valid search direction from a given PPT state ρ ∈ P 1 is a nonzero matrix σ such that ρ+ǫσ ∈ P 1 for every ǫ in some finite interval, ǫ 1 ≤ ǫ ≤ ǫ 2 , with ǫ 1 < 0 < ǫ 2 . This means that σ is a traceless Hermitian matrix satisfying the four equations where P i is the orthogonal projection onto Img ρ T i , and where we write ρ T 0 = ρ. The real vector space H of Hermitian N × N matrices has dimension N 2 and is a Hilbert space with the scalar product A, B = B, A = Tr(AB) .
The constraints in equation (9) can now be written as eigenvalue equations P i σ = σ for linear projection operators P i on H, symmetric with respect to this scalar product, defined as follows, The four eigenvalue equations may be combined into one single eigenvalue equation This equation always has the trivial solution σ = ρ, but ρ is not traceless. If ρ is the only solution it is an extremal point of P 1 . The number of linearly independent traceless solutions to this eigenvalue equation is the dimension of the restricting face of P 1 . We find a random search direction by finding a complete set of eigenvectors solving equation (12), then choosing σ 0 as a random linear combination of these, and finally obtaining a traceless solution as The bipartite version of this algorithm is described in more detail in [19].

Separability test for low rank states
A modified version of this algorithm can be used as a separability test for low rank states. If a PPT state is not extremal it can be written as a convex combination of extremal PPT states of lower rank. A set of such states can be found by searching in both directions σ and −σ every time a search direction σ is chosen. The original state is separable if all the extremal states found in this process are pure. In general it may be possible to write a separable state as a convex combination involving some mixed extremal states. However, for low rank states few mixed extremal states of lower rank are available and these "false negatives" are therfore less likely to occur.

Upper limit on the ranks
We derive an upper limit on the ranks of extremal PPT states by counting the number of constraints on σ. Since each operator P i is a projection with eigenvalues 0 and 1, the constraint equation P i σ = σ places a number of constraints on σ equal to the dimension of the kernel of P i . If ρ T i has rank m i , then the rank of P i is m 2 i , and the dimension of its kernel is N 2 − m 2 i . Thus, for a given ρ the total number of constraints on σ, apart from the zero trace condition, is These constraints need not be linearly independent, hence this is an upper limit on the number of linearly independent constraints. Since ρ is an extremal PPT state if and only if it is the only Hermitian matrix satisfying the constraints, the number of linearly independent constraints when ρ is extremal must be N 2 − 1. The total number of constraints must be at least as large. This implies the following limit on the ranks of extremal PPT states and their partial transposes in a three qubit system, The derivation of the formula involves a sum over the independent partial transposes, the number of which is 2 n−1 for an n-partite system. The formula can therefore be generalized to for an n-partite system of total dimension N . For large n this inequality is a very mild restriction on the ranks.

Search for PPT states of specified low ranks
The second algorithm we describe take as input the desired ranks m 0 , m 1 , m 2 , m 3 of ρ and its partial transposes. After we find a ρ with the specified ranks we may check, by the algorithm just described, whether or not it is extremal. The equations defining a Hermitian matrix ρ as a PPT state of the given ranks are of the form µ = 0, where µ is a list of all the lowest eigenvalues of ρ and its partial transposes. We include the N − m i lowest eigenvalues of ρ T i (recall our definition ρ T 0 = ρ). When we write ρ as a linear combination of matrices M j forming a basis for the real vector space of Hermitian matrices, the equations are of the form µ i (x) = 0. There are N 2 real variables x j , and the number of equations is

A minimization problem
One way to solve the equations µ i (x) = 0 is to minimize the function and obtain a minimum value of zero. We may use a random search algorithm for finding the minimum. This approach may be good enough for many purposes, although not very efficient, and it has the advantage of being easy to program. Note that a minimum point of a function is in general not very precisely determined numerically, because the function varies quadratically around its minimum. This problem does not arise here because the desired minimal value is zero.

The conjugate gradient method
A more refined and efficient iterative approach is to linearize the equations and solve in each iteration the approximate, linearized equations using the conjugate gradient method. A similar algorithm for the bipartite case is described in [20].
Given an approximate solution x we try to find a better solution x + ∆x using the linearized equation where the matrix elements of B(x) are B ij (x) = ∂µ i (x)/∂x j . We multiply this equation by B T , and get another equation where A = B T B and b = −B T µ. The matrix A is real and symmetric as well as positive semidefinite. It is likely to be singular, but the last equation may anyway be solved by the conjugate gradient method. We calculate the components of B using first order perturbation theory. If for example µ i is the eigenvalue λ k of ρ then we use the formula Here ψ k is the eigenvector corresponding to the eigenvalue λ k of ρ. This formula is based on nondegenerate perturbation theory, which is strictly speaking not valid in the present case where we try to make many eigenvalues simultaneously equal to zero. The formula nevertheless works well in practice. We use similar formulas for the eigenvalues of the partial transposes of ρ.

Ranks of PPT states found numerically
Numerical searches for PPT states of specified ranks were successful for most of the rank combinations. Note however that our search algorithm might occasionally produce PPT states of lower ranks than the ones specified. Note also that for the purpose of classification of states there is full symmetry between ρ and all its seven partial transposes (including the double partial transposes and the total transpose). Furthermore, with three identical subsystems there is full symmetry between states related by an interchange of the subsystems. Therefore when we list the ranks m 0 , m 1 , m 2 , m 3 we use the convention that m 0 ≤ m 1 ≤ m 2 ≤ m 3 . The PPT states found with specified ranks might be separable or entangled, extremal or nonextremal. We used them as starting points in searches for extremal PPT states. This would either show them to be extremal, or return extremal PPT states of lower ranks.
A list of all confirmed rank combinations and the corresponding square sums of ranks is given in Table 1. By equation (15), for an extremal PPT state the square sum of ranks can not be larger than 193. We found no PPT states at all with the very asymmetric rank combinations 5568, 5588, or 5888. Otherwise, the rank combinations 2222, 3333, and 5688 are the only cases where the inequality allows extremal PPT states to exist, but we did not find any.  Table 1: A list of all confirmed rank combinations and the corresponding sums of squared ranks. The upper limit on the sum of squared ranks of extremal PPT states, as given by equation (15), is 193. There are three cases, each marked by an asterisk, where this limit allows extremal states to exist but only nonextremal states were found.
A similar study on the ranks of extremal states in low dimensional bipartite systems is presented in [20]. The structure of the rank combinations of our three qubit states closely resembles the structure in the 3×3, 3×4, 3×5, and 4×4 systems. For the lowest rank combinations only states where all partial transposes have the same rank are found, and below some threshold rank, four in our case, all these states are separable. Above the threshold extremal states are found for essentially all rank combinations where PPT states exist and where the upper limit given in equation (15) holds. The sole exception is the rank combination 5688 where extremal states are allowed by the limit, but we were unable to confirm their existence. We consider it likely that extremal states of this rank combination do exist, but are difficult to find numerically by our methods.

Rank 4444 PPT states
The rank four extremal PPT states are of special interest because they are the entangled PPT states of lowest rank. Thus, any rank four PPT state is either extremal or separable. Another special property of the rank four PPT states is that all their partial transposes have the same rank, we find no PPT states of rank 4445, for example.
Nine out of 196 extremal rank four states found in a random numerical search stand out because their quadratic Lorentz invariant is zero, which means in practice of order 10 −14 to 10 −16 . We find a qualitative difference between the states with nonzero quadratic invariant and those with the invariant equal to zero. In particular, in our numerical sample we see no continuous transition from invariants of order one to those of order 10 −14 .
We will discuss the rank 4444 entangled PPT states in much more detail in Section 4.

PPT states with special symmetries
An obvious way to construct a PPT state is to construct a positive Hermitian matrix which is symmetric under all three partial transpositions T 1 , T 2 , T 3 . These symmetries together imply symmetry under the total transposition T = T 1 T 2 T 3 , which makes the matrix real. The most general Hermitian matrix symmetric under all three partial transpositions contains 27 real parameters. See Appendix A.
Completely symmetric PPT states of rank eight are easily generated as follows. Generate a random complex matrix, and multiply it with its Hermitian conjugate to obtain a Hermitian and positive matrix. Add to this matrix all its seven partial transposes. The resulting matrix is completely symmetric and has a high probability of being positive.
Another way to generate a completely symmetric PPT state of any rank m is the following. Start with any matrix in the 27 dimensional space of completely symmetric matrices. Then take the square sum of its 8 − m lowest eigenvalues, and minimize this by varying the matrix so as to obtain the minimum value of zero. We used this method to create states of rank four, five, six, and seven. The rank four states were all extremal, while all the higher rank states were nonextremal.
It is also possible to relax the symmetry requirements and only require symmetry under two partial transpositions, say T 1 and T 2 . This raises the number of free parameters in a Hermitian matrix from 27 to 36 and allows it to be complex. A state ρ with the symmetries ρ = ρ T 1 = ρ T 2 is still guaranteed to be a PPT state, because ρ T 3 = ρ T has the same eigenvalues as ρ. Only about 75% of the rank four states generated with these symmetries were extremal.

Decomposition
The separability test described in Subsection 2.1.1 is of particular interest for the nonextremal rank 5555 states, which are the states of lowest rank where the test is nontrivial. Among the states of this kind found in random searches for states of specified ranks we found both separable and entangled states. Every entangled state was found to be an interior point of a one dimensional face of P 1 , this face is then a line segment with a pure state at one end and a rank 4444 extremal PPT state at the other end. Every separable state was found to be an interior point of a four dimensional face of P 1 , as expected when it is a convex combination of five pure states.
We also tested the rank five states that were constructed so as to be symmetric under all partial transpositions. The result is that every such state is an interior point of a face of dimension five, and decomposes always into two rank four extremal states having the same symmetry.
We found the last result surprising, and therefore applied the same test to symmetric states with higher ranks. The rank six states behave in a very similar manner, decomposing into four symmetric rank four states that are always extremal. The rank seven states, on the other hand, decompose mainly into rank 6677 states which obviously do not possess the same symmetry. Every rank eight state decomposes into 32 rank 6777 states that also do not possess any symmetries.

Rank four entangled PPT states as biseparable states
We find empirically that the entangled and extremal rank 4444 PPT states in dimension 2 × 2 × 2 fall into two qualitatively different classes, which we choose to call types I and II. By definition, a state ρ of type I has a nonzero (always positive) value of the quadratic Lorentz invariant whereas this invariant vanishes for a state ρ of type II. See Appendix B. As we shall see, another characteristic difference is that a state of type I can be transformed by a product transformation to a standard form where it is real and symmetric under all partial transpositions, whereas a state of type II can not in general be transformed to a real form. Type I includes as a special case the so called UPB states constructed from unextendible product bases (UPBs), introduced in [15,16], which we will also describe here. The basic fact leading to an understanding of these states is that if ρ is a density matrix of rank four in dimension 2 × 2 × 2, and if ρ T 1 ≥ 0, then ρ is separable in dimension 2 × 4. Similar results hold of course for the cases ρ T 2 ≥ 0 and ρ T 3 ≥ 0. This has been proved analytically [17] and also verified numerically [20]. We will describe in this section how to use the threefold biseparability in order to understand our rank 4444 PPT states.
We first discuss some common properties of the two types of states, and a new method for constructing them numerically, before we turn to the more detailed mathematical understanding of the two types separately, including the UPB states.

General considerations
Consider a generic four dimensional subspace U ⊂ 8 . It contains exactly four product vectors of dimension 2 × 4, which form a basis for U , and which may be calculated by the method described in Appendix E. It contains also exactly four product vectors where y j ∈ 2 , v j ∈ 4 , and where we use the split tensor product defined in Appendix F. These product vectors form another basis for U . And it contains exactly four product vectors of dimension 4 × 2, forming a third basis for U . A density matrix ρ of rank four with Img ρ = U and ρ T 1 ≥ 0 must be biseparable in dimension 2 × 4, and hence of the form Note that although ρ is separable as a bipartite state, it is necessarily entangled as a tripartite state, except in the special case when the four vectors e i are 2 × 2 × 2 product vectors, which means that they are identical to the vectors f j , and also identical to the vectors g k . If ρ T 2 ≥ 0, then ρ must be biseparable by the split tensor product, and we must have If ρ T 3 ≥ 0, then ρ must be biseparable in dimension 4 × 2, and

Constructing PPT states numerically as biseparable states
We find that with a generic subspace U any two of the three equations (27), (28), and (29) are incompatible, they have no common nonzero solution for ρ. If we require that two of the three equations, for example (27) and (28), should be compatible, then this restricts the subspace U . It is easy to find subspaces where two of the three equations are compatible, but then the third equation will in general be incompatible with the two. This means that we may have, for example, ρ ≥ 0, ρ T 1 ≥ 0, and ρ T 2 ≥ 0, but ρ T 3 ≥ 0.
When we want a PPT state with ρ ≥ 0 and ρ T i ≥ 0 for i = 1, 2, 3, we need to find a subspace U where the equations (27) and (28) are compatible, and at the same time the equations (27) and (29). We will describe next how to solve these compatibility problems numerically.
The equations (27) and (28), apart from the positivity conditions on the coefficients λ i and µ j , are compatible when there is at least one linear dependence between the eight Hermitian 8 × 8 matrices e i e † i and f j f † j . We do the calculation numerically by converting the eight matrices to eight vectors in Ê 64 and making a singular value decomposition of the 64 × 8 matrix. The eight singular values are non-negative, by definition, and we want the smallest one to be zero. This is a minimization problem where we vary the subspace U in order to minimize the smallest singular value. We solve this in practice by a random search method. An output of the singular value decomposition, after minimization, is the two sets of coefficients λ i and µ j corresponding to the singular value zero. However, the singular value decomposition puts no restriction on the signs of the coefficients λ i and µ j . If we get a solution with all λ i < 0, then we simply switch all four signs in order to get all λ i > 0. If, on the other hand, we get four coefficients λ i with both signs, then we have a matrix ρ which is biseparable in two ways, but we can not make it positive semidefinite. This is then a subspace U which we can not use.
The procedure for making the equations (27) and (29) compatible is the same, we formulate it as a second minimization problem for a smallest singular value. In order to obtain a subspace U solving both problems simultaneously we simply minimize the sum of the smallest singular value from the first problem and the smallest singular value from the second problem.
The numerical procedure described here for constructing PPT states of rank four works well, although the random search for a minimum may take some time. In practice, we always get one singular value zero with a corresponding unique solution for the coefficients λ i , µ j , and ν k . However, it seems that each coefficient λ i for i = 1, 2, 3, 4 comes with an essentially random sign. We need the same sign for all four coefficients in order to satisfy the positivity condition ρ ≥ 0, and on the average we have to try about eight times before we succeed once.
The majority of the rank four PPT states found numerically by this method are of type I, having a nonzero value of the quadratic invariant. We find however also a significant fraction of states of type II, with vanishing quadratic invariant. As already remarked, these are two distinct classes of rank four PPT states. We will proceed next to discuss them separately, after discussing the so called UPB states.

Standard form
In the generic case there always exists a product transformation uniquely defined up to normalization, transforming the two dimensional vectors x, y, z into V 1 x, V 2 y, V 3 z such that the transformed vectors after suitable normalizations have the standard form defined in Appendix C, equation (100), with complex parameters t 1 , t 2 , t 3 , If ρ is a density matrix of rank four with Img ρ = U , then the product transformation V transforms ρ into V ρV † , which in our terminology is SL ⊗ SL ⊗ SL equivalent to ρ. This transformed matrix will then be a standard form for ρ. Note that the standard form is not completely unique, because it depends on the ordering within each of the three sets of product vectors e i , f j , and g k . Thus there is a discrete ambiguity.

States constructed from unextendible product bases
The construction of rank four entangled PPT states in dimension 2 × 2 × 2 from unextendible product bases (UPBs) is due to Bennett et al. [15,16]. It has been shown, both numerically [21,23] and analytically [24,25], that the UPB construction in dimension 3×3 produces all the rank four entangled PPT states there. In contrast, we find that the UPB states in dimension 2 × 2 × 2 are only a subclass of the rank four PPT states of type I, with nonvanishing quadratic invariant.
An unextendible product basis (a UPB) is a maximal set of orthogonal product vectors which is not a complete basis for the vector space. In other words, it is an orthogonal basis of product vectors for a subspace, with the property that the orthogonal subspace contains no product vectors. As shown in Appendix D a UPB in dimension 2 × 2 × 2 consists of four product vectors From these vectors, normalized to unit length, we construct a PPT state ρ as a normalized projection operator onto the orthogonal subspace, where ½ is the 8 × 8 unit matrix. This is obviously a PPT state, since ρ T 1 , for example, is again an orthogonal projection, with The state ρ constructed in this way is obviously entangled, since Img ρ contains no product vector, and a characteristic property of a separable state σ is that Img σ is spanned by product vectors. From ρ we get other rank four entangled PPT states of the form ρ = aV ρV † where V is a product matrix, V = V 1 ⊗ V 2 ⊗ V 3 , and a is a normalization constant. We have then that The kernel of ρ is spanned by the four orthogonal product vectors ψ i . The kernel of ρ is spanned by the product vectors which are orthogonal if and only if V is unitary. The existence of a basis of product vectors in the kernel, orthogonal or not, is the characteristic property of the rank four PPT states of the UPB type and their SL ⊗ SL ⊗ SL transforms. These states are not generic, because a generic subspace of dimension four contains no product vector. In fact, a generic and suitably normalized product vector contains three complex parameters a, b, c. In order to lie in a given four dimensional subspace it has to satisfy 8 − 4 = 4 linear constraints, which it can not do in general with only three parameters. Using numerical methods, we have searched for product vectors in the kernels of our extremal rank four PPT states, but found none. We conclude that there are no states equivalent to UPB states in our numerical sample. The UPB states exist, but are not sufficiently generic to be easily found in random numerical searches.
The UPB has a standard form where all the vectors x, y, z are real and given by three continuous real parameters (three angles), as in equation (109). Therefore ρ has a standard form in which it is symmetric under all partial transpositions, ρ = ρ T 1 = ρ T 2 = ρ T 3 .

Rank four PPT states with nonzero quadratic invariant
We turn now from the UPB states to the more general class of rank 4444 extremal PPT states with nonvanishing quadratic invariant. The key to understanding these states comes from a study of our numerical examples. It turns out that all the states we have found numerically may be transformed to a standard form where they are symmetric under all partial transformations, and hence real.
Let ρ be a state with nonvanishing quadratic invariant. We compute the three sets of four product vectors e i = x i ⊗ u i , f j = y j ⊗ s v j , and g k = w k ⊗ z k in Img ρ, with x i , y j , z k ∈ 2 and u i , v j , w k ∈ 4 , as already defined in Subsection 4.1. Then we use a product transformation V = V 1 ⊗ V 1 ⊗ V 3 to transform the two dimensional vectors x, y, z to the standard form given in equation (31).
We find empirically that we always get real values for the parameters t 1 , t 2 , t 3 , and that the same transformation V always makes u, v, w real. Moreover, the transformed matrix V ρV † always becomes symmetric under all partial transpositions. These are the key observations. We have no mathematical proof that ρ must have these properties if its quadratic invariant is nonzero, but we take it as a working hypothesis. As such it is very powerful, and we will see next how to construct the most general PPT state with these properties.

Explicit construction
In equation (27) we may absorb a factor √ λ i into the vector e i = x i ⊗ u i and write The problem we face is to construct a density matrix ρ which has this form with x and u real, and is symmetric under all partial transpositions, ρ = ρ T 1 = ρ T 2 = ρ T 3 . We also want x to have the standard form given in equation (31), so that Then ρ has the following form, Let us see what happens if we simply take u to be a random 4 × 4 real matrix. As shown in Appendix A, ρ as given in equation (41) is symmetric under all partial transpositions if and only if the 4 × 4 matrix A is real and has the following general form, a 1 a 5 a 6 a 7 a 5 a 2 a 7 a 8 a 6 a 7 a 3 a 9 a 7 a 8 a 9 a 4 and B and C are also real and have similar forms. In the case of A the only condition which is not automatically satisfied as a consequence of equation (42) Since the condition that ensures the right form for B, equation (44), may be rewritten as We solve it simply by choosing Next, the condition that ensures the right form for A is the following, Define now If α 2 > 0, then we replace u 1 by αu 1 and have a solution of equation (48). If instead α 2 < 0, then we change the sign of u T 1 (ǫ ⊗ ǫ) u 1 , for example by changing the signs of the first two components of u 1 . With this new vector u 1 we have α 2 > 0, and again we solve equation (48) by substituting αu 1 for u 1 .
The condition that ensures the right form for C is the following, We solve it for u 2 in the same way as we solved equation (48) for u 1 . This completes the construction of ρ, except that we should normalize ρ to have unit trace.

Parameter counting
A random 4×4 real matrix u contains 16 freely variable real parameters. In the process of constructing ρ to be symmetric under partial transformations we had to redefine the normalizations of u 1 and u 2 , this reduces the number of free parameters to 14. The final normalization of ρ reduces the number of parameters to 13. Note that when we take x to have its standard form and choose u randomly, the resulting vectors y and z will most likely not have their standard forms. In order to transform them to standard forms we will need two SL(2, Ê) transformations, containing altogether six real parameters. After we fix the standard forms of x, y, z there are no continuous degrees of freedom left in the product transformation V . There is however some discrete freedom left, because we may permute the product vectors e i , as well as the vectors f j and the vectors g k .
We conclude that there is a family of SL ⊗ SL ⊗ SL equivalence classes of rank 4444 extremal PPT states with nonzero quadratic invariant described by 13 − 6 = 7 continuous real parameters.

States with vanishing quadratic invariant
In our numerical random searches for PPT states of rank four we found a special class of rank 4444 PPT states with the property that the quadratic Lorentz invariant of each state ρ vanishes, that is, See Appendix B. In this subsection we present our understanding of these states. We find that such a state has a standard form described by one continuous complex parameter. The spectral representation of ρ, with four eigenvalues λ i > 0 and corresponding orthonormal eigenvectors η i , gives that Since the four eigenvalues are positive, the only way to have Tr(ρ T EρE) = 0 is to have all the eigenvectors orthogonal in the scalar product defined by E, Since E T = −E, this scalar product is antisymmetric, ψ T Eφ = −φ T Eψ, hence every vector is orthogonal to itself, ψ T Eψ = 0. Equation (54) means that the antisymmetric scalar product vanishes identically in the four dimensional subspace Img ρ spanned by the eigenvectors η i with i = 1, 2, 3, 4. We may use this observation in order to search numerically for PPT states of this type.

A random search method
It is a straightforward exercise to construct a random four dimensional subspace U ⊂ 8 where the antisymmetric scalar product vanishes identically. We may start with a random normalized vector ψ 1 . The two conditions ψ † 1 ψ = 0, ψ T 1 Eψ = 0 restrict ψ to a six dimensional subspace, and we choose ψ 2 as a random normalized vector in this subspace. Next, the conditions ψ † i ψ = 0, ψ T i Eψ = 0 for i = 1, 2 restrict ψ to a four dimensional subspace, and we choose ψ 3 as a random normalized vector in this subspace. Finally, we choose ψ 4 as a random normalized vector in a two dimensional subspace. In this way we get vectors ψ 1 , ψ 2 , ψ 3 , ψ 4 such that The vectors φ i = Eψ * i for i = 1, 2, 3, 4 lie in the orthogonal subspace U ⊥ , since and they satisfy the same relations as the vectors ψ i , For a given subspace U we may try to construct numerically a PPT state ρ with Img ρ = U , by the method described in Subsection 4.1.1. Again we find that the construction fails in general, but we may select a suitable subspace with vanishing antisymmetric scalar product where the construction succeeds.
As remarked in Subsection 4.1.1, when we search for a random subspace U where we can find a matrix ρ satisfying all three biseparability conditions (27), (28), and (29), disregarding the positivity conditions, then it happens only in about one case out of eight that the ρ we find is positive. Surprisingly, when we now do the same search restricted to subspaces with the special property that the antisymmetric scalar product vanishes identically, then the positivity conditions hold, not every time, but almost every time.

Explicit construction
By studying our numerical examples of such states we find empirically that they have certain properties which enable us to construct them explicitly. Again we have no strict proof that these properties are necessary, but it is a very powerful working hypothesis to assume that they hold.
Like before, we introduce three sets of basis vectors for the subspace Img ρ, The vanishing of the antisymmetric scalar product in Img ρ means that We have that x T i ǫx j = 0 for i = j, but x T i ǫx j = 0 for i = j. Hence we must have Note that ǫ is antisymmetric, but ǫ ⊗ ǫ is symmetric, In our numerical examples we see that equation (59) is solved in the following remarkable way. The vectors u i with i = 1, 2, 3, 4 can not all be product vectors, because then ρ would be separable. Instead, each u i is a linear combination of two product vectors, as follows, There is no sum here over the indices k, l, m, n, and we have always k = m and l = n. For each value of i, six different combinations occur for klmn, and we may order the vectors y and z in such a way that we get the index combinations listed in Table 2. Since we are guaranteed that u i and u j satisfy the orthogonality relation in equation (59) if We see from Table 2 that two vectors u i and u j are always orthogonal in two ways. For example, u 1 = a 112 y 1 ⊗ z 2 + a 121 y 2 ⊗ z 1 = a 134 y 3 ⊗ z 4 + a 143 y 4 ⊗ z 3 , u 2 = a 211 y 1 ⊗ z 1 + a 222 y 2 ⊗ z 2 = a 233 y 3 ⊗ z 3 + a 244 y 4 ⊗ z 4 .
We now transform x, y, z to the standard form defined in equation (31), with complex parameters t 1 , t 2 , t 3 . We find that the linear dependencies listed in Table 2 require that t 2 = t 3 , and when this  relation holds they give a unique solution for u depending on the single complex parameter t = t 2 = t 3 . The solution is The condition that the four vectors g k = w k ⊗ z k must be linear combinations of the vectors e i = x i ⊗ u i requires that also t 1 = t. We find the overall solution v = w = u, giving the vectors Thus we arrive at the following explicit standard form for ρ, depending on the single complex parameter t, with and a = 1 5|t| 4 + 10|t| 2 + 1 + (3|t| 2 + 1) |1 + t| 2 . (69)

The effect of partial transposition
The effect of the partial transposition T 1 , for example, on the PPT state ρ given in equation (67) is that with e i = x * i ⊗ u i . The two dimensional vectors x are complex conjugated while the four dimensional vectors u are unchanged. If x and u have their standard forms as given in equation (31) (with t 1 = t) and in equation (65), it means that the parameter t is complex conjugated in x but not in u.
We find that the complex conjugation of the standard form of u is equivalent to a linear transformation on product form. That is, there exists a 2 × 2 matrix W such that with four complex normalization constants C i having different phases but equal absolute values, |C i | = |C j | for i, j = 1, 2, 3, 4. There are four solutions for W , since the general solution contains two arbitrary signs ǫ 1 = ±1 and ǫ 2 = ±1.
It follows when we define V = ½ ⊗ W ⊗ W , with the 2 × 2 unit matrix ½, and introduce a This shows that the partial transpose ρ T 1 has a standard form like the standard form of ρ, but with the complex conjugated parameter value t * . The same observation applies of course also to the partial transposes ρ T 2 and ρ T 3 , and to the total transpose ρ T = ρ * .

Summary and outlook
We have presented here a numerical survey of entangled PPT states in the system of three qubits, with an emphasis on the extremal PPT states. Important characteristics of a state are the ranks of the state itself and of its three single partial transposes. For ranks equal to five or higher very few rank combinations are missing, the variety could not be much larger than what we see. This is true also for the extremal PPT states, where there is an upper bound of 193 for the square sum of ranks. The existence of a large variety of very different states indicates that it will not be easy to understand all the extremal PPT states analytically, in the same way as we understand the pure states.
A natural place to start a project to gain some analytical understanding is with the rank 4444 states, which are the lowest rank entangled PPT states. Another good reason for our interest in these states is that they have genuine tripartite entanglement, since they are separable with respect to any bipartition. An unanswered question about the rank four states is why the rank combination 4444 is the only one observed.
The most important advance reported here is that we have uncovered the mathematical structure behind the rank 4444 entangled PPT states, and shown how to reproduce analytically all the states found numerically. They fall into two distinct classes with very different analytical representations. A sobering fact is that the known states of the UPB type were not found numerically, although they appear as just a special case in our classification scheme. We take this as a reminder that there may exist unknown types of rank four PPT states that we have missed because they are not generic.
We intend to continue our project and try to understand for example the rank 5555 extremal PPT states in the three qubit system. The numerical study of higher dimensional tripartite systems quickly becomes impractical, simply because the dimensions grow rapidly.

A Partial transpositions
Write an 8 × 8 complex matrix as a 4 × 4 matrix of 2 × 2 matrices, The partial transposition T 1 moves 4 × 4 submatrices, T 2 moves 2 × 2 submatrices within 4 × 4 submatrices, T 3 transposes the 2 × 2 submatrices, If X is symmetric under all three partial transpositions, then it has the form with A T = A, B T = B, and so on. In particular, X is symmetric, X T = X T 1 T 2 T 3 = X. If in addition X is Hermitian, then it is real.

B SL(2, ), Lorentz transformations, and Lorentz invariants
Let ǫ be the two dimensional Levi-Civita symbol, It is a square root of −1, ǫ 2 = −½. The Lie group SL(2, ) consists of all 2 × 2 complex matrices with unit determinant. If V ∈ SL(2, ), then Thus V ǫV T = ǫ, in this sense ǫ is an invariant tensor under SL(2, ) transformations. A general 2 × 2 Hermitian matrix may be written as where x µ is a real fourvector, σ 0 = ½ is the unit matrix, and σ j for j = 1, 2, 3 are the Pauli matrices.
The determinant of X is where g µν is the metric tensor, More generally, for X = x µ σ µ and Y = y µ σ µ we have The transformation X → X = V XV † with det V = 1 is a continuous Lorentz transformation. It leaves the determinant invariant, and leaves the scalar product between two fourvectors invariant because det V T = det V † = 1, hence V T ǫV = V † ǫV * = ǫ and Tr( X T ǫ Y ǫ) = Tr(X T (V T ǫV )Y (V † ǫV * )) = Tr(X T ǫY ǫ) .
The parity inversion x 2 = −x 2 takes the form X = X T and leaves the scalar product invariant, although it is not of the form X = V XV † . In 8 = 2 ⊗ 2 ⊗ 2 the antisymmetric tensor is invariant under SL⊗SL⊗SL transformations, in the sense that V EV T = E when V = V 1 ⊗V 2 ⊗V 3 and V 1 , V 2 , V 3 ∈ SL(2, ).
A general 8 × 8 Hermitian matrix may be written as with 4 × 4 × 4 = 64 real coefficients A product transformation of the form A = V AV † with V = V 1 ⊗ V 2 ⊗ V 3 , as above, acts as three independent continuous Lorentz transformations on the three Lorentz indices. Note that the partial transpositions are discrete Lorentz transformations, since they are parity inversions. For Hermitian matrices A, B the quantity is real and invariant under the product transformations A = V AV † , B = V BV † . It is also invariant under all three partial transpositions. Note that each Lorentz index on a tensor a λµν represents its own subsystem and can therefore only be contracted against the corresponding index on the tensor b λµν . A density matrix ρ in dimension 2 × 2 × 2 has one quadratic Lorentz invariant For pure states ρ i = ψ i ψ † i we have that The inequality I 2 ≥ 0 follows because ρ is always a convex combination of pure states, and hence There are five different fourth order invariants obtained by different combinations of index contractions, but one of these is simply the square of the second order invariant. The four new invariants can be written as I 41 = ρ µνλ ρ µνγ ρ αβγ ρ αβλ , I 42 = ρ µνλ ρ µβλ ρ αβγ ρ ανγ , Note that all these Lorentz invariants are invariant under SL⊗SL⊗SL transformations of a density matrix without subsequent normalization to unit trace. Division of the fourth order invariants by the square of the second order invariant gives true invariants that are also independent of the normalization of the density matrix. They may be used in order to test whether two density matrices belong to the same SL ⊗ SL ⊗ SL equivalence class.

C Standard forms of sets of vectors
Given four vectors x i ∈ 2 , i = 1, 2, 3, 4. We write We consider here the generic case with det(x i , x j ) = 0 for i = j. Multiplication by the matrix and a subsequent normalization of the vectors gives the form with det(y i , y j ) = 0 for i = j, which means that t 1 t 2 = 0, 1. Multiplication by and normalization now gives the standard form with one variable parameter t = −t 1 t 2 = 0, −1. Since and this ratio of determinants is invariant under nonsingular linear transformations and vector normalizations, we have that We see that t will be complex in the generic case. In the special case where t is real and positive, we may multiply by and normalize so as to obtain the standard form where the vectors are real and pairwise orthogonal, w † i w j = w T i w j = 0 for i, j = 1, 2 and i, j = 3, 4. Instead of equation (100) we might have chosen one of two alternative standard forms, or z ′′ = 1 1 t ′′ 0 0 −1 1 1 .

D Unextendible product bases
An unextendible product basis (a UPB) in a subspace U of 2 ⊗ 2 ⊗ 2 is a set of orthogonal product vectors e i = x i ⊗ y i ⊗ z i spanning U such that no product vector is orthogonal to them all. We assume that the vectors x i , x j ∈ 2 are linearly independent when i = j, and similarly with y i , y j and z i , z j . Obviously, in order to have e i ⊥ e j we must have either x i ⊥ x j , y i ⊥ y j , or z i ⊥ z j . These conditions have a solution with four product vectors, unique up to permutations, as illustrated in Figure 1. In a similar way as discussed in Appendix C the vectors may always be transformed to the real standard forms where c i = cos θ i , s i = sin θ i , and θ 1 , θ 2 , θ 3 are three angular parameters.

E A problem of finding product vectors in a subspace
Given four vectors ψ j ∈ 8 with components ψ ij , i = 1, . . . , 8, j = 1, . . . , 4. We write a linear combination of them as a matrix product, with α ∈ 4 . Assume that φ is a tensor product The presence of the first factor in the tensor product implies the equation where µ = a/b, and A and B are the following 4 × 4 matrices, The presence of the second factor implies the equation where µ = c/d, and where The presence of the third factor implies the equation Equation (112) is a generalized eigenvalue equation, having in the generic case four different complex eigenvalues µ i with corresponding eigenvectors α i , defining four vectors that are product vectors in dimension 2 × 4 of the form Similarly, equation (114) gives four eigenvalues µ i with corresponding eigenvectors α i , defining four vectors that are product vectors when we use the split tensor product defined in Appendix F, Finally, equation (116) gives four eigenvalues µ i with corresponding eigenvectors α i , defining four vectors that are product vectors in dimension 4 × 2, If the two equations (112) and (116) have a common eigenvector α, then φ = ψα is a product vector in two ways, both in dimension 2 × 4 and in 4 × 2. This means that it is a product vector in dimension 2 × 2 × 2, and the same α is an eigenvector of equation (114).
Note that the standard method for solving the generalized eigenvalue equations (112), (114), and (116) depends on the nonsingularity of the matrices B, D, and F . If one or more of these matrices are singular, we may usually avoid the problem simply by making a random product transformation ψ → ψ = V ψ with V = V 1 ⊗ V 2 ⊗ V 3 , then solving the problem with ψ instead of ψ and transforming back in the end.

F The split tensor product
We find it useful to define a split tensor product so as to be able to take out the middle factor in a tensor product of three factors. Thus we define x ⊗ y ⊗ z = y ⊗ s (x ⊗ z) .