How to construct — based on physics-type notation — a magical simplex $\mathcal{W}$ of bipartite qutrits?

Question

I have a short Mathematica program:

S[0] = {1, 0, 0}; S[1] = {0, 1, 0}; S[2] = {0, 0, 1};

Omega[0, 0] = (1/Sqrt[3]) Sum[KroneckerProduct[S[s], S[s]], {s, 0, 2}]

w = Exp[2 Pi I/3];

Do[
 If[{k, l} != {0, 0}, 
  Omega[k, l] = 
   (1/Sqrt[3]) 
    Sum[KroneckerProduct[w^(k (s - l)) (S[s] - S[l]), S[s]], {s, 0, 2}] 
     Omega[0, 0]], 
 {k, 0, 2}, {l, 0, 2}]

Do[
 P[k, l] = 
  ArrayReshape[Outer[Times, Omega[k, l], ConjugateTranspose[Omega[k, l]]], {9, 9}], 
 {k, 0, 2}, {l, 0, 2}];

Do[
 Print[{k, l, FullSimplify[Sum[P[k, l][[i, i]], {i, 1, 9}]]}], 
 {k, 0, 2}, {l, 0, 2}]

by which I'm attempting to implement the construction (quoting from the beginning of sec. IV of https://arxiv.org/abs/0705.1403):

We start with a maximally entangled pure state, this is a Bell type state, in a chosen basis $\{0,1,2\}$ \begin{equation*} \Omega _{0 0} = \frac{1}{\sqrt{3}} \sum_{s=0}^2 |s\rangle \otimes |s\rangle\,. \end{equation*} On the first subspace, the system of Alice, we act with the Weyl operators, defined by $W_{k,l}|s\rangle = w^{k(s-l)}|s-l\rangle$ with $w=e^{2\pi i/3}$, while Bob's subsystem is always left inert. The indexes $k$ and $l$ run from $0$ to $2$. The other eight Bell states are constructed by acting with the Weyl operators onto the chosen Bell state \begin{equation*} \Omega_{k,l}=W_{k,l}\otimes I_3\; \Omega_{0,0}\,. \end{equation*} With that we can construct nine Bell projectors $P_{k,l}=|\Omega_{k, l}\rangle\langle\Omega_{k,l}|$. The mixtures of these pure states form our object of interest, the magic simplex $W$: \begin{equation*} W\;=\; \{ \quad \sum c_{k l}\;P_{k,l} \; |\; c_{kl}\geq 0 , \quad \sum c_{k l}=1 \quad \}\end{equation*}

However, I'm not presently sure the given program is fully successful in the intended implementation. The trace of $P_{k,l}$ is 1 for $k=0,l=0$ and $\frac{2}{9}$ in the other eight cases (as indicated in the output above), while for a pure state I would expect them to be all 1. Also, I'm not sure as to the use of the ConjugateTranspose in the 2nd Do expression. Further, I also looked into the use of TensorProduct, rather than KroneckerProduct. Also, is the ArrayReshape command in the 2nd Do expression really necessary?

Obviously, this all strongly pertains to the use of physics (Dirac) notation — which I note has been the subject of a number of questions on this site.

Answer 1

I'll try to express this problem along the lines of what I teach in my book Using Mathematica for Quantum Mechanics: A Student's Manual.

First, the basis states can be defined in many ways, for example as you do

S[0] = {1, 0, 0}; S[1] = {0, 1, 0}; S[2] = {0, 0, 1};

When taking tensor products of these, it's important to flatten out the result so that you end up with 9-vectors instead of 3x3-matrices:

Ω[0, 0] = 1/Sqrt[3] * Sum[Flatten[KroneckerProduct[S[s], S[s]]], {s, 0, 2}]
(*    {1/Sqrt[3], 0, 0, 0, 1/Sqrt[3], 0, 0, 0, 1/Sqrt[3]}    *)

This means that the result contains only the first, fifth, and ninth basis function of the combined 9-dimensional vector space. These basis vectors are given by

Tuples[{0, 1, 2}, 2]
(*    {{0,0}, {0,1}, {0,2}, {1,0}, {1,1}, {1,2}, {2,0}, {2,1}, {2,2}}    *)

so we see that $\Omega_{0,0}$ indeed only contains the basis functions $\lvert0,0\rangle=\lvert0\rangle\otimes\lvert0\rangle$, $\lvert1,1\rangle=\lvert1\rangle\otimes\lvert1\rangle$, and $\lvert2,2\rangle=\lvert2\rangle\otimes\lvert2\rangle$, as desired.

Now we define the Weyl operators with

w = E^(2*π*I/3);
W[k_, l_] := SparseArray[Table[{Mod[s-l, 3], s} + 1 -> w^(k*(s-l)), {s, 0, 2}], {3, 3}]

