How do I implement $SU(2)$ rotation on the Bloch Sphere unsing qutip?

Question

Using qutip I am trying to implement a qubit rotation according the formula $(25)$ provided in this document "Lecture notes: Qubit representations and rotations":

$$ U_{\hat n}(\theta)=\sigma_0\cos\frac{\theta}{2}-i(\hat n\cdot \sigma)\sin\frac{\theta}{2}\\ $$

$$ M_{q'}=U_{\hat n}(\theta)\cdot M_q\cdot U_{\hat n}^A(\theta) $$

import qutip as qt
from qutip.qip.operations import rx
import numpy as np

def to_spherical(state):
    r0 = np.abs(state[0])
    ϕ0 = np.angle(state[0])
    r1 = np.abs(state[1])
    ϕ1 = np.angle(state[1])
    r = np.sqrt(r0 ** 2 + r1 ** 2)
    θ = 2 * np.arccos(r0 / r)
    ϕ = ϕ1 - ϕ0
    return [r, θ, ϕ]

def rn_su2(θ, state, nx, ny, nz):
    Ψ = [state.data[0,0], state.data[1,0]]
    arr = to_spherical(Ψ)
    s_θ = arr[1]
    s_ϕ = arr[2]
    M_q = np.sin(s_θ)*np.cos(s_ϕ)*qt.sigmax() + np.sin(s_θ)*np.sin(s_ϕ)*qt.sigmay() + np.cos(s_θ)*qt.sigmaz()
    U_n = qt.qeye(2)*np.cos(θ/2) -1j*(nx*qt.sigmax()+ny*qt.sigmay()+nz*qt.sigmaz())*np.sin(θ/2)
    r_state = U_n*M_q*U_n.dag()
    return r_state

b = qt.Bloch()
b.clear()
b.make_sphere()

states = []
alpha = 1/np.sqrt(2)
beta = 1/np.sqrt(2)
s = np.array([alpha,beta])
state = qt.Qobj(s)
states.append(state)

rotated = rn_su2(np.pi/8, state, 0, 0, 1)
states.append(rotated)

rotated = rn_su2(np.pi/8, rotated, 0, 0, 1)
states.append(rotated)

b.add_states(states)
b.show()

Unfortunatelly the rotated qubits lie not on the Bloch Sphere anymore:

I would appreciate any help in getting this rotation working. I uploaded the complete notebook publicly at GitHub.

Answer 1

It looks like you haven't converted a normal vector into a density matrix quite right. Recall that given a normal vector $\hat{n} = (n_x, n_y, n_z)$ you can write the corresponding density matrix on the Bloch sphere as

\begin{align} \rho = \frac{1}{2} \left(I + \sigma_x n_x + \sigma_y n_y + \sigma_z n_z \right) \tag{1} \end{align}

(see Exercise 2.72 in Nielsen and Chuang, for example). So you've correctly converted to spherical coordinates (e.g. $n_x = \sin \theta \cos \phi$) but have forgotten to add the factor of $1/2$ and identity component. Here's the fixed line:

M_q = np.eye(2) / 2 + 1/2 * (np.sin(s_θ)*np.cos(s_ϕ)*qt.sigmax() + np.sin(s_θ)*np.sin(s_ϕ)*qt.sigmay() + np.cos(s_θ)*qt.sigmaz())

For future reference, when you run into this kind of issue its good to check numerically that the matrices and arrays you're working with make physical sense. It should always be true that $\text{tr}(\rho)=1$ and $\rho = \rho^\dagger$ for any state $\rho$, $U^\dagger U = I$ for any unitary $U$, $\lVert \hat{n}\rVert^2 = 1$ for any normal vector $\hat{n}$, and so on. Checking the first condition immediately revealed the source of the bug.