Robust Hadamard matrices, unistochastic rays in Birkhoff polytope and equi-entangled bases in composite spaces

We study a special class of (real or complex) robust Hadamard matrices, distinguished by the property that their projection onto a $2$-dimensional subspace forms a Hadamard matrix. It is shown that such a matrix of order $n$ exists, if there exists a skew Hadamard matrix of this size. This is the case for any even dimension $n\le 20$, and for these dimensions we demonstrate that a bistochastic matrix $B$ located at any ray of the Birkhoff polytope, (which joins the center of this body with any permutation matrix), is unistochastic. An explicit form of the corresponding unitary matrix $U$, such that $B_{ij}=|U_{ij}|^2$, is determined by a robust Hadamard matrix. These unitary matrices allow us to construct a family of orthogonal bases in the composed Hilbert space of order $n \times n$. Each basis consists of vectors with the same degree of entanglement and the constructed family interpolates between the product basis and the maximally entangled basis.


Introduction
Hadamard matrices with a particular structure attract a lot of attention, as their existence is related to several problems in combinatorics and mathematical physics. Usually one poses a question, whether for a given size n there exists a Hadamard matrix with a certain symmetry or satisfying some additional conditions. Such a search for structured Hadamard matrices can be posed as well in the standard case of real Hadamard matrices [16], or in the more general complex case [10,29,28].
To simplify investigations one studies equivalence classes of Hadamard matrices, which differ only by permutations of rows and columns and by multiplication by diagonal unitary matrices. It was shown by Haagerup [15] that for dimensions n = 2, 3 and n = 5 all complex Hadamard matrices are equivalent while for n = 4 there exists a one parameter family of inequivalent matrices.
In the recent paper by Karlsson [19], aimed to construct new families of complex Hadamard matrices of size n = 6, the author studied H 2 -reducible matrices defined as Hadamard matrices with the property that all their blocks of size 2 form Hadamard matrices.
In this work we introduce a related, but different notion of robust Hadamard matrices. This class contains a Hadamard matrix H n , such that each of its principal minors of order two forms a Hadamard matrix H 2 . The name refers to the fact that the Hadamard property is robust with respect to projection, as any projection of H n onto a two dimensional subspace spanned by the elements of the standard basis is a Hadamard matrix.
The notion of robust Hadamard matrices will be useful to broaden our understanding of the problem of unistochasticity inside the Birkhoff polytope [6]. For any bistochastic matrix B one asks whether there exist a unitary matrix U , such that B ij = |U ij | 2 . In the simplest case n = 2 any bistochastic matrix is unistochastic, for n = 3 the necessary and sufficient conditions for this property are known [3,17,23], in order n = 4 we provide an explicit construction in Appendix A, while for n ≥ 5 the unistochasticity problem remains open [1,5].
A robust Hadamard matrix H n will be used to to find unitary matrices of order n, designed to prove unistochasticity of bistochastic matrices located at any ray of the Birkhoff polytope. The key idea behind this construction -a decomposition of bistochastic matrix of an even dimension into square blocks of size two -was used in the algorithm of Haagerup to analyze the unistochasticity problem for n = 4. A family of unistochastic matrices joining a permutation matrix with the flat bistochastic matrix W n of van der Waerden allows us to construct a family of orthogonal bases in the composed Hilbert space H n ⊗ H n , which contains equi-entangled vectors and interpolates between the product basis and the maximally entangled basis [30]. Although such "equi-entangled bases" have already appeared in the literature [14,18], the present construction is simpler as it corresponds to straight lines in the Birkhoff polytope B n .
This work is organized as follows. In Section 2 we review the necessary notions, introduce robust Hadamard matrices and demonstrate their existence if skew Hadamard matrices or symmetric conference matrices exist. In Section 3 we present the unistochasticity problem and show how robust Hadamard matrices can be used to prove unistochasticity of bistochastic matrices at a ray joining the center of the Birkhoff polytope with any of its corners. The notion of complementary permutation matrices is used in Section 4 to demonstrate unistochasticity of certain subsets of the Birkhoff polytope of an even dimension n for which robust Hadamard matrices exist. Unitary matrices of an even order n ≤ 20 corresponding to rays in the Birkhoff polytope are applied to construct in Section 5 a family of equi-entangled bases for the composite systems of size n × n. In Section 6 we analyze the unistochasticity problem for n = 4: using the Haagerup procedure, presented in Appendix A, we investigate the geometry of the set U 4 of unistochastic matrices of order n = 4. An estimation of the relative volume of the set U 4 with respect to the standard Lebesgue measure is presented in Appendix B. In Appendix C it is shown that robust Hadamard matrices do not exists for n = 3 and n = 5, while extensions of these results for any odd n, suggested by the referee, are presented in Appendix D.

