Basic partial trace and matrix multiplication questions

[Andrew Valentini]

Hello, I'm new to QuTiP and do not fully understand how to use its partial trace function or why I'm seemingly not able to multiply matrices of the same dimension. Here is what I am trying to do:

state0 = qt.basis(2, 0)
state1 = qt.basis(2, 1)

density_matrixA_0 = .25 * state0 * state0.dag()
density_matrixA_1 = .75 * state1 * state1.dag()

density_matrixA = density_matrix0 + density_matrix1

density_matrixB_0 = .75 * state0 * state0.dag()
density_matrixB_1 = .25 * state1 * state1.dag()

density_matrixB = density_matrix0 + density_matrix1

# I would then convert density_matrixA and B into numpy arrays here so that the multiplication below can be carried out

unitary = np.array([[1,0,0,0],[0,.5,.5,0],[0,-.5,.5,0],[0,0,0,1]])
unitary_dag = unitary.conj().T

entangled_12 = unitary*(np.kron(density_matrixA, density_matrixB))*unitary_dag

entangled_12 = qt.Qobj(entangled_12)
qubit1_trace = entangled_12.ptrace(0)
qubit2_trace = entangled_12.ptrace(1)

I am allowed to take entangled_12.ptrace(0) but the entangled_12.ptrace(1) line throws the error "Invalid selection index in ptrace" and I don't understand why. I'm guessing that I just don't understand what ptrace() is doing enough but I can't find any documentation online besides a single example of what I tried with entangled_12.ptrace(0). 

I also don't understand why I'm not able to compute:

 unitary*qt.tensor(density_matrixA, density_matrixB)*unitary_dag

when I instead initialize unitary as a QuTiP object. I get the error "incompatible dimensions [[4], [4]] and [[2, 2], [2, 2]]" when I try this. I'm assuming that this is because Python isn't actually taking the tensor product of density_matrixA and density_matrixB before multiplying unitary by it but, again, I'm not sure and there isn't good documentation on this online that I've found. 

I'd be very grateful if anyone was able to explain what I'm doing wrong here. I haven't been able to find help from documentation online. 

[Paul Menczel]

Hello,

The problem here is that qutip keeps track of the "dimensions" of quantum objects. In these dimensions, qutip remembers whether the Hilbert spaces are product spaces or not. For example, the dimensions of a 2-qubit space are `[2, 2]` and the dimension of a single 4-level system are `[4]` (even though both Hilbert spaces are 4-dimensional). An operator on the 2-qubit space (which is a map from the 2-qubit space to the 2-qubit space) therefore has dimensions `[[2, 2], [2, 2]]` and an operator on a 4-level system space has dimensions `[[4], [4]]`.

When you multiply two operators, qutip will check that the dimensions match. Because if they don't, there is a good chance that the user is doing something wrong.

When you want to take a partial trace, the dimensions must be set correctly so that qutip knows about the structure of the Hilbert space and can figure out where to take the trace.

In your case, the dimensions of `unitary` are probably not set correctly. When you create the Qobj from the numpy array, you must do `Qobj(np.array(...), dims=[[2, 2], [2, 2]])` explicitly.

This might seem like annoying extra effort, but it is also a good reminder to check that your 4x4 matrix representation uses the same basis that qutip uses. Remember that there is no canonical ordering of the basis of a product space. That is, you should check whether qutip uses the basis |-->, |-+>, |+->, |++> or the basis   |-->, |+->, |-+>, |++>, and which of those options your 4x4 matrix corresponds to. (Without checking, I do not know which one it is.)

Best

Paul