What are the conditions ensuring a two-qubit density matrix is positive semidefinite?

Question

I've seen some papers writing $$\rho=\frac{1}{4}\left(\mathbb{I} \otimes \mathbb{I}+\sum_{k=1}^{3} a_{k} \sigma_{k} \otimes \mathbb{I}+\sum_{l=1}^{3} b_{l} \mathbb{I} \otimes \sigma_{l}+\sum_{k, l=1}^{3} E_{k l} \sigma_{k} \otimes \sigma_{l}\right).\tag1$$ I wonder what condition should the matrix $E$ obey?

For a one-qubit state, the density matrix satisfies $\newcommand{\tr}{{\operatorname{tr}}} \tr(\rho)=1$ and is positive semidefinite. And the general form is $$\frac12\begin{pmatrix}1+z & x-iy\\x+iy & 1-z\end{pmatrix},$$ satisfying the trace condition. As for positive semidefinite conditions, it is well known that a Hermitian matrix is positive semidefinite iff its eigenvalues are non-negative. So I calculate the eigenvalues of it and get the restriction that $x^2+y^2+z^2 \le 1$, which can be seen as exactly the Bloch sphere.

Then I want to see the same thing happens in two qubits case. But I can only mimic the same reasoning and get the general form from Eq. (1). And when I try to calculate the eigenvalues of this matrix, even Mathematica showed a complicated result. But if we think about the separable case, it easy to see that the vector $a$ and vector $b$ should have length less than 1. But I can't find the restriction on matrix $E$. Hence my question:

What's the restriction on the parameters in Eq. (1) to make $\rho$ a valid density matrix, i.e., trace condition and positive semidefinite condition?

Answer 1

1. (General framework) The general form of this problem is the following. Take a generic Hermitian, unit-trace operator, $X\in\operatorname{Herm}(\mathcal X), \operatorname{Tr}(X)=1$, acting on some $N$-dimensional complex vector space $\mathcal X$. The set of such operators is an affine space of dimension $N^2-1$, embedded in the $N^2$-dimensional real vector space $\operatorname{Herm}(\mathcal X)$ of Hermitian operators. Any operator can be decomposed linearly with an operatorial basis. In particular, any unit-trace Hermitian operator can be decomposed linearly in an orthonormal basis of Hermitian traceless operators. This is the "generalised Bloch representation", which we can write as $$X = \frac{1}{N}\left(I + \sum_{k=1}^{N^2-1} c_k\sigma_k\right).$$ More precisely, we should write $c_k=c_k(X)$. These coefficients $c_k$ implement a (linear) isomorphism between the set of unit-trace Hermitian operators and the vector space $\mathbb R^{N^2-1}$.

2. (Positivity constraint) If we want a positive semidefinite unit-trace Hermitian operator, write it $\rho$, we need to impose further restrictions on the coefficients $c_k$. More specifically, rather than these spanning a full linear space, they will now only cover a convex subset of $\mathbb R^{N^2-1}$. This means, in particular, that you can describe the corresponding set of states via its boundary.

3. (Boundary ain't nice for $N>2$) This boundary is, in general, not very nice. For a single qubit you get a sphere, but as soon as you increase the dimensions, things get messier. For example, this boundary contains non-pure states (which doesn't happen for $N=2$). Similarly, even though all pure states still lie on a hypersphere (assuming a suitable parametrisation), they only cover some subsets of it.

4. (Radial characterisation is doable) Nonetheless, as shown here, it is possible to get a decent characterisation of this boundary, by using "radial" coordinates. This requires computing the eigenvalues corresponding to a given direction, that is, the eigenvalues corresponding to a given linear combination of basis operators. This is numerically trivial to do for any given direction (assuming $N$ is not too big, of course), but finding an analytical expression working for all directions would be very complicated.

   It is worth noting that this amounts to describing the (convex) set of states via its support function. The conversion between this description and a more standard analytical characterisation is quite complicated, if doable at all.

5. (So... what are the constraints?) Getting back to the question at hand: what does it really mean to find "restrictions for the coefficients $E_{jk}$"? Generally speaking, the boundary of the set of coefficients $c_k(\rho)$ corresponding to valid states $\rho$ is the solution set of $$f(c_1,...,c_{N^2-1})=0$$ for some function $f$. E.g. for $N=2$, we have $f(c_1,c_2,c_3)=c_1^2+c_2^2+c_3^2$. In higher dimensions, this $f$ is more complicated. For the case at hand, the coefficients are related to each other by some implicit relation $$f(a_1,a_2,a_3,b_1,b_2,b_3,E_{11},...,E_{33})=0.$$ This means that the viable range of any $E_{jk}$ depends on the value of all the other parameters. What you can do easily is find things like max and min of any parameter. This amounts to asking, e.g. what is the largest value of $E_{ij}$ that corresponds to any quantum state. Using the radial characterisation above, this can be seen to equal $1/|\lambda_{\rm min}(\sigma_i\otimes\sigma_j)|=1$. In fact, by similar reasoning, you can see that, using your operatorial basis built with products of Pauli matrices, every parameter lies in the interval $[-1,1]$.

   Note that this might seem to contradict the simple condition you get imposing $\operatorname{Tr}(\rho^2)\le1$, given in the other answer. That this is not a contradiction, is an interesting point per se. In fact, $\operatorname{Tr}(\rho^2)=1$ traces a spherical boundary that (1) surrounds the full set of states, and (2) is only touched by pure states. In other words, the answer $\max_\rho |E_{jk}|=1$ does't contradict the boundary given in the other answer, because the former is achieved by non-pure states (which is also easy to verify directly).

Answer 2

Since Pauli matrices are traceless and Hermitian, the trace condition $\text{Tr}(\rho)=1$ is satisfied automatically and the numbers $a_k, b_l, E_{kl}$ should be real if we want $\rho$ to be Hermitian. The hard part is to verify the positive semidefinite condition $\rho \ge 0$.

The straightforward way is to use Sylvester’s criterion.

The positivity condition $\rho > 0$ is equivalent to just 4 strict inequalities, representing the condition that leading principal minors must be positive. In our case, the first minor is $\frac{1}{4}(1 + a_z + b_z + E_{zz})$, the last one is the determinant of $\rho$. It's hard to make a simple meaning of obtained inequalities, in contrast to the single qubit case.

The non-strict positivity condition $\rho \ge 0$ is even harder to express since by Sylvester’s criterion we need all principal minors (4+6+4+1=15 in our case) to be non-negative.

A simple necessary condition you can get from the inequality $\text{Tr}(\rho^2)\le 1$. It gives $$\sum a_k^2 + \sum b_l^2 + \sum E_{kl}^2 \le 3.$$ Though, it should be deducible from those 15 inequalities anyway.