Robust Hadamard Matrices
We are going to discuss various subclasses of the set of Hadamard matrices. Let us start with basic definitions. Let M n (R) and M n (C) denote set of all real (complex) matrices of order n.
In analogy to the standard notion of (real) Hadamard matrices one defines [10,29] complex Hadamard matrices.
The name refers to the fact that the Hadamard property is robust with respect to projections: H R remains Hadamard after a projection Π 2 onto a subspace spanned by two vectors of the basis used, Π 2 H R Π 2 ∈ H(2).
Robust Hadamard matrices of order n will be denoted by H R n or H R . Equivalent definition is: Robust Hadamard matrix is a Hadamard matrix, in which all principal minors of order two are extremal so that their modulus is equal to 2. The notion of Hadamard matrices robust with respect to projection can also be used in the complex case.

Equivalence relations
Let P(n) or P denote the set of all permutation matrices P of order n. Permutation matrices are used to introduce the following equivalence relation in the set of Hadamard matrices. Definition 2.3. Two Hadamard matrices are called equivalent, written H 1 ∼ H 2 , if there exist unitary diagonal matrices D 1 , D 2 and permutation matrices P 1 , P 2 ∈ P such that (2.1) In dimensions 1, 2, 4, 8 and 12 all real Hadamard matrices are equivalent. For n = 16 there exist five equivalence classes [16] and this number grows fast [24] with n. It is also convenient to distinguish a finer notion of equivalence with respect to signs. Definition 2.4. Two Hadamard matrices are equivalent with respect to signs, written H 1 ≈ H 2 , if one is obtained from the other by negating rows or columns, Proof. To show this observe that by multiplying the columns of a robust Hadamard matrix H R which contain negative entries at the diagonal by (−1) we obtain a skew matrix. Note that multiplication of columns does not take H R out of the set of robust Hadamard matrices.

Basic properties of robust Hadamard matrices
One can easily see that the following Hadamard matrix of order four is robust and is equivalent to a skew Hadamard: To show that the set of H 2 -reducible matrices introduced by Karlsson [19] differs from the set H R (n) of robust Hadamard matrices discussed here, consider the following matrix of size n = 4, It is easy to see that the above matrix is H 2 -reducible but is not robust, so these two sets do differ.
2.4. Symmetric conference matrices and robust Hadamard matrices Definition 2.8. A symmetric matrix with entries ±1 outside the diagonal and 0 at the diagonal, which satisfy orthogonality relations: is called a symmetric conference matrix. Lemma 2.9. Matrices of the structure: where C is symmetric conference matrix, are complex robust Hadamard matrices.
This statement holds as the determinant of any principal submatrix of H R or order two is equal to −2.
To establish whether for a given order n there exists a robust Hadamard matrix it is therefore sufficient to find a skew Hadamard matrix or a symmetric conference matrix of size n. Concerning odd dimensions we show in Appendix C and D that robust Hadamard matrices do not exist for any odd n.

Robust Hadamard matrices and unistochasticity
Robust Hadamard matrices, introduced in the previous Section, will be used to investigate the problem, whether a given bistochastic matrix is unistochastic. Before demonstrating such an application let us recall the necessary notions.

Definitions
Definition 3.1. The set B n of bistochastic matrices of size n consists of matrices with non-negative entries which satisfy two sum conditions: Due to the Birkhoff theorem [6], the set B n is equal to the convex hull of all permutation matrices of order n. If matrix U is orthogonal the matrix B is called orthostochastic. The sets of all unistochastic and orthostochastic matrices of order n will be denoted by U n and O n , respectively. Above definitions imply:

