Manipulating a Haar random unitary matrix and computing probabilities

Question

I am new to Mathematica and its many features. Here is a problem that I was trying to model.

Consider a $4 \times 4$ Haar random unitary $U$.

Define by $|00\rangle \langle 00|$ the following matrix:

\begin{equation}|00\rangle \langle 00| = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 &0 &0 \\ 0 & 0 &0 &0 \end{bmatrix}.\end{equation}

Additionally, define \begin{equation}X = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix},\end{equation} and define a two-qubit operator $\mathsf{X}_2 = X \otimes I,$ where $\otimes$ indicates tensor product and \begin{equation}I = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix},\end{equation} is the $2 \times 2$ identity matrix. So, explicitly,

\begin{equation}\mathsf{X}_2 = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 &0 &0 \\ 0 & 1 &0 &0 \end{bmatrix}.\end{equation}

---

Now, define an operator $\Phi$ as follows:

\begin{equation}\Phi(\color{red}{M}) = \frac{1}{2} \left[(U X_2 U) \color{red}{M} (U^* X_2 U^*) \right] + \frac{1}{2} \left[(U^2) \color{red}{M} ({U^*}^2) \right]\end{equation}

Let us consider the following final matrix:

\begin{equation}F = \Phi\big(\Phi\big(\cdots\Phi\big(|00\rangle \langle 00\big)\cdots\big)\big),\end{equation} where $\Phi$ is applied $k$ times.

---

Now, define a probability distribution $\mathsf{p}$ as follows. For every $x \in \{00, 01, 10, 11 \}$,

\begin{equation}\mathsf{p}(x) = \text{Tr}(E_x F).\end{equation}

where \begin{equation}E_{00} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 &0 &0 \\ 0 & 0 &0 &0 \end{bmatrix}, ~E_{01} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 &0 &0 \\ 0 & 0 &0 &0 \end{bmatrix}, \\~E_{10} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1 &0 \\ 0 & 0 &0 &0 \end{bmatrix},~\\E_{11} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 &0 \\ 0 & 0 &0 &1 \end{bmatrix}.\end{equation}

---

For example, when $k = 1$,

\begin{equation}\mathsf{p}(x) = \frac{1}{2}\text{Tr}\bigg[E_x (U~\mathsf{X}_2~U)|00\rangle \langle 00|~(U^{*} ~\mathsf{X}_2 ~U^{*})\bigg] + \frac{1}{2}\text{Tr}\bigg[E_x U^2~|00\rangle \langle 00|~\left(U^2\right)^{*}\bigg],\end{equation}

---

I want to arrange the elements of $p$ in descending order, and then plot them. The $y$ axis will be the sorted elements, the $x$ axis will be their indices (element $1$ has index $1$ and so on.)

I also want to write a code to do this for a general $k$, for $2^{n} \times 2^{n}$ matrices, and probability distributions with $n$ elements, where $\mathsf{X}_n = X \otimes I_{n-1}$, where $I_{n-1}$ is the $2^{n-1} \times 2^{n-1}$ identity matrix, and $E_x$ is extended naturally. I want help on how to automate the code instead of just brute force computing these values, for a particular $k$, for every possible $n$ and every possible $x$ and then plotting the result by manually creating a list.

---

For $k=1$, here is a brief sketch of what I tried.

I generated a unitary $U$ using the following command:

InitialMatrix = {{1, 0, 0, 0},
    {0, 0, 0, 0},
    {0, 0, 0, 0},
    {0, 0, 0, 0}};

U = RandomVariate[CircularUnitaryMatrixDistribution[4]];

X2 = {{0, 0, 1, 0},
    {0, 0, 0, 1},
    {1, 0, 0, 0},
    {0, 1, 0, 0}};

E00 =  {{1, 0, 0, 0},
    {0, 0, 0, 0},
    {0, 0, 0, 0},
    {0, 0, 0, 0}};

E01 = {{0, 0, 0, 0},
    {0, 1, 0, 0},
    {0, 0, 0, 0},
    {0, 0, 0, 0}};

E10 = {{0, 0, 0, 0},
    {0, 0, 0, 0},
    {0, 0, 1, 0},
    {0, 0, 0, 0}};

E11 = {{0, 0, 0, 0},
    {0, 0, 0, 0},
    {0, 0, 0, 0},
    {0, 0, 0, 1}};

p0 = Abs[1/2 Tr[E00.(U.X2.U).InitialMatrix.(ConjugateTranspose[U].X2.ConjugateTranspose[U])] + 1/2 Tr[E00.(U.U).InitialMatrix.(ConjugateTranspose[U].ConjugateTranspose[U])]];

p1 = Abs[1/2 Tr[E01.(U.X2.U).InitialMatrix.(ConjugateTranspose[U].X2.ConjugateTranspose[U])] + 1/2 Tr[E01.(U.U).InitialMatrix.(ConjugateTranspose[U].ConjugateTranspose[U])]];

p2 = Abs[1/2 Tr[E10.(U.X2.U).InitialMatrix.(ConjugateTranspose[U].X2.ConjugateTranspose[U])] + 1/2 Tr[E10.(U.U).InitialMatrix.(ConjugateTranspose[U].ConjugateTranspose[U])]];

p3 = Abs[1/2 Tr[E11.(U.X2.U).InitialMatrix.(ConjugateTranspose[U].X2.ConjugateTranspose[U])] + 1/2 Tr[E11.(U.U).InitialMatrix.(ConjugateTranspose[U].ConjugateTranspose[U])]];

Then, I sort $\{p_0, p_1, p_2, p_3\}$ in descending order: I call the new list $\{q_0, q_1, q_2, q_3\}$. Finally I plot $i$ in the horizontal axis and $q_i$ in the vertical axis.

What I want is a faster and more efficient way to do this and to generalize it to higher dimensions beyond $4$. I also want to generalize to $k$ more than $1$.

---

EDIT: I tried to generalize @user293787's answer below to $k > 1$. For $k=2$, the code can be brute-force modified as:

calculatep[U_] := With[{x = U[[;; , 1]]}, 1/4*(Abs[U.(U.x[[fliplist]])[[fliplist]]]^2 + Abs[U.U.x[[fliplist]]]^2 + Abs[U.(U.x)[[fliplist]]]^2 + Abs[U.U.U.x]^2)];

Is there a more elegant way to do this?

Answer 1

Code to calculate the probabilities when $U$ is a $2^n \times 2^n$ unitary matrix:

n=3;
fliplist=Join[Range[2^(n-1)+1,2^n],Range[1,2^(n-1)]];
calculatep[U_]:=With[{x=U[[;;,1]]},
                  1/2*(Abs[U.x[[fliplist]]]^2+Abs[U.x]^2)];

Example:

SeedRandom[1];
p = RandomVariate[CircularUnitaryMatrixDistribution[2^n]] // calculatep
(* {0.146664,0.0553558,0.173761,0.178616,0.0484086,0.084664,0.148739,0.163792} *)

To plot use ListPlot[ReverseSort[p]].