Is there another parameterization of a qutrit 3-level system, besides Gell-Mann?

Question

Question: Is there a parameterization of a a general qutrit 3-level system similar to: $$\rho = \begin{bmatrix} p_1 & r_{12}e^{-i\phi_{12}} & r_{13}e^{-i\phi_{13}}\\ \cdot & p_2 & r_{23}e^{-i\phi_{23}}\\ \cdot & \cdot & 1-p_1-p_2 \end{bmatrix} \tag{1.1}$$ With the constraints: $$\begin{align} & 0\leq p_j \\ & p_1+p_2\leq1 \\ & 0\leq r_{ij}\leq \sqrt{p_ip_j} \\ & 0 \leq \phi_{ij} < 2\pi \end{align} \tag{1.2}$$ ?

Today I learned a 3-level system can be parameterized, from Gell-Mann matrices, as:

$$ \rho = \begin{bmatrix} \frac13+a_3+\frac{a_8}{\sqrt 3} & a_1-ia_2 &a_4-ia_5 \\ \cdot & \frac13-a_3+\frac{a_8}{\sqrt 3} & a_6-ia_7 \\ \cdot &\cdot & \frac13-\frac{2a_8}{\sqrt 3} \end{bmatrix} \tag{2.1}$$ Subject to the constraints: $$\begin{align} & 0\leq \frac13\pm a_3+\frac{a_8}{\sqrt 3}\leq1 \\ & 0\leq \frac13-\frac{2a_8}{\sqrt 3}\leq1 \\ & \mathbf a\cdot \mathbf a\leq\frac13 \\ & 0 \leq \det(\rho) \end{align} \tag{2.2}$$

Disclaimer: I know parameterization $(1)$ is not correct, hence I am asking if there is anything similar. In particular, my problem with the Gell-Mann parameterization $(2)$ is that I have to compute a determinant each time I draw 8 random numbers... I would like something faster to generate and validate.

---

Background: I am doing numerical simulations, where I need to quickly generate a large number of random 3-level density operators. I was using parameterization $(2)$, with which you can easily draw the eight degrees of freedom of the qutrit.

However, the parameterization is wrong since it allows non-physical states, such as: $$\rho = \begin{bmatrix} p_1 &\sqrt{p_1 p_2} &\sqrt{p_1 p_3}\\ \cdot &p_2 &\sqrt{p_2 p_3}e^{-i\pi}\\ \cdot &\cdot &p_3 \end{bmatrix}$$ which is not positive semidefinite.

Answer 1