(using the mod function as suggested by @egwenesedai in the comments). They are 3x3-matrices. From here we get the Bell states

Ω[k_, l_] := KroneckerProduct[W[k, l], IdentityMatrix[3]] . Ω[0, 0]

where each $\Omega_{k,l}$ is a 9-vector, found by multiplying $\Omega_{0,0}$ by a 9x9-matrix, which in turn is the Kronecker product of two 3x3-matrices.

Finally, the Bell projectors are 9x9-matrices,

P[k_, l_] := KroneckerProduct[Ω[k, l], Conjugate[Ω[k, l]]]

and you can check that their traces are all 1:

Table[Tr[P[k, l]], {k, 0, 2}, {l, 0, 2}]
(*    {{1, 1, 1}, {1, 1, 1}, {1, 1, 1}}    *)

For deeper info I'd recommend you have a look at chapter 2 of my book referenced above.

Answer 2

Consulting the preceding paper, "The state space for two qutrits has a phase space structure in its core", of the three authors (B. Baumgartner, B. C. Hiesmayr and H. Narnhofer) https://arxiv.org/abs/quant-ph/0606083, we find in their sec. 3, the explicit representations of the nine Weyl operators $W_{k,l}$, (k, l =0, 1, 2):

[W[0, 0] = IdentityMatrix[3]; W[0, 1] = {{0, 1, 0}, {0, 0, 1}, {1, 0, 0}}; W[0, 2] = {{0, 0, 1}, {1, 0, 0}, {0, 1, 0}}; W[1, 0] = {{1, 0, 0}, {0, w, 0}, {0, 0,Conjugate[w]}};  W[1, 1] = {{0, 1, 0}, {0, 0, w}, {Conjugate[w], 0, 0}}; W[1, 2] = {{0, 0, 1}, {w, 0, 0},{0, Conjugate[w], 0}}; W[2, 0] = {{1, 0, 0}, {0, Conjugate[w], 0}, {0, 0, w}}; W[2, 1] = {{0, 1, 0}, {0, 0, Conjugate[w]}, {w, 0, 0}}; W[2, 2] = {{0, 0, 1}, {Conjugate[w], 0, 0}, {0, w, 0}}]

where $w = e^{\frac{2 i \pi }{3}}$.

Together, with the previously-given assertion, \begin{equation*} \Omega_{k,l}=W_{k,l}\otimes I_3\; \Omega_{0,0}\,, \end{equation*} it appears that the following code achieves the stated goal of constructing the "magical simplex $\mathcal{W}$ of bipartite qutrits"

S[0] = {1, 0, 0}; S[1] = {0, 1, 0}; S[2] = {0, 0, 1}; Omega[0, 0] = (1/Sqrt[3]) Sum[TensorProduct[S[s], S[s]], {s, 0, 2}]; Do[Omega1[k, l] = ArrayReshape[TensorProduct[W[k, l], IdentityMatrix[3] Omega[0, 0]], {9, 9}]/Sqrt[3], {k, 0, 2}, {l, 0, 2}]; Do[P[k, l] = Outer[Times, Omega1[k, l].ConjugateTranspose[Omega1[k, l]]], {k, 0, 2}, {l, 0, 2}]; Do[Print[{k, l, Sum[P[k, l][[i, i]], {i, 1, 9}]}], {k, 0, 2}, {l, 0, 2}]

The last command Do[Print[{k, l, Sum[P[k, l][[i, i]], {i, 1, 9}]}], {k, 0, 2}, {l, 0, 2}] confirms that the traces (sums of diagonal entries) of the nine Bell projectors $P_{k,l}$ are all unity.

(In light of the earlier comments of egwene sedai and me, let us also note the authors' statement that "Throughout this paper the letters $\{s,t,j,k,\ell,m,n,p,q\}$ denote the numbers $0, 1, 2$. Calculations with them are to be understood as ``modulo 3''.")

With this construction of the 8-dimensional simplex, we hope to obtain the Hilbert-Schmidt probability that a member (two-qutrit [$9 \times 9$ density matrix]) of the simplex, \begin{equation*} W\;=\; \{ \quad \sum c_{k l}\;P_{k,l} \; |\; c_{kl}\geq 0 , \quad \sum c_{k l}=1 \quad \}\end{equation*} has a positive partial transpose (PPT) and, further, the (sub)-probability that any member with a PPT is also separable--that is, the convex combination of product (qutrit) states. (Contrastingly, by the famous Peres-Horodecki result, for lower-dimensional $6 \times 6$ and $4 \times 4$ density matrices, the PPT-property is equivalent--not just necessary--for separability to hold. PPT-states that are not separable, are termed "bound entangled".)