In search for unistochasticity
It is easy to see that for n = 2 all three sets coincide, O 2 = U 2 = B 2 . While conditions for unistochasticity are known [3,11,12] for n = 3, the case n ≥ 4 remains open -see [5]. Given an arbitrary bistochastic matrix B of order n, it is easy to verify whether certain necessary conditions for unistochasticity are satisfied, but in general, no universal sufficient conditions are known and one has to rely on numerical techniques [5,13].
Since we are looking for a unitary U such that B = U •Ū , the absolute value of each entry is fixed, |U ij | = B ij , and one can investigate constraints implied by the unitarity. The scalar product of any two different columns of U is given by a sum of n complex numbers of the form U ijŪij . This sum can vanish if the largest modulus is not larger than the half of the sum of all n moduli. For n = 3 this condition is equivalent to the triangle inequality. The above observation, allows one to obtain a set of necessary conditions, a bistochastic B must satisfy to be unistochastic [26], related to orthogonality of columns, and rows of U , We shall refer to the above relations as polygon inequalities or "chain" conditions: the longest link of a closed chain cannot be longer than the sum of all other links. For matrices of size n = 3 these necessary conditions are known to be sufficient for unistochasticity [3,23], while B is orthostochastic if the bounds are saturated. Furthermore, in this dimension all chain conditions are equivalent -if inequality (3.1) is satisfied for any pair of two columns it will also be satisfied for the remaining pairs of columns and vectors [17,11]. The case n = 4 occurs to be more complicated as the chain inequalities are not sufficient [1,5], and examples of bistochastic matrices are known for which these condition are satisfied for some pairs of columns or rows and not satisfied by the other pairs [25,31]. An algorithm allowing one to establish numerically whether a given bistochastic matrix B of order n = 4 is unistochastic is described in the Appendix A.
Not being able to solve the unistochasticity problem in its full extent, we shall consider particular subsets of the Birkhoff polytope B n of an arbitrary dimension n. Figure 1 presents a sketch of the set of bistochastic matrices visualizing the problems considered. FIGURE 1. Sketch of the Birkhoff polytope B n -a set of (n − 1) 2 dimensions: the flat matrix W n in the center, permutation matrices P i at the corners. The corresponding rays R and counter-rays R meet at W n . The counter-ray R ends at the counter-permutation matrixP 1 Bold lines denote complementary edges, while dashed lines represent non-complementary edges. Region in gray -the triangle ∆P 1 , P 2 , W n -represents sets proved to be unistochastic. Definition 3.3. Bistochastic matrix W n , in which every element is equal to 1/n, is called the flat matrix.
Definition 3.4. A one-dimensional set of bistochastic matrices obtained by a convex combination of the flat matrix W n and any permutation matrix P is called a bistochastic ray, Definition 3.5. A set R n of bistochastic matrices belonging to the line joining the flat matrix W n and a permutation matrix P , outside the segment P W n is called a counter-ray. These matrices can be expressed as pseudo-mixtures (3.3) with a negative weight α ∈ [− 1 n−1 , 0). For P = I, the family of matrices belonging to the ray is as follows 3.3. Robust Hadamard matrices imply unistochasticity Lemma 3.6. If there exists a robust complex Hadamard matrix H R of order n then all rays and counter-rays of the Birkhoff polytope B n of order n are unistochastic.
Proof. First let us show, that for any robust Hadamard matrix H R the following property holds: where D is the diagonal matrix containing diagonal entries of H R . Diagonal elements of the lefthand side of eq. (3.5) are equal to 2H ii H * ii = 2, while off-diagonal entries of this sum, i = j, read H ij H jj * + H ji * H ii . In order to evaluate these terms we shall use the fact that the matrix Let M 2 be the principal submatrix of H R D * of order two spanned by the rows i and j. Since H R D * is robust then M 2 ∈ H C (2), so that M 2 M * 2 = 2I 2 . Writing down the entries of this matrix explicitly one obtains Since this matrix is proportional to identity its off-diagonal entries vanish, so that For any bistochastic matrix R of order n belonging to the ray R I or counter ray R let us now construct a matrix U , where real parameters a and b defined by eq. (3.4) satisfy the condition (n − 1)b + a = 1. Making use of these relations, definition (2.1) and eq. (3.5) we get This shows that the matrix U is unitary. Since the matrix R satisfies the relation R ij = |U ij | 2 , it is unistochastic. Hence any matrix R at any ray R of the Birkhoff polytope B n or any counter-ray R is unistochastic for any dimension n, for which a robust Hadamard matrix exists.
In particular, if the robust Hadamard matrix is real, so that U becomes orthogonal, then the matrix R is orthostochastic. Since every skew Hadamard matrix is robust (lemma 2.6), we arrive at the following statements: Proposition 3.7. For every order n, for which there exists a symmetric conference matrix, every matrix belonging to any ray R or any counter-ray R of the Birkhoff polytope B n is unistochastic.
Proposition 3.8. For any order n, for which there exists a skew Hadamard matrix H S , every matrix belonging to any ray R or any counter-ray R of the Birkhoff polytope B n is orthostochastic.
Existence of skew Hadamard matrices for orders n = 4k is proved for k < 69, proper construction was done by Paley. There are infinitely many cases of skew Hadamard matrices of higher orders [21].
It is known that for dimensions n = 6, 10, 14, 18 there exists a symmetric conference matrix [22]. However, for order n = 22 there are no such matrices, since 21 is not the sum of two squares.
Those facts imply the main result of this work: Proof. It is a simple conclusion from proposition 3.7 and proposition 3.8 which allow one for an explicit construction of the corresponding orthogonal and unitary matrices.
Note that the above statement holds also for infinitely many dimensions, for which symmetric conference matrices are known (first found by Paley -see e.g. [16]).

