How to generate all permutation matrices for 4 qubits?

Question

In quantum mechanics, a qubit can be understood as a 2 by 1 vector denoted by Dirac notation as $$|0\rangle \equiv \left( \begin{array}{c} 1\\ 0\\ \end{array} \right) ,|1\rangle \equiv \left( \begin{array}{c} 0\\ 1\\ \end{array} \right) .$$ The above notation stands for one qubit. For more qubits, we simply tensor them. E.g. there are four possible two qubits $$\left( \begin{array}{c} 1\\ 0\\ 0\\ 0\\ \end{array} \right) =|0\rangle \otimes |0\rangle \equiv |00\rangle ,|01\rangle ,|10\rangle ,|11\rangle .$$ I want to consider the permutation between different qubits, that is we permute the first qubit with the second qubit will lead $|01\rangle $ into $|10\rangle $ and $|10\rangle $ into $|01\rangle $ and keep $|00\rangle ,|11\rangle $ unchanged. The matrix stands for the above operation can be represented as $$|00\rangle \langle 00|+|01\rangle \langle 10|+|10\rangle \langle 01|+|11\rangle \langle 11|=\left( \begin{matrix} 1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\\ \end{matrix} \right) \tag{1}$$ where $\langle 10|$ stands for transpose of $|10\rangle $(cases we are considering can be limited to real numbers). There are only two permutation matrices for two qubits. They are identity matrix and eq(1) corresponding to do not change and change the first qubit with the second qubit respectively.

I want to know, how can I use Mathematica to generate all the permutation matrices (16 by 16) similar to eq(1) to 4 qubits? There are $4!=24$ permutation matrices for 4 qubits.

Answer 1

Assuming I understand what you're asking, you want a Mathematica implementation of the operator $P_\pi$ with action $$P_\pi = \sum_{i_1,..., i_n} |\pi(i_1)\cdots \pi(i_n)\rangle\!\langle i_1,...,i_n|,$$ for arbitrary permutations $\pi\in S_n$. The only minor issue to do this is that if you want this to act as a standard matrix, you'll have to convert indices of input vectors into lists of indices corresponding to each qubit. You can do this for example doing a base-2 decomposition (assuming we're talking qubits; otherwise just use the suitable basis). Here's a possible implementation:

indexToQubitIndices[n_, numQubits_] := PadLeft[IntegerDigits[n - 1, 2], numQubits];
permutationMatrix[permutation_] := With[{numQubits = Length @ permutation},
    SparseArray[
        {
            {i_, j_} :> 1 /; Equal[
                Permute[indexToQubitIndices[i, numQubits], permutation],
                indexToQubitIndices[j, numQubits]
            ]
        },
        2^numQubits {1,1}
    ]
];
permutationMatrix[{1, 3, 2}] // Normal // MatrixForm

Note how you retrieve identity and Swap with permutationMatrix@{1, 2} and permutationMatrix@{2, 1}, respectively.

If you now want all permutations, just use Permutations@{1, 2, 3, 4} to generate all permutations and use them as the second element for permutationMatrix above. Eg something like permutationMatrix /@ Permutations @ Range[4].

Answer 2

Looks like in the bra ket tensor notation you search for something like

mm[16] = Normal /@ PermutationMatrix /@ Permutations[{1, 2, 3, 4}]

TableForm@Partition[mm[16],4]