Jaynes-Cummings model

In this tutorial we demonstrate how to get a CQED effective Hamiltonian using Pymablock with bosonic operators. As an example, we use the Jaynes-Cummings model, which describes a spin coupled to a boson. This tutorial shows how to use Pymablock with second-quantized operators, without the need to transform the Hamiltonian to a matrix representation.

Let’s start by importing the sympy functions we need to define the Hamiltonian. We will make use of sympy’s quantum mechanics module and its matrices.

from sympy import Matrix, Symbol, symbols, Eq, simplify
from sympy.physics.quantum.boson import BosonOp
from sympy.physics.quantum import Dagger

Define a second quantization Hamiltonian

We define the onsite energy \(\omega_r\), the energy gap \(\omega_q\), the perturbative parameter \(g\), and \(a\), the bosonic annihilation operator.

# resonator frequency, qubit frequency, Rabi coupling
wr, wq, g = symbols(r'\omega_r \omega_q g', real=True)

# resonator photon annihilation operator
a = BosonOp("a")

The Hamiltonian reads

H_0 = Matrix([[wr * Dagger(a) * a + wq / 2, 0], [0, wr * Dagger(a) * a - wq / 2]])
H_p = Matrix([[0,  g * a], [g * Dagger(a), 0]])

Eq(Symbol('H'), H_0 + H_p, evaluate=False)

\[\begin{split}\displaystyle H = \left[\begin{matrix}\frac{\omega_{q}}{2} + \omega_{r} {{a}^\dagger} {a} & g {a}\\g {{a}^\dagger} & - \frac{\omega_{q}}{2} + \omega_{r} {{a}^\dagger} {a}\end{matrix}\right]\end{split}\]

where the basis corresponds to the two spin states.

Get the Hamiltonian corrections

We can now define the block-diagonalization routine by calling block_diagonalize

%%time

import numpy as np
from pymablock import block_diagonalize

H_tilde, U, U_adjoint = block_diagonalize([H_0, H_p], symbols=[g])

CPU times: user 35.3 ms, sys: 1.08 ms, total: 36.4 ms
Wall time: 36.6 ms

The function block_diagonalize takes the Hamiltonian and the perturbative parameter as input. Differently from the rest of the tutorials, here we do not provide susbpace_vectors or subspace_indices. Pymablock treats the Hamiltonian as a single block, where the goal is to remove all terms that are not diagonal. The output therefore is a \(2 \times 2\) diagonal Hamiltonian that only contains one block with number operators.

Note

Pymablock only supports diagonal unperturbed Hamiltonians when using bosonic operators. This means that \(H_0\) must be block-diagonal and its entries need to be convertible to functions of the number operator, without single boson terms. This limitation may be lifted using advanced functionality, by providing a custom solve_sylvester input to the block_diagonalize function.

For example, to compute the 2nd order correction of the Hamiltonian of the \(↑, ↓\) subspaces we use

%%time

Eq(Symbol(r'\tilde{H}_2'), simplify(H_tilde[0, 0, 2]), evaluate=False)

CPU times: user 36.5 ms, sys: 123 μs, total: 36.6 ms
Wall time: 36.3 ms

\[\begin{split}\displaystyle \tilde{H}_2 = \left[\begin{matrix}\frac{g^{2} \left(1 + {N_{a}}\right)}{\omega_{q} - \omega_{r}} & 0\\0 & - \frac{g^{2} {N_{a}}}{\omega_{q} - \omega_{r}}\end{matrix}\right]\end{split}\]

Pymablock’s number operator

Pymablock’s output contains \(N_a = a^\dagger a\), the number operator for the bosonic mode, which is a NumberOperator object. Furthermore, the output is stored in the NumberOrderedForm class that Pymablock uses for efficient manipulation of second quantized expressions. In the example above, the diagonal entries of the Hamiltonian are NumberOrderedForm objects.

To see the result in terms of individual bosonic operators, we may use the doit method:

simplify(H_tilde[0, 0, 2])[0, 0].doit()

Check out the documentation for more information on how to use number operators and simplify expressions that contain them.

Higher order corrections to the Hamiltonian work exactly the same:

%%time

Eq(Symbol(r'\tilde{H}_4'), simplify(H_tilde[0, 0, 4]), evaluate=False)

CPU times: user 105 ms, sys: 1.8 ms, total: 107 ms
Wall time: 107 ms

\[\begin{split}\displaystyle \tilde{H}_4 = \left[\begin{matrix}- \frac{g^{4} \left(1 + {N_{a}}\right)^{2}}{\left(\omega_{q} - \omega_{r}\right)^{3}} & 0\\0 & \frac{g^{4} {N_{a}}^{2}}{\left(\omega_{q} - \omega_{r}\right)^{3}}\end{matrix}\right]\end{split}\]

%%time

Eq(Symbol(r'\tilde{H}_6'), simplify(H_tilde[0, 0, 6]), evaluate=False)

CPU times: user 1.59 s, sys: 30 ms, total: 1.62 s
Wall time: 1.62 s

\[\begin{split}\displaystyle \tilde{H}_6 = \left[\begin{matrix}\frac{2 g^{6} \left(1 + 3 {N_{a}} + 3 {N_{a}}^{2} + {N_{a}}^{3}\right)}{\omega_{q}^{5} - 5 \omega_{q}^{4} \omega_{r} + 10 \omega_{q}^{3} \omega_{r}^{2} - 10 \omega_{q}^{2} \omega_{r}^{3} + 5 \omega_{q} \omega_{r}^{4} - \omega_{r}^{5}} & 0\\0 & - \frac{2 g^{6} {N_{a}}^{3}}{\omega_{q}^{5} - 5 \omega_{q}^{4} \omega_{r} + 10 \omega_{q}^{3} \omega_{r}^{2} - 10 \omega_{q}^{2} \omega_{r}^{3} + 5 \omega_{q} \omega_{r}^{4} - \omega_{r}^{5}}\end{matrix}\right]\end{split}\]

We see that also computing the 6th order correction takes effectively no time.