Unistochasticity of certain triangles embedded inside Birkhoff polytope
A convex combination of any two permutation matrices forms an edge (or a diagonal) of the Birkhoff polytope B n , E = αP + (1 − α)Q Au-Yeng and Cheng introduced the notion of complementary permutations [2]. Definition 4.1. Let P = (P ij ) and Q = (Q ij ) be two permutation matrices. Then the matrices P and Q are called complementary if equality P ij = P hk = Q ik = 1 implies Q hj = 1 for all i, j, h, k ∈ {1, ..., n}; (consequently, if Q ij = Q hk = P ik = 1 then P hj = 1).
Proposition 4.2. If P and Q are two complementary permutation matrices then the entire edge P Q of the Birkhoff polytope B n connecting P and Q is orthostochastic. If P and Q are not complementary, then the edge is not unistochastic, besides the orthostochastic corners P and Q.
Below we generalize above proposition, proved by Au-Yeng and Cheng [2] to establish unistochasticity of some sets of larger dimension. Let us first distinguish the following notion: In other words their non-zero elements are put in different places, so their Hadamard product vanishes, P • Q = 0. Due to symmetry of B n one can take the identity matrix for the permutation matrix P without loosing generality. A matrix Q is strongly complementary to P = I if and only if Q is an involution, Q 2 = I, and every diagonal element of Q is 0. Due to complementarity, the dimension of Q is even. Making use of the flat bistochastic matrix W n we can now formulate statements concerning triangles ∆(A, B, C) belonging to B n and spanned by vertices A, B and C. This statement is visualized by a gray triangle of unistochastic matrices shown in Fig. 1. The above result can be generalized for a larger set of permutation matrices {P 1 , P 2 , ..., P k } of dimension n, in which each pair of matrices is strongly complementary. Then an analogous reasoning shows that the bistochastic matrices belonging to 2-faces of the polytope defined by the convex hull of k these permutation matrices and the flat matrix W n are unistochastic.
Due to the definition of strong complementarity each pair of matrices has non-zero entries in different places, so that k ≤ n. We are not able to determine how the maximal number k of such matrices depends on the dimension n, nor whether the bistochastic matrices belonging to the interior of this polytope are unistochastic.

