Jaynes-Cummings ladder

The Jaynes-Cummings model was introduced in quantum optics to describe interaction of a two-level atom with a single quantized mode of an optical cavity [JC63]. Phonons in the mode are created/destroyed by bosonic ladder operators \(a^\dagger\)/\(a\). Atom’s Hilbert space is spanned by the ground state \(|g\rangle\) and the excited state \(|e\rangle\), and is equivalent to that of a spin-half. The operators acting in this isospin space are raising and lowering operators \(S_+ = |e\rangle\langle g|\), \(S_- = |g\rangle\langle e|\) and the atomic inversion operator \(2S_z = |e\rangle\langle e| - |g\rangle\langle g|\).

Under the assumpion that the photon frequency \(\omega_c\) and the atomic transition frequency \(\omega_a\) satisfy \(|\omega_c - \omega_a| \ll \omega_c + \omega_a\) (small detuning \(\delta = \omega_a - \omega_c\)), one can write down the Jaynes-Cummings Hamiltonian in the rotating wave approximation as follows,

\[\hat H = \hbar \omega_c a^\dagger a + \hbar \omega_a S_z + \hbar \frac{\Omega}{2}(a S_+ + a^\dagger S_-),\]

where \(\Omega\) is the atom-field interaction strength. The Jaynes-Cummings Hamiltonian commutes with operator \(\hat N = a^\dagger a + S_z\), which means \(\hat H\) takes on a block-digonal structure. The ground state \(|n = 0, g\rangle\) is within its own one-dimensional block, while the rest of the blocks are spanned by a pair \(|n-1, e\rangle\), \(|n, g\rangle\) each (\(n = \overline{1,\infty}\)). Accordingly, the ground state energy is \(E_g = -(1/2)\hbar\omega_a\), while all excited states form a series of doublets \(E_\pm(n)\) – the so called Jaynes-Cummings ladder.

\[E_\pm(n) = \hbar\omega_c \left(n - \frac{1}{2}\right) \pm \frac{\hbar}{2} \sqrt{\delta^2 + n\Omega^2}.\]

The respective eigenstates in each block (dressed states) are

\[\begin{split}|n, +\rangle &= \cos\left(\frac{\alpha_n}{2}\right) |n - 1, e\rangle + \sin\left(\frac{\alpha_n}{2}\right) |n, g\rangle,\\ |n, -\rangle &= \sin\left(\frac{\alpha_n}{2}\right) |n - 1, e\rangle - \cos\left(\frac{\alpha_n}{2}\right) |n, g\rangle\\\end{split}\]

with mixing angle

\[\alpha_n = \tan^{-1}\left(\frac{\Omega\sqrt{n}}{\delta}\right).\]