First, note that the direct answer to the question "Is there another parameterization of a qutrit 3-level system, besides Gell-Mann?" is a pretty obvious yes. There isn't anything intrinsically special about using Gell-Mann matrices to decompose qutrits, it's a basis of operators like many others (well, it has some symmetries that make it interesting in some areas, as discussed e.g. here, but this doesn't matter much for the purpose of this question). But of course, parametrising state space in dimensions larger than two is tricky due to the boundary not having a simple spherical shape, which makes the constraints on any such parametrisation a potentially delicate matter. Some of this has already been covered in What are the conditions ensuring a two-qubit density matrix is positive semidefinite?, and links therein. Here I'll try to focus a bit more on the parametrisation you mention and the associated semipositivity constraints.

Let's consider $d$-dimensional states. A very general answer to the question is that you can take an arbitrary orthonormal basis of Hermitian operators $\sigma_i\in\operatorname{Herm}(\mathbb{C}^d)$. Any state $\rho$ is then represented (or "parametrised", if you prefer) by the vector of coefficients in the decomposition $\rho=\sum_{i=1}^{d^2-1}c_i(\rho)\sigma_i$. Note that here the "orthonormality" means we're requiring $\operatorname{tr}(\sigma_i\sigma_j)=\delta_{ij}$. It is also common to use operators normalised as $\operatorname{tr}(\sigma_i\sigma_j)=d\delta_{ij}$, but this doesn't change much in the discussion, so I'll stick with actual orthonormal operator bases here.

Given the normalisation constraint on density matrices, it is always convenient to include the identity operator in the basis. Rescaling to ensure normalisation, we can thus assume without (much) loss of generality that $\sigma_0=I/\sqrt d$, and then orthogonality implies that all $\sigma_i$ with $i>0$ are traceless. Then, for any such basis, you can characterise the boundary of state space in "hyperspherical coordinates" with the constraint $$\|\vec c(\rho)\| \le \frac{1}{|\min_k \lambda_k(\rho)|},$$ where $\lambda_k(\rho)$ is the $k$-th eigenvalue of $\rho$. You can find a proof of this e.g. here.

Regarding a parametrisation like your (1.1), in terms of probabilities and phases, consider rewriting it as $$\rho = \begin{bmatrix} p_1 & \sqrt{p_1 p_2} r_{12}e^{-i\phi_{12}} & \sqrt{p_1 p_3 }r_{13}e^{-i\phi_{13}}\\ \cdot & p_2 & \sqrt{p_2 p_3}r_{23}e^{-i\phi_{23}}\\ \cdot & \cdot & 1-p_1-p_2 \end{bmatrix}, $$ where all parameters are real. Sylvester's criterion now tells you that for positive semidefiniteness is sufficient to check all the principal minors. The conditions associated to 1x1 submatrices amount to $(p_1,p_2,1-p_1-p_2)$ being a probability vector, that is, $0\le p_1,p_2,p_1+p_2\le 1$. The conditions associated to 2x2 submatrices give you $0\le r_{ij}\le 1$. Finally, the condition on the full determinant can be put in the form $$\|\vec r\|^2 \le 1 + 2r_{12} r_{13} r_{23} \cos(\phi_{13}-\phi_{12}-\phi_{23}).$$ Here $\|\vec r\|^2\equiv r_{12}^2+r_{23}^2+r_{13}^2$. Note how this choice of parameterization with the $\sqrt{p_i p_j}$ factors gives conditions on the rest of the parameters which are independent on the probabilities themselves. This kind of thing probably works also in higher dimensions, and maybe you can write the associated general expressions, but it might get quite a bit trickier, and I'd guess you'll have to keep track of how different permutations classes contribute to the different determinants. I haven't actually tried.

Another advantage of this type of representation is that it gives a particularly simple expression for the purity: $$\operatorname{tr}(\rho^2) = \sum_{j=1}^3 p_j^2 + 2\sum_{j<k} p_j p_k r_{jk}^2 = 1 + 2\sum_{j<k} p_j p_k (r_{jk}^2-1),$$ and thus the phases $\phi_{jk}$ do not affect the purity at all, while $r_{jk}=1$ ensures pure states, because you can then complete the square.

Answer 2

Updated answer after DaftWullie's comment:

I think I over-complicated; if $U$ is a random unitary in $SU(3)$, $U\vert0\rangle$ is good enough to get a random pure state, so you could get a random pure density matrix with

$$\rho = U\vert0\rangle\langle0\vert U^\dagger$$

and then convexly combine pure density matrices to get mixed ones.

On the other hand, what DaftWullie was really suggesting seems even less complicated, and so more reasonable; just pick a random basis that diagonalizes the matrix.

$$\rho = U (a \vert{0}\rangle\langle{0}\vert + b\vert{1}\rangle\langle{1}\vert + (1-a-b)\vert{2}\rangle\langle{2}\vert) U^\dagger$$

and the comments below about the parametrization of random unitaries still apply.

Old answer:

I don't have enough reputation to comment, unfortunately, so consider what's below as an answer stub, rather than a fully fledged answer, though I think it might be a good first step.

As DaftWullie mentioned, you can get a pure density matrix more-or-less easily, by choosing suitable amplitudes in a random basis; thus if you have the ability to generate a Haar-random unitary matrix $U$, you can do

$$ \rho = U (a \vert0\rangle + be^{i\phi_1}\vert1\rangle + e^{i\phi_2}\sqrt{1-a^2-b^2}\vert2\rangle) \cdot \mbox{c.c.} $$

for $a,b \in [-1,1]$ satisfying $a^2 + b^2 \leq 1$.

Parametrizing the space of unitaries seems doable, but not easy; see this answer.

Once you can generate random pure density matrices, you can make use of the fact that the density matrix space is convex; choose two matrices and interpolate between them. I doubt this will give you a uniform distribution though, simply based on the intuition that getting the fully mixed matrix will require a very specific choice of vertices for the interpolation.

I hope someone can build on this and help you further.