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.