The program below verifies these analytical results numerically.

  1#include <algorithm>
  2#include <cmath>
  3#include <iostream>
  4#include <vector>
  5
  6#include <libcommute/libcommute.hpp>
  7
  8using namespace libcommute;
  9
 10int main() {
 11
 12  // Planck constant
 13  double const hbar = 1.0;
 14  // Angular frequency of the cavity mode
 15  double const w_c = 1.0;
 16  // Atomic transition frequency
 17  double const w_a = 1.1;
 18  // Atom-field coupling constant
 19  double const Omega = 0.2;
 20  // Detuning
 21  double const delta = w_a - w_c;
 22
 23  // For a_dag(), a(), S_p(), S_m() and S_z()
 24  using namespace static_indices;
 25
 26  // The Jaynes-Cummings Hamiltonian
 27  auto H = hbar * w_c * a_dag() * a() + hbar * w_a * S_z() +
 28           hbar * 0.5 * Omega * (a() * S_p() + a_dag() * S_m());
 29
 30  // Check that the number operator commutes with the Hamiltonian
 31  auto N = a_dag() * a() + S_z();
 32  std::cout << "[H, N] = " << (H * N - N * H) << '\n';
 33
 34  // Complete Hilbert space of the system is a product of the isospin-1/2 space
 35  // and the bosonic space truncated to 2^10 - 1 excitations.
 36  auto es_spin = static_indices::make_space_spin(0.5);
 37  auto es_boson = static_indices::make_space_boson(10);
 38
 39  hilbert_space</* Our operators have no indices */> hs(es_spin, es_boson);
 40
 41  // Ground state energy
 42  double const E_gs = -hbar * w_a / 2;
 43
 44  // Excited energy levels E_{\pm}(n)
 45  auto energy = [=](int pm, int n) {
 46    return hbar * w_c * (n - 0.5) +
 47           pm * 0.5 * hbar * std::sqrt(delta * delta + n * Omega * Omega);
 48  };
 49
 50  // Mixing angle \alpha_n
 51  auto alpha = [&](int n) { return std::atan(Omega * std::sqrt(n) / delta); };
 52
 53  // Two state vectors in the complete Hilbert space
 54  std::vector<double> phi(hs.dim()), psi(hs.dim());
 55
 56  // L^2 norm of |\spi> - E|\phi>
 57  // Will be used to verify correctness of the eigenpairs
 58  auto diff = [&](std::vector<double> const& psi,
 59                  std::vector<double> const& phi,
 60                  double E) {
 61    double d = 0;
 62    for(unsigned int i = 0; i < hs.dim(); ++i) {
 63      d += std::pow(psi[i] - E * phi[i], 2);
 64    }
 65    return std::sqrt(d);
 66  };
 67
 68  // loperator object corresponding to H
 69  auto Hop = make_loperator(H, hs);
 70
 71  //
 72  // Verify that |0, g> is the ground state
 73  //
 74
 75  // Index of the |0, g> basis vector
 76  auto gs_index = hs.basis_state_index(es_boson, 0) + // n = 0
 77                  hs.basis_state_index(es_spin, 0);   // 0 -> g
 78
 79  // Set |\phi> = |0, g>
 80  std::fill(phi.begin(), phi.end(), 0);
 81  phi[gs_index] = 1;
 82
 83  psi = Hop(phi);
 84
 85  // Check |\psi> \approx E_{gs} |\phi>
 86  std::cout << "|||0,g> - E_g|0,g>||_{L^2} = " << diff(psi, phi, E_gs) << '\n';
 87
 88  //
 89  // Now, do a similar check for the Jaynes-Cummings ladder doublets.
 90  //
 91
 92  // Check the first 20 doublets
 93  for(int n = 1; n <= 20; ++n) {
 94
 95    // Index of basis state |n - 1, e>
 96    auto n1_e_index = hs.basis_state_index(es_boson, n - 1) + // n - 1
 97                      hs.basis_state_index(es_spin, 1);       // 1 -> e
 98    // Index of basis state |n, g>
 99    auto n_g_index = hs.basis_state_index(es_boson, n) + // n
100                     hs.basis_state_index(es_spin, 0);   // 0 -> g
101
102    // Set \phi to be the "+" dressed state
103    std::fill(phi.begin(), phi.end(), 0);
104    phi[n1_e_index] = std::cos(alpha(n) / 2);
105    phi[n_g_index] = std::sin(alpha(n) / 2);
106
107    psi = Hop(phi);
108
109    std::cout << "|||" << n << ",+> - E_+(" << n << ")|" << n
110              << ",+>||_{L^2} = " << diff(psi, phi, energy(1, n)) << '\n';
111
112    // Set \phi to be the "-" dressed state
113    std::fill(phi.begin(), phi.end(), 0);
114    phi[n1_e_index] = std::sin(alpha(n) / 2);
115    phi[n_g_index] = -std::cos(alpha(n) / 2);
116
117    psi = Hop(phi);
118
119    std::cout << "|||" << n << ",-> - E_-(" << n << ")|" << n
120              << ",->||_{L^2} = " << diff(psi, phi, energy(-1, n)) << '\n';
121  }
122
123  return 0;
124}
[JC63]

“Comparison of quantum and semiclassical radiation theories with application to the beam maser”, E. T. Jaynes and F.W. Cummings, Proc. IEEE, Vol. 51, issue 1, pp. 89-109 (1963), https://doi.org/10.1109/PROC.1963.1664