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?