Equi-entangled bases
Any quantum system composed from two subsystems, with n levels each, can be described in a Hilbert space with a tensor product structure. It is often natural to use the standard, product basis, usually denoted by |k, m = |k ⊗ |m with k, m = 1, . . . , n. However, for certain problems it is advantageous to use bases consisting of maximally entangled states. In the simplest case of n = 2 one uses the Bell basis consisting of four orthogonal states of size four, |Ψ ± = (|0, 0 ± |1, 1 )/ √ 2 and |Φ ± = (|0, 1 ± |1, 0 )/ √ 2. These Bell states are also called maximally entangled, as their partial traces are maximally mixed. In this Section we follow the notation of [30], often used in quantum theory.
For any higher dimension n such entangled bases were constructed by Werner [30]. A slightly more general variant of this construction discussed in [4] allows one to write a set of n 2 orthogonal states related to a given unitary matrix U of order n. Making use of unistochastic matrices belonging at a ray of the Birkhoff polytope and determined by robust Hadamard matrices, we shall construct a family of bases interpolating between separable and maximally entangled basis.
Consider a unitary matrix U of order n, associated with the unistochastic matrix B = U •Ū , which belongs to the ray R of the Birkhoff polytope B n -see Fig. 1. It exists for all dimensions, for which a robust Hadamard matrix exists (i.e. for all even dimensions n ≤ 20) and allows us to construct an equi-entangled orthogonal basis in the Hilbert space describing a composed system of size n × n. Proof. a) To show this property it is sufficient to check the scalar product Due to unitarity of U the orthogonality holds for any k, k , m, m = 1, . . . , n, which implies that the vectors (5.1) form an orthonormal basis. b) To analyze the degree of entanglement we specify our considerations to a family U (a) of unitary matrices given in (3.4) and parameterized by a single parameter a. A state |ψ m,k defined in eq. (5.1) can be rewritten in a slightly different way, where |j 1 is an element of the computational basis in the first system, while |j 2 = |j + k |mod n is obtained by relabeling the second basis. As any quantum state is defined up to a complex phase we are allowed to replace the complex prefactor U mj by its modulus |U mj |. Note that the above form can be interpreted as the Schmidt decomposition of the bipartite state, |ψ AB = n j=1 λ j |j A ⊗ |j B . Observe that the components of the Schmidt vector, λ j = |U mj | 2 , form a row of the bistochastic matrix B a belonging to a ray of the Birkhoff polytope B n . For each basis state |ψ m,k its ordered Schmidt vector is the same, λ = (a, b, b, b, . . . , b) with b = (1 − a)/(n − 1), so that all measures of entanglement of the states |ψ m,k are equal. Therefore the orthonormal basis parameterized by a number a ∈ [1/n, 1] forms a family of equi-entangled bases, which interpolate between a maximally entangled basis (a = 1/n) and a separable basis (a = 1).
To characterize the degree of entanglement one can apply the entanglement entropy of a basis state, equal to the Shannon entropy of the corresponding Schmidt vector and to the von Neumann entropy of the reduced density matrix. As the Schmidt vector of each state has the structure λ = (a, b, b, b, . . . , b) its entropy reads and yields the entanglement entropy of the basis states (5.3). It is equal to zero for a = 1 and a separable basis and it achieves the maximum equal to log n for a = 1/n, which corresponds to the maximally entangled basis, analyzed in [30].
6. Set U 4 of unistochastic matrices of order n = 4 For even dimensions n ≤ 20 we have shown that all rays of the Birkhoff polytope are unistochastic. This statement holds also in the case n = 4 -the smallest dimension for which the unistochasticity problem is still open [1]. In this case the chain conditions (3.1) and (3.2) are necessary, but in contrast to the case n = 3 they are not sufficient to assure unistochasticity [5]. Although analytic form of such conditions remains still not known, we shall apply a numerical procedure proposed by Haagerup -see Appendix A -to study the properties of the set U 4 of unistochastic matrices of order n = 4. Generating random bistochastic matrices according to the flat measure in the Birkhoff polytope B 4 according to the algorithm described in [7] we found that the relative volume of the set C 4 of bistochastic matrices satisfying all chain conditions (3.1) and (3.2) is V C ≈ 0.71, while the volume of its subset U 4 containing unistochastic matrices is V U ≈ 0.61 in comparison to the total volume of B 4 .
However, since not much is known about geometric properties of the subsets U 4 ⊂ C 4 ⊂ B 4 of the Birkhoff polytope, we shall study various cross-sections which include the flat matrix W 4 located at its center.
Any such cross-section is determined by three matrices which do not belong to a single line. Since one of these matrices is selected to be W 4 , we have to specify only two other matrices. Crosssections shown in Fig. 2 are specified by two permutation matrices. Due to the symmetry of the polytope B 4 we may take the first one as identity, I 4 , without loosing the generality.
There Cross-sections determined by unistochastic edges of length √ 8 are simple as all bistochastic matrices are unistochastic, which follows from Section 4. Below we present cross-sections determined by edges from remaining three classes, created by the following permutation matrices (the numbers in the subscripts denote the cycles):  a) unistochastic short edge I ↔ P A of length 2 b) not unistochastic middle edge I ↔ P B of length √ 6 c) not unistochastic long edge I ↔ P C of length √ 8 Darkest gray (in color: dark blue) represents unistochastic set U 4 , medium gray (red) represents the set C 4 of matrices that satisfy the chain conditions, and the lightest gray (light blue) denotes bistochastic matrices, that do not satisfy chain conditions. At the section showed in panel c) the sets U 4 and C 4 do coincide, while the bistochastic matrix at the edge reads B c = 2W 4 − 1 2 P c − 1 2 I 4 . Note the dark sets U 4 in all panels containing both counter-rays are not convex but are star-shaped.
To produce each plot we generated a lattice of around 10 4 bistochastic matrices belonging to a given 2D cross-section and verified if conditions (3.1) are satisfied and whether numerical procedure described in Appendix A returns the corresponding unitary matrix. The Birkhoff polytope B 4 possesses an interesting property -in every neighborhood of the flat unistochastic matrix W 4 localized at its center there are non-unistochastic matrices. Hence there is a direction, in which deviation by an arbitrary small leads to a matrix which is not unistochastic [5]. Selecting one of those directions according to [5] we find a cross-section shown in Fig. 3, such that non-unistochastic matrices are located arbitrary close to W 4 : FIGURE 3. Cross-section of the Birkhoff polytope B 4 , determined by the directions V 1 and V 2 described in [5]. The flat matrix W 4 in the center of the B 4 is also localized at the boundary of the set U 4 . Colors have the same meaning as in Fig.  2.
Our findings do not contradict the conjecture, that the set U 4 of unistochastic matrices is starshaped with respect to W 4 . There exist 2D cross-sections, in which the sets U 4 and C 4 coincide, but it is easy to find a cross-section which reveals that the inclusion relation U 4 ⊂ C 4 is proper. Compare the discussion of relative volumes of these sets presented in Appendix B.

