Extract single qubit state from combined state in QuTiP using a Hamiltonian

Question

This questions is very similar to another question asked a couple of months ago. I would like to extract the state of single qubits from combined states, using QuTip (the Quantum Toolbox in Python) using a Hamiltonian simulator.

Example

For example, suppose I begin with two qubits, |0⟩ and |1⟩. I can create these in QuTiP (as q1 and q2) as follows:

q1 = basis(2,1)
q2 = basis(2,0)

$q_1 = \begin{bmatrix}0\\1\end{bmatrix} \quad q_2 = \begin{bmatrix}1\\0\end{bmatrix}$

I can then combine q1 and q2 using tensor multiplication:

q1q2 = tensor(q1, q2)

$q_1 \otimes q_2 = \begin{bmatrix}0\\0\\1\\0\end{bmatrix}$

Now, suppose that I apply a CNOT gate to $q_1\otimes q_2$. We can define CNOT as:

circ =  QubitCircuit(2)
circ.add_gate("CNOT",targets=[1],controls=[0])
U = gate_sequence_product(circ.propagators())

$U=\begin{bmatrix} 1&0&0&0 \\ 0&1&0&0 \\ 0&0&0&1 \\ 0&0&1&0 \\ \end{bmatrix}$

Probability of $q_2$

I now want to figure what q2 would look like after application of this gate. We can do this using the partial trace of the density matrix:

M = ket2dm(U*q1q2).ptrace(0)
q_2' = np.diag(np.sqrt(M.full()))

$M = \begin{bmatrix}0&0\\0&1 \end{bmatrix}\\ q_2' = \begin{bmatrix}0\\1 \end{bmatrix}$

Probability of $q_2$ version 2

I would like to do the same now, but using a Hamiltonian simulator. For this, I use the built-in master equation solver which solves:

$iℏ\frac{d}{d}t|ψ⟩ = H |ψ⟩$

The code for this is:

times = np.linspace(0.0, 20.0, 100)
data = mesolve(U, q1q2, times)
M2 = ket2dm(data.states[-1]).ptrace(0)
q_2'' = np.diag(np.sqrt(M2.full()))

$M_2 = \begin{bmatrix}0.167&0.373j\\-0.373j&0.833 \end{bmatrix}\\ q_2'' = \begin{bmatrix}0.40809641\\0.91293884 \end{bmatrix}$

Which is not what I expected. Question: can I simply use my unitary matrix here or do I need to do some transformation before/after running the simulator? Are there other ways to come up with the probabilities of one qubit, besides using ptrace?

Answer 1

Say $U$ is a unitary evolution (the CNOT in your case), and $\vert\psi\rangle$ an initial state.

Evolving $|\psi\rangle$ with $U$ amounts to computing $U\lvert\psi\rangle$.

If on the other hand you want to study the dynamics of this evolution, then you need a generator for the evolution. In other words, you want an Hamiltonian $H$ such that after a certain time $t$ you get $$e^{-i t H}\lvert\psi\rangle=U\lvert\psi\rangle.$$ Such an Hamiltonian is what you need to feed to mesolve. Here is a working version of the code doing just this:

import numpy as np
import scipy
import qutip

cnot_generator = qutip.Qobj(1j * scipy.linalg.logm(qutip.cnot().full()), dims=[[2] * 2] * 2)
initial_state = qutip.tensor(qutip.basis(2, 1), qutip.basis(2, 0))
times = np.arange(0, 1.1, 0.1)
qutip.mesolve(cnot_generator, initial_state, times).states

Note how the last state given by the above is (in good approximation) the same you obtained yourself. I used ket notation here because, being everything unitary, there is no need to use density matrices. You may adapt the code appropriately to use density matrices.