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

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)
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?

• I don't know QuTip, but I'd be surprised if its master-equation solver took the unitary time-evolution operator as input. If you know $U$ then you have already solved the master equation. Typically you would only know the generator of time evolution (in this case $\propto \ln U$). Mar 2 '18 at 9:48
• Thank you @MarkMitchison! I assume this is because any Hamiltonian can be described as a unitary like so: $U=e^{−iHt}{\hbar}$, right? That would be swell, because then we could state the input as $H(t) = \frac{i \hbar}{t} \ln(U)$
– ciri
Mar 3 '18 at 8:00
• Yes, you need to input the Hamiltonian that generates the unitary dynamics if you want to solve the equations of motion. If you already know the unitary, then you simply apply it directly to the initial state to obtain the final state. Mar 4 '18 at 11:19

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.

• This is exactly what I was looking for, thank you @glS!
– ciri
Mar 6 '18 at 10:35