Jaynes-Cummings model

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

def display_eq(title, expr):
    return Eq(Symbol(title), expr, evaluate=False)

display_eq('H', H_0 + H_p)
\[\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 31.6 ms, sys: 0 ns, total: 31.6 ms
Wall time: 31.8 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

display_eq(r'\tilde{H}_2', H_tilde[0, 0, 2])

CPU times: user 24 ms, sys: 0 ns, total: 24 ms
Wall time: 23.8 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

display_eq(r'\tilde{H}_4', H_tilde[0, 0, 4])

CPU times: user 105 ms, sys: 177 μs, total: 105 ms
Wall time: 105 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

display_eq(r'\tilde{H}_6', H_tilde[0, 0, 6])

CPU times: user 211 ms, sys: 215 μs, total: 211 ms
Wall time: 211 ms
\[\begin{split}\displaystyle \tilde{H}_6 = \left[\begin{matrix}\frac{2 g^{6} \left(1 + {N_{a}}\right)^{3}}{\left(\omega_{q} - \omega_{r}\right)^{5}} & 0\\0 & - \frac{2 g^{6} {N_{a}}^{3}}{\left(\omega_{q} - \omega_{r}\right)^{5}}\end{matrix}\right]\end{split}\]

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

Spin operators#

An alternative to defining the Hamiltonian as a 2x2 matrix is to use spin operators:

from sympy.physics.quantum import pauli

H_0 = wr * Dagger(a) * a + wq * pauli.SigmaZ("s") / 2
H_p = g * (pauli.SigmaPlus("s") * a + pauli.SigmaMinus("s") * Dagger(a))

display_eq('H', H_0 + H_p)
\[\displaystyle H = \frac{\omega_{q} {\sigma_z^{(s)}}}{2} + \omega_{r} {{a}^\dagger} {a} + g \left({\sigma_-^{(s)}} {{a}^\dagger} + {\sigma_+^{(s)}} {a}\right)\]

Similarly as before, we provide the Hamiltonian to block_diagonalize:

H_tilde, *_ = block_diagonalize([H_0, H_p], symbols=[g])

display_eq(r'\tilde{H}_4', simplify(H_tilde[0, 0, 4]))
\[\displaystyle \tilde{H}_4 = \frac{g^{4} \left({N_{a}}^{2} \left(1 - {N_{s}}\right) - {N_{s}} \left(1 + {N_{a}}\right)^{2}\right)}{\left(\omega_{q} - \omega_{r}\right)^{3}}\]

The fourth order correction to the Hamiltonian depends on the spin occupation number \(N_s\), such that for \(N_s = 0, 1\) we recover the diagonal entries of the previous example.

Contents