Simple and efficient construction of tensor product operators

Question

Problem Description

I'm using mathematica to study a problem in quantum mechanics, which is naturally understood in terms of vector spaces. In the simplest example, suppose we have a 16-dimensional vector space built out of the tensor product of four two-dimensional spaces. I need a simple (and ideally, efficient) code to construct operators which act on different subsystems. I think this is doable, but I'm struggling with correctly organizing arrays.

Here is an example. Suppose I have a random unitary matrix $U_{12}$ which acts on the first two two-dimensional spaces, but leaves the others alone. So we can write it as $U_{12}\otimes I_{4}$, where $I_{4}$ is the $4\times4$ identity matrix. This is easy to construct in mathematica:

U12 = RandomVariate[CircularUnitaryMatrixDistribution[4];
UBig = KroneckerProduct[U12, IdentityMatrix[4]];

This does what I want it to. What I'm in need of is a code that generates a matrix representation of the analogous operator $U_{14}$, which acts non-trivially on subsystems 1 and 4. More generally, it would be nice to have a code which can generate matrix representations of operators that act on subsystems $i,j$ when they aren't neighbors (i.e., KroneckerProduct would only work if $i = j\pm 1$). Also, my case of interest would involve roughly ten such two-dimensional spaces, so it would be nice if the solution scaled decently.

My attempt at a solution

Here is my thinking as I have tried to find a solution. Nothing has worked yet so it would be great if someone could tell me where I've gone wrong or suggest an alternative. We could start by generating a random unitary matrix that acts on subspaces 1 and 2, and reshaping the array to have the correct index structure:

U12 = RandomVariate[CircularUnitaryMatrixDistribution[4];
UBig = KroneckerProduct[U12, IdentityMatrix[4]];
UBig2 = ArrayReshape[UBig, {2,2,2,2,4,4}]

The thinking here is that I'll want to re-arrange indices such that we end up with U12 acting on subspaces (1,4). This can only happen if we have 2-valued indices to rearrange; the 4-valued indices should be blocks left alone to parameterize the identity for the subspace we aren't touching. Now let's rearrange the indices:

UBig3 = TensorTranspose[UBig2, Cycles[{{2, 5}, {4, 6}}]]

I think this is the correct cycle based on a diagram I've drawn - the dimensions of this tensor are {2,4,2,4,2,2}, where I understand the first and third indices index subspace 1, while indices 5 and 6 index subspace 4, and the 4-valued indices are the identity space (subspaces 2,3). Now I think we want to flatten these indices such that the there is a single "incoming" and "outgoing index for the (1,4) subspace, i.e.,

UBig4 = Flatten[UBig3, {{1,5},{3,6}};

Now we have a 16x16 matrix, which is indeed what I want. Sadly, as you can check, it is not unitary, so something has gone wrong. This operation should preserve unitarity.

Answer 1

This type of computation can be done in two different ways, using TensorProduct or using KroneckerProduct. I find easier to use TensorProduct, but I give here how I'd address this with both methods:

Define the following objects and note how similar the inputs are in the last two lines:

Id2 = IdentityMatrix[2];
U12 = RandomVariate[CircularUnitaryMatrixDistribution[4]];
U12T = Transpose[ArrayReshape[U12, {2, 2, 2, 2}], {1, 3, 2, 4}];
UBigK = KroneckerProduct[U12, Id2, Id2];
UBigT = TensorProduct[U12T, Id2, Id2];

UBigK and UBigT contain the same inofmation, and in fact we can convert each one into the other:

In[]:= Flatten[UBigT, {{1, 3, 5, 7}, {2, 4, 6, 8}}] === UBigK
Out[]= True

In[]:= Transpose[ArrayReshape[UBigK, {2, 2, 2, 2, 2, 2, 2, 2}], {1, 3, 5, 7, 2, 4, 6, 8}] === UBigT
Out[]= True

Let us first discuss UBigT. It is an operator with 8 indices, organized by four consecutive pairs corresponding to operation in your for vector spaces. If we want to exchange the roles of spaces 2 and 4 then we need to change indices 3 and 4 with indices 7 and 8 respectively. Then we flatten again as we did before to get the 16x16 unitary matrix you want:

In[]:= UnitaryMatrixQ[res1 = Flatten[Transpose[UBigT, Cycles[{{3, 7}, {4, 8}}]], {{1, 3, 5, 7}, {2, 4, 6, 8}}]]
Out[]= True

To perform the computation directly on UBigK seems more difficult to me because it requires a nontrivial permutation. A possible way to get it is as follows. Construct a representation of the basis states:

In[]:= states = Ket @@@ Tuples[{0, 1}, 4]
Out[]= {Ket[0, 0, 0, 0], Ket[0, 0, 0, 1], Ket[0, 0, 1, 0], Ket[0, 0, 1, 1],
        Ket[0, 1, 0, 0], Ket[0, 1, 0, 1], Ket[0, 1, 1, 0], Ket[0, 1, 1, 1],
        Ket[1, 0, 0, 0], Ket[1, 0, 0, 1], Ket[1, 0, 1, 0], Ket[1, 0, 1, 1],
        Ket[1, 1, 0, 0], Ket[1, 1, 0, 1], Ket[1, 1, 1, 0], Ket[1, 1, 1, 1]}

Now find the permutation that exchanges slots 2 and 4:

In[]:= perm24 = PermutationList[FindPermutation[states, states /. Ket[a1_, a2_, a3_, a4_] :> Ket[a1, a4, a3, a2]], 16]
Out[]= {1, 5, 3, 7, 2, 6, 4, 8, 9, 13, 11, 15, 10, 14, 12, 16}

Apply to both indices of the matrix:

In[]:= UnitaryMatrixQ[res2 = Part[UBigK, perm24, perm24]]
Out[]= True

Check that the two results coincide:

In[]:= res1 === res2
Out[]= True