Concluding remarks
A notion of Hadamard matrices robust with respect to a projection onto a subspace formed by any two of the basis vectors is introduced. In the case of a double even dimension, n = 4k, existence of such matrices follows from existence of skew real Hadamard matrices. Furthermore, existence of symmetric conference matrices of dimension n = 4k + 2 with k ≥ 1 yields a construction of a robust complex Hadamard matrix of this size. This implies that robust Hadamard matrices exist for all even dimensions n ≤ 20.
Existence of robust Hadamard matrices of order n allows us to show that all the rays of the Birkhoff polytope B n of bistochastic matrices are unistochastic. Hence for any point B a at the line joining an edge (a permutation matrix) with the flat matrix W n at the center of the polytope, there exists a corresponding unitary matrix U (a) of size n such that B a = U (a) •Ū (a). Furthermore, if two permutation matrices P and Q are strongly complementary [2] -see Definition 4.3 -the entire triangle ∆(P, Q, W n ) is unistochastic.
The family of unitary matrices U (a) of order n obtained with help of robust Hadamard matrices of size n, allows us to construct a family of equi-entangled bases in the composed Hilbert space of n × n system. This family interpolates between the separable basis and the maximally entangled basis. Each interpolating basis consists of n 2 normalized vectors with the same Schmidt vectors, which determine the entropy of entanglement (5.4). In contrast to the earlier constructions given in [18] and extended in [14] our construction is based on a straight line in the Birkhoff polytope B n , so the Schmidt vectors contain only two different entries.
This work contributes to investigations [5] of the set U 4 of unistochastic matrices of order n = 4. Making use of the Haagerup algorithm to verify numerically whether a given bistochastic matrix B ∈ B 4 is unistochastic -see Appendix A -we provided a first estimation of the relative volume of the set U 4 of unistochastic matrices and the larger set C 4 of these bistochastic matrices, which satisfy all chain conditions (3.1) and (3.2). Furthermore, we studied the geometry of crosssections of the set U 4 along planes determined by selected permutation matrices -corners of the Birkhoff polytope In this appendix we present a method due to the late Uffe Haagerup of searching for a unitary matrix corresponding to any unistochastic matrix of order n = 4. Given a bistochastic matrix B ∈ B 4 we wish to find U such that B ij = |U ij | 2 .
Any matrix of order four can be decomposed into four blocks of order two. Let us then represent in this way the unitary matrix U we are looking for, Note that all the moduli are determined by the bistochastic matrix B, so only the phases remain unknown. Unitarity of U implies that Due to existence of the Hadamard matrix H 4 we know that the flat matrix W 4 is orthostochastic, so we can assume that B = W 4 . Permuting rows and columns of B one can rearrange the matrix in such a way that the following relation holds After this is done the norm of the block A is bounded, ||A|| 2 HS < 1 ⇒ eigenvalue(AA * ) < 1 ⇒ eigenvalue(XX * ) > 0 , so the neighboring block X is invertible. Hence we can use the matrix X −1 to transform the second unitarity condition in (7.2) into an explicit expression for the lower diagonal block D, The next step is to find the phases in the blocks X and Y . To take into account orthogonality between the first two rows of the matrix U , it is convenient to introduce four auxiliary variables, They allow us to rewrite this orthogonality condition, which can be treated as an equation for the unknown phases. Orthogonality requires that a chain formed out of four links of lengths l i has to be closed, so the longest link is not longer than the sum of remaining links, If this condition is not satisfied, the matrix B is clearly not unistochastic. Observe that these conditions can be interpreted as particular cases of the general conditions (3.1) and (3.2).
Making use of the orthogonality condition between two first columns of U we arrive at an equation analogous to (7.4), which for a given phase φ can be solved for unknown phases β 1 and β 2 . This determines the blocks X and Y and allows us to obtain the remaining block D of U .
The catch is that the explicit formula (7.3) produces a matrix D, which needs not to be compatible with the structure imposed by the initial bistochastic matrix B. To find the desired unitary matrix U it is sufficient to check that a single element of D has the correct norm. Hence we arrive at the following criterion for unistochasticity for a matrix B of order four: A bistochastic matrix B ∈ B 4 is unistochastic if there exist a phase φ entering eq. (7.1), (which determines phases α i and β i and thus blocks X and Y ), such that the block D obtained by eq. (7.3) satisfies the constraint |D 11 | 2 = B 33 .
If this is the case the unitary matrix U given by the form (7.1) satisfies the unistochasticity condition, conveniently written by the Hadamard product, B = U •Ū . Note that the above procedure can be easily implemented numerically for any given bistochastic matrix B = W 4 of order four.

