Quantised Light-Atom Interactions in the Semi-classical Limit: Recovering the Periodically Driven Atom

This example shows how to code up the Hamiltonians:

\[\begin{split}H&=& \Omega a^\dagger a + \frac{A}{2}\frac{1}{\sqrt{N_\mathrm{ph}}}\left(a^\dagger + a\right)\sigma^x + \Delta\sigma^z, \nonumber\\ H_\mathrm{sc}(t) &=& A\cos\Omega t\;\sigma^x + \Delta\sigma^z.\end{split}\]

Details about the code below can be found in SciPost Phys. 2, 003 (2017).

Script

download script

  1########################################################################################
  2#                                   example 3                                          #
  3#    In this example we show how to use the photon_basis class to study spin chains    #
  4#    coupled to a single photon mode. To demonstrate this we simulate a single spin    #
  5#    and show how the semi-classical limit emerges in the limit that the number of     #
  6#    photons goes to infinity.                                                         #
  7########################################################################################
  8from quspin.basis import spin_basis_1d, photon_basis  # Hilbert space bases
  9from quspin.operators import hamiltonian  # Hamiltonian and observables
 10from quspin.tools.measurements import obs_vs_time  # t_dep measurements
 11from quspin.tools.Floquet import Floquet, Floquet_t_vec  # Floquet Hamiltonian
 12from quspin.basis.photon import coherent_state  # HO coherent state
 13import numpy as np  # generic math functions
 14
 15#
 16##### define model parameters #####
 17Nph_tot = 60  # maximum photon occupation
 18Nph = Nph_tot / 2  # mean number of photons in initial coherent state
 19Omega = 3.5  # drive frequency
 20A = 0.8  # spin-photon coupling strength (drive amplitude)
 21Delta = 1.0  # difference between atom energy levels
 22#
 23##### set up photon-atom Hamiltonian #####
 24# define operator site-coupling lists
 25ph_energy = [[Omega]]  # photon energy
 26at_energy = [[Delta, 0]]  # atom energy
 27absorb = [[A / (2.0 * np.sqrt(Nph)), 0]]  # absorption term
 28emit = [[A / (2.0 * np.sqrt(Nph)), 0]]  # emission term
 29# define static and dynamics lists
 30static = [["|n", ph_energy], ["x|-", absorb], ["x|+", emit], ["z|", at_energy]]
 31dynamic = []
 32# compute atom-photon basis
 33basis = photon_basis(spin_basis_1d, L=1, Nph=Nph_tot)
 34# compute atom-photon Hamiltonian H
 35H = hamiltonian(static, dynamic, dtype=np.float64, basis=basis)
 36#
 37##### set up semi-classical Hamiltonian #####
 38# define operators
 39dipole_op = [[A, 0]]
 40
 41
 42# define periodic drive and its parameters
 43def drive(t, Omega):
 44    return np.cos(Omega * t)
 45
 46
 47drive_args = [Omega]
 48# define semi-classical static and dynamic lists
 49static_sc = [["z", at_energy]]
 50dynamic_sc = [["x", dipole_op, drive, drive_args]]
 51# compute semi-classical basis
 52basis_sc = spin_basis_1d(L=1)
 53# compute semi-classical Hamiltonian H_{sc}(t)
 54H_sc = hamiltonian(static_sc, dynamic_sc, dtype=np.float64, basis=basis_sc)
 55#
 56##### define initial state #####
 57# define atom ground state
 58# psi_at_i=np.array([1.0,0.0]) # spin-down eigenstate of \sigma^z in QuSpin 0.2.3 or older
 59psi_at_i = np.array(
 60    [0.0, 1.0]
 61)  # spin-down eigenstate of \sigma^z in QuSpin 0.2.6 or newer
 62# define photon coherent state with mean photon number Nph
 63psi_ph_i = coherent_state(np.sqrt(Nph), Nph_tot + 1)
 64# compute atom-photon initial state as a tensor product
 65psi_i = np.kron(psi_at_i, psi_ph_i)
 66#
 67##### calculate time evolution #####
 68# define time vector over 30 driving cycles with 100 points per period
 69t = Floquet_t_vec(Omega, 30)  # t.i = initial time, t.T = driving period
 70# evolve atom-photon state with Hamiltonian H
 71psi_t = H.evolve(psi_i, t.i, t.vals, iterate=True, rtol=1e-9, atol=1e-9)
 72# evolve atom GS with semi-classical Hamiltonian H_sc
 73psi_sc_t = H_sc.evolve(psi_at_i, t.i, t.vals, iterate=True, rtol=1e-9, atol=1e-9)
 74#
 75##### define observables #####
 76# define observables parameters
 77obs_args = {"basis": basis, "check_herm": False, "check_symm": False}
 78obs_args_sc = {"basis": basis_sc, "check_herm": False, "check_symm": False}
 79# in atom-photon Hilbert space
 80n = hamiltonian([["|n", [[1.0]]]], [], dtype=np.float64, **obs_args)
 81sz = hamiltonian([["z|", [[1.0, 0]]]], [], dtype=np.float64, **obs_args)
 82sy = hamiltonian([["y|", [[1.0, 0]]]], [], dtype=np.complex128, **obs_args)
 83# in the semi-classical Hilbert space
 84sz_sc = hamiltonian([["z", [[1.0, 0]]]], [], dtype=np.float64, **obs_args_sc)
 85sy_sc = hamiltonian([["y", [[1.0, 0]]]], [], dtype=np.complex128, **obs_args_sc)
 86#
 87##### calculate expectation values #####
 88# in atom-photon Hilbert space
 89Obs_t = obs_vs_time(psi_t, t.vals, {"n": n, "sz": sz, "sy": sy})
 90O_n, O_sz, O_sy = Obs_t["n"].real, Obs_t["sz"].real, Obs_t["sy"].real
 91# in the semi-classical Hilbert space
 92Obs_sc_t = obs_vs_time(psi_sc_t, t.vals, {"sz_sc": sz_sc, "sy_sc": sy_sc})
 93O_sz_sc, O_sy_sc = Obs_sc_t["sz_sc"].real, Obs_sc_t["sy_sc"].real
 94##### plot results #####
 95import matplotlib.pyplot as plt
 96import pylab
 97
 98# define legend labels
 99str_n = "$\\langle n\\rangle,$"
