I would like to study the effects of non-computational levels of two coupled qubits by a resonator using QuTiP. The idea is to consider the following Hamiltonians
Qubits Hamiltonian $$ H_{q_i} = \sum_n 4E_C (n-n_g)^2|n\rangle\langle n| + \frac{E_J}{2}(|n+1\rangle\langle n|+|n\rangle\langle n+1|) $$
Resonator Hamiltonian $$ H_r = \omega_r a^\dagger a $$
Free Hamiltonian $$ H_0 = H_{q_1} \otimes H_{q_2} \otimes H_r $$
Interacting Hamiltonian $$ H_{1r} = Q_1\otimes \mathbb{I} \otimes Q_r/C = -2i\beta_1 e V_\text{rms} N_{q_1} \otimes \mathbb{I} \otimes (a - a^\dagger) $$ $$ H_{2r} = \mathbb{I} \otimes Q_2 \otimes Q_r/C = \mathbb{I} \otimes -2i\beta_2 e V_\text{rms} N_{q_2} \otimes (a - a^\dagger) $$ $$ H_\text{int} = H_{1r} + H_{2r} $$
where $N_{q_i}$ is the number operator in the energy basis, $\beta_i = \frac{C_{g_i}}{C_{g_i} + C_{s_i}}$ which is proportional to the coupling capacitance between the qubit and the resonator, and the shunt capacitance.
Total Hamiltonian $$ H = H_0 + H_\text{int} $$
After instantiating these Hamiltonians I want to use the Schrieffer-Wolff Transformation to get an effective Hamiltonian limited to the computational levels and decoupled from the resonator.
The basic idea is to implement the Schrieffer-Wolff transformation to reduce the transmon Hamiltonian (multiple levels) to an effective Hamiltonian (two levels).
$$ H = H_0 + \varepsilon V$$
introducing the unitary transformation
$$ e^{\varepsilon S} H e^{-\varepsilon S}$$
where $\varepsilon$ is small, $\exp{(\varepsilon S)}$ is unitary and $S$ anti-hermitian. By expanding the expression using the Baker-Campbell-Hausdorff formula
$$ e^{\varepsilon S} H e^{-\varepsilon S} = H + \varepsilon[S, H] + \frac{\varepsilon^2}{2}[S, [S, H]] + O(\varepsilon^2) = H_0 + \varepsilon(V + [S, H_0]) + \frac{\varepsilon^2}{2}(2[S,V] + [S,[S,H_0]]) + O(\varepsilon^2)$$
The Schrieffer-Wolff transformation chooses $S$ such that
$$V + [S, H_0] = 0$$
which is known as Sylvester equation. What remains is
$$e^{\varepsilon S} H e^{-\varepsilon S} = H_0 + \frac{\varepsilon^2}{2}[S,V] + O(\varepsilon^2)$$
By defining the projectors $P$ on the low-energy space. The low-energy effective Hamiltonian becomes
$$H_\text{eff} = P(H_0 + \frac{\varepsilon^2}{2}[S,V] )P$$
# Eigenvectors
_, eigenvectors_q1 = q1.h.eigenstates()
_, eigenvectors_q2 = q2.h.eigenstates()
_, eigenvectors_r = r.h.eigenstates()
# Projectors (|0><0|+|1><1|) ⊗ (|0><0|+|1><1|) ⊗ Identity
# Each tensor product is an energy eigenstate of qubits' Hamiltonian, for the resonator an identity is used
# Qubit 1 ⊗ Qubit 2 ⊗ Resonator
p = tensor(eigenvectors_q1[0] * eigenvectors_q1[0].dag() + eigenvectors_q1[1] * eigenvectors_q1[1].dag(), eigenvectors_q2[0] * eigenvectors_q2[0].dag() + eigenvectors_q2[1] * eigenvectors_q2[1].dag(), qeye(r.dim))
# S matrix from Sylvester equation (https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_sylvester.html)
s = Qobj(solve_sylvester(-h0.full(), h0.full(), -hint.full()))
# Effective Hamiltonian
heff = Qobj(h0.full() + 0.5 * commutator(s.full(), hint.full()))
# Reduced Hamiltonian
hred = Qobj(p * heff.full() * p)
When I use hred as the Hamiltonian for the Master Equation and I use psi100 as the initial state
psi100 = tensor(eigenvectors_q1[1], eigenvectors_q2[0], eigenvectors_r[0])
psi010 = tensor(eigenvectors_q1[1], eigenvectors_q2[0], eigenvectors_r[0])
t_list = np.linspace(0, 20, 1000)
result = mesolve(H=hred, rho0=psi[1,0,0], tlist=t_list, c_ops=[], e_ops=[], args={})
There is no exchange between states psi100 and psi010, when I would expect, in an ideal two-level system an iSWAP implementation. Or at least some kind of interaction.
As for the Hamiltonian definitions, I am quite sure that they are okay. Where I have some doubts is about the Schrieffer-Wolff transformation, I might have wrong implemented some steps or I am missing some diagonalizations, wrong states, and wrong projector definition.