Appendix B. Unistochasticity and chain link conditions in small dimensions -numerical results
Random points in a cube. The simplest method to generate a random bistochastic matrix B from the Birkhoff polytope B n with respect to the uniform measure is to take random numbers, uniformly distributed in [0, 1] for each entry of the core -the principal submatrix of order n − 1, which determines the matrix B. If the sum of any row or column of the core is greater than 1 or the sum of all elements of the core is smaller than n−2, then the generated matrix is not bistochastic and will not be considered. If all these conditions are met, we generate a bistochastic matrix by filling the missing row and column with values which will add up all rows and columns to 1. The entry B nn is equal to the sum of all elements of the core reduced by n − 1.
To optimize computational procedure of generating the cores, we use row discrimination. Drawing the core row-by-row (rows are independent) allows one to check the sum of the row at every step and stop the procedure if it exceeds unity. This method occurs to work efficiently for dimensions n ≤ 6, so for higher dimensions one needs to apply other techniques [7,9].
Sinkhorn algorithm To generate a bistochastic matrix for n ≥ 7 we also used a method described in [7], analogous to the Sinkhorn algorithm [27], which normalizes rows and columns of a given square matrix in a following sequence: 1. input a random stochastic matrix of dimension n with positive elements, 2. normalize every row by dividing it by the sum of its elements, 3. normalize every column of the matrix by dividing it by the sum of its elements, 4. go to point 2., unless the matrix is bistochastic up to a certain accuracy with respect to the chosen norm.
In practice we stopped the procedure if all the sums of rows and vectors are close to unity up to the sixth decimal place (0.9999999 ≤ Σrows, Σcolumns ≤ 1.0000001). This procedure occurs to be much faster than taking points at random from the core for dimensions n > 6.
Dirichlet distribution. It is natural to study the measure induced in B n by the Sinkhorn algorithm applied to random stochastic matrices with entries described by the Dirichlet distribution [20], P s (x 1 , . . . , x n ) ∝ n i=1 x 1−s i , with the constraint i x i = 1. For s = 1 this distribution coincides with the uniform distribution in the probability simplex. For a certain value of the parameter s * = 1 2n 2 (n 2 − 2 + √ n 4 − 4) the measure induced in the set of bistochastic matrices becomes uniform [7] in the limit of large n.
To generate a sequence of n numbers distributed according to Dirichlet distribution with a chosen parameter s we used the following sequence: 1. generate n independent random numbers from the gamma distribution of rate 1 and shape s.
Explicitly, every element of (x 1 , x 2 , . . . , x n ), is drawn with probability density 1 Γ(s) x s−1 e −x . 2. normalize every element by dividing it by the sum of all elements, Applying the above procedure n times we obtain random stochastic matrix of n independent rows each distributed according to Dirichlet distribution. Making use of the Sinkhorn algorithm we to obtain an ensemble of bistochastic matrices which is uniformly distributed at the center of the Birkhoff polytope.
Numerical computations show that the first method (generating random points in a cube) is reliable for n ≤ 6, while for n > 6 the Sinkhorn method, used with the Dirichlet parameter s * , becomes more efficient. A numerical estimation of the relative volume of the set of unistochastic matrices is shown in Fig. 4. Note that the relative volume of the set C n of matrices satisfying the chain conditions approaches unity.
In this Appendix we provide results concerning existence of robust Hadamard matrices which were suggested by the referee and go beyond the statements formulated in Appendix C.
Let C be a complex conference matrix of order n, that is C jj = 0, |C jk | = 1 and CC * = (n − 1)I. Lemma 7.1. If R is a robust Hadamard matrix of order n ≥ 2, then R is equivalent in the sense of eq. (2.4) to a matrix H = C + iI, where C is a self-adjoint complex conference matrix.
Proof. The relation between R and H is explicitly: where D is the diagonal matrix containing diagonal entries of R and i is imaginary unit. All of the diagonal elements of H are equal to i. Because D is a unitary matrix, then R ∼ H and thus H is a robust Hadamard matrix. This implies that any principal submatrix of H of order two, has to satisfy |det(M 2 )| = 2. Because |h jk | = 1 for any j, k, then h jk = h * kj . It implies that the matrix consisting of off-diagonal elements of H, namely C = H − iI is self-adjoint (C = C * ) and |C jk | = 1 for j = k. Self-adjointness of C, it implies that H − H * = 2iI. Above result shows that robust complex Hadamard matrices are essentially a family of matrices equivalent to self-adjoint complex conference matrices with the diagonal elements filled with imaginary unit i. Proof. Applying Lemma 7.1 we can restrict ourselves to the case H = C + iI. Since C = C * and C is unitary up to scalar, then its spectrum is real and the eigenvalues are equal to ± √ n − 1. Note also that C is traceless. Therefore, for n ≥ 2, the number of positive and negative eigenvalues must be equal. It follows that the number of eigenvalues, which is equal to n, must be even.