100str_z = "$\\langle\\sigma^z\\rangle,$"
101str_x = "$\\langle\\sigma^x\\rangle,$"
102str_z_sc = "$\\langle\\sigma^z\\rangle_\\mathrm{sc},$"
103str_x_sc = "$\\langle\\sigma^x\\rangle_\\mathrm{sc}$"
104# plot spin-photon data
105fig = plt.figure()
106plt.plot(t.vals / t.T, O_n / Nph, "k", linewidth=1, label=str_n)
107plt.plot(t.vals / t.T, O_sz, "c", linewidth=1, label=str_z)
108plt.plot(t.vals / t.T, O_sy, "tan", linewidth=1, label=str_x)
109# plot semi-classical data
110plt.plot(t.vals / t.T, O_sz_sc, "b", marker=".", markersize=1.8, label=str_z_sc)
111plt.plot(t.vals / t.T, O_sy_sc, "r", marker=".", markersize=2.0, label=str_x_sc)
112# label axes
113plt.xlabel("$t/T$", fontsize=18)
114# set y axis limits
115plt.ylim([-1.1, 1.4])
116# display legend horizontally
117plt.legend(loc="upper right", ncol=5, columnspacing=0.6, numpoints=4)
118# update axis font size
119plt.tick_params(labelsize=16)
120# turn on grid
121plt.grid(True)
122# save figure
123plt.tight_layout()
124plt.savefig("example3.pdf", bbox_inches="tight")
125# show plot
126# plt.show()
127plt.close()