Homework on atom-photon interaction: the Jaynes-Cummings model ¶
General information¶
- the deadline for this homework: 23/12/18 @10:00pm Paris time (date of email)
- must be sent by mail to guillaume.roux@u-psud.fr in .html format only (single file)
- contact me in case of trouble, I should answer within 24h.
- the notebook source for this file can be downloaded at: http://lptms.u-psud.fr/membres/groux/Test/Homework-HO.ipynb
The subject looks terribly long because there actually is a lot of reading. There are a few questions to check your understanding and a few (mostly at the end) on which you have to write a bit of code. All the classes are prefilled so that you have use them for the last question.
All codes are written in Python 3 so check your Kernel.
import numpy as np
In this homework, we would like to write an exact diagonalization code that handles two-level systems (this is exactly equivalent of what has been done for spin 1/2) and harmonic oscillators and their product states. For instance, this would allow us to describe a two-level atom $\{|\sigma\rangle \}=\{|g\rangle ,|e\rangle \}$ coupled to a quantized electromagnetic mode $\{|n\rangle\}$ of a given frequency $\omega$, through the tensor product basis $\{|\sigma,n\rangle\}$. More generally, we would like to account for several ''spin-string'' systems $\{|\sigma_1,n_1,\sigma_2,n_2,\ldots\rangle\}$ or a single atom coupled to many modes $\{|\sigma,n_1,n_2,n_3,\ldots\rangle\}$.
The main work is to reconsider how to write the Hilbert space (configuration space) in a more generic way and how to apply operators on systems with inequivalent degrees of freedom.
We apply this to the Jaynes-Cummings Hamiltonian and benchmark the results with the exactly known time-evolution. See https://en.wikipedia.org/wiki/Jaynes%E2%80%93Cummings_model for the physics of the Jaynes-Cummings model.
After that, we will see how to perform time-evolution with sparse matrices in order to tackle larger Hilbert spaces.
Local Hilbert spaces and local operators¶
We start by describing a single degree of freedom by defining a generic localHilbert space object that contains all the list of possible quantum numbers corresponding to this degree of freedom (hence we derive it from the python list container) and add explicit names in the dictionary self.states for printing purposes.
Two-level atom¶
We give below the example of the two-level local degree of freedom.
class localHilbert(list):
"describing a local degree of freedom"
def __init__(self,*args):
list.__init__(self, *args)
self.states = {}
def __repr__(self):
return self.states.__repr__()
class TwoLevel(localHilbert):
"Two level degree of freedom"
def __init__(self):
localHilbert.__init__(self, [0,1])
self.states[0] = 'g'
self.states[1] = 'e'
print(TwoLevel())
We now would like to apply a local operator $O_i$ on the degree of freedom at site $i$. We use exactly the same general definition as in the lecture for the application of $O_i$ on a wave-vector $|\psi\rangle$, with the object psi containing a reference to the Hilbert object that provides the dictionnary of indices for each configuration.
class Operator:
"Generic class to perform O|psi>"
def __mul__(self,psi):
res = zeroWf(psi.hilbert,dtype=psi.dtype)
for conf, i in psi.hilbert.Confs.items():
coef, newconf = self.Apply(conf)
j = psi.hilbert.Confs[newconf]
res[j] += coef * psi[i]
return res
def __str__(self):
return self.name
Furthermore, we assume that the specialization of such a local operator also contains a reference to an hilbert object that has the two following methods:
- hilbert.readSite(i, conf): reads the value of the quantum number at site i in conf
- hilbert.changeSite(n, i, conf): changes the quantum number to value n at site i in conf
We give as an example the specialization for the two level atoms using the following definitions of operators in the $\{|g\rangle ,|e\rangle \}$ basis:
$$\sigma^Z = \begin{pmatrix} -1 & 0 \\ 0 & 1 \end{pmatrix} \qquad \sigma^+ = \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} \qquad \sigma^- = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} \qquad $$ We recall that the Apply retruns the matrix coefficient and the new generated configuration (if not diagonal).
class sigmaZ(Operator):
def __init__(self, hilbert, site):
self.site = site
self.hilbert = hilbert
self.name = "sigZ_" + str(self.site)
def Apply(self, conf):
s = self.hilbert.readSite(self.site, conf)
return 2*s-1, conf
class sigmaPlus(Operator):
def __init__(self, hilbert, site):
self.site = site
self.hilbert = hilbert
self.name = "sig+_" + str(self.site)
def Apply(self, conf):
s = self.hilbert.readSite(self.site, conf)
if s == 0 :
newconf = self.hilbert.changeSite(1, self.site, conf)
return 1, newconf
else :
return 0.0, conf
class sigmaMinus(Operator):
def __init__(self, hilbert, site):
self.site = site
self.hilbert = hilbert
self.name = "sig-_" + str(self.site)
def Apply(self, conf):
s = self.hilbert.readSite(self.site, conf)
if s == 1 :
newconf = self.hilbert.changeSite(0, self.site, conf)
return 1, newconf
else :
return 0.0, conf
Harmonic oscillator¶
The local hilbert space corresponding to an harmonic oscillator. As the dimension of the local hilbert space is in principle infinite, we introduce a truncation by setting $N_{\text{max}}$ the maximum value of the occupation number $n$.
class Harmonic(localHilbert):
"Harmonic oscillator degree of freedom"
def __init__(self,Nmax=2):
localHilbert.__init__(self,range(Nmax+1))
self.states = { i:str(i) for i in range(Nmax+1) }
ho = Harmonic(5)
print(ho)
N.-B.: We assume that the hilbert object provides the dimension of the local hilbert space at site i through the access to hilbert.shape[i] from which one can extract the local value of Nmax at site i.
1/ Fill the basic definitions of operators for the harmonic oscillator. We recall
$$\hat{a}^{\dagger}|n\rangle = \sqrt{n+1}|n+1\rangle \qquad \hat{a}|n\rangle = \sqrt{n}|n-1\rangle \qquad \hat{N}|n\rangle = n|n\rangle $$
class Nop(Operator):
def __init__(self,hilbert,site):
self.site = site
self.hilbert = hilbert
self.name = "N_"+str(site)
def Apply(self,conf):
return ???, conf
class Aop(Operator):
def __init__(self, hilbert, site):
self.site = site
self.hilbert = hilbert
self.minn = 0
self.name = "A_" + str(self.site)
def Apply(self, conf):
n = self.hilbert.readSite(self.site,conf)
if n > self.minn :
newconf = ???
return np.sqrt(n), newconf
else :
return 0.0, conf
class Adagop(Operator):
def __init__(self, hilbert, site):
self.site = site
self.hilbert = hilbert
self.maxn = self.hilbert.shape[site] - 1
self.name = "Adag_" + str(self.site)
def Apply(self, conf):
n = self.hilbert.readSite(self.site, conf)
if n < self.maxn :
newconf = ???
return np.sqrt(n+1), newconf
else :
return 0.0, conf
Product of local operators¶
We will have to deal with operators that are products of local operators, for instance of the form $\hat{a}\sigma^+$
class OpTerm(list):
"product of local operators"
def __mul__(self,psi):
res = zeroWf(psi.hilbert,dtype=psi.dtype)
res[:] = psi[:]
for op in self[::-1]: res = op*res
return res
def __str__(self):
s = self[0].name
for op in self[1:]: s += "*"+op.name
return s
The example below will work after you have declared the Hilbert class just below.
H = Hilbert( [ TwoLevel(), Harmonic(5) ])
print(OpTerm([Aop(H,1),sigmaPlus(H,0)]))
print(OpTerm([Adagop(H,1),Adagop(H,1),Aop(H,1),sigmaPlus(H,0)]))
Full Hilbert space construction¶
We now turn to combining various degrees of freedom by making the tensor product of local Hilbert spaces. As in the lecture, we are using integers to store the configurations but describing a local site now requires more than one bit if the local dimension is greater than two.
2/ Let $d$ be the dimension of a local hilbert space. How many bits are required to describe it? Explain how the method NumBitsForASite(self,dimension) works in the Hilbert class below
The idea is then to construct an Hilbert space by giving the list of local hilbert spaces (for instance [TwoLevel(), Harmonic(5),...]) and allocating the required number of bits to each degree of freedom from the right to the left in the memory word (integer) corresponding to a configuration.
3/ Several useful objects are stored in the Hilbert class. Please, explain what are and what are their use for the following attributes by looking at the setConfigurationMask() method and changeSite() and readSite() methods (also have a look at the cell below that contains an example to play with):
masksize, shifts, sitemasks, sitenegativs
4/ Explain how the Hilbert space is created in the create method. N.-B.: the conf2int method converts a configuration in list format into the binary configuration.
class Hilbert:
def __init__(self,description=[ range(2) ]):
"Description is a list containing the local hilbert space at site i"
self.description = description
self.Nsites = len(self.description)
self.shape = tuple();
for loc in self.description: self.shape += (len(loc),)
self.hilbertsize = 1
for h in self.shape: self.hilbertsize *= h
self.setConfigurationMask()
self.create()
def NumBitsForASite(self,dimension):
import math
"number of bits required for local hilbert space dimension"
if dimension < 2:
print ("Warning, number of bits will be zero")
return int(math.ceil(math.log(dimension,2)))
def setConfigurationMask(self):
self.masksize = tuple()
for loc in self.description:
self.masksize += (self.NumBitsForASite(len(loc)),)
self.wordsize = sum(self.masksize)
# shifts
self.shifts = tuple(); shift_ = 0
for m in self.masksize:
self.shifts += (shift_,)
shift_ += m
# masks
self.sitemasks = tuple();
for i in range(self.Nsites):
self.sitemasks += ( (2**self.masksize[i]-1)<<self.shifts[i], )
# negative masks
fullmask = (2**63-1)<<1|1
self.sitenegativs = tuple()
for i in range(self.Nsites):
self.sitenegativs += ( fullmask^self.sitemasks[i], )
def binconf(self,c):
"binary representation of a conf integer"
return np.binary_repr(c,self.wordsize)
def ket(self,conf):
if isinstance(conf,int):
L = len(self.description)
c = self.int2conf(conf)
elif isinstance(conf,list):
c = conf
else:
print ("Error in ket")
s = "|"
#if not len(c) == len(self.description): return "Error"
for l in range(len(c)):
s += self.description[l].states[c[l]]+","
return s[:-1]+">"
def conf2int(self,conf):
"conversion from list representation to binary [1,0,3,0] --> 10110"
r = 0
for i in range(self.Nsites-1, -1, -1):
r = (r<<self.masksize[i]) | conf[i]
return r
def int2conf(self,intconf):
"conversion from binary representation to list 10110 --> [1,0,3,0]"
return [ self.readSite(i,intconf) for i in range(self.Nsites) ]
def readSite(self,i,conf):
"read quantum number on site i"
return (conf&self.sitemasks[i])>>self.shifts[i]
def changeSite(self,qn,i,conf):
"set quantum number qn on site i, does not modify conf"
r = (qn)<<self.shifts[i]
return ( conf & self.sitenegativs[i] ) | r
def create(self, constraint = None):
"builds up the Hilbert space"
from copy import deepcopy
conf_list = [ (loc,) for loc in self.description[0] ]
for local in self.description[1:]:
tmp = []
for conf in conf_list:
for s in local:
tmp.append(conf+(s,))
conf_list = deepcopy(tmp)
self.Confs = { self.conf2int(conf_list[i]):i for i in range(len(conf_list)) }
def index(self,conf):
"returns the index of a configuration"
if isinstance(conf,list):
return self.Confs[self.conf2int(conf)]
else:
return self.Confs[conf]
def __repr__(self):
"prints the first 30 confs of hilbert space"
from itertools import islice
num = min(30,self.hilbertsize)
s = "Printing first "+str(num)+" configurations among "+str(self.hilbertsize)+"\n"
s += ", ".join([self.ket(c) for c in list(islice(self.Confs.keys(),num))])
return s
hilbert = Hilbert( [TwoLevel(),Harmonic(5)] )
print ("Description", hilbert.description)
print ("Shape", hilbert.shape)
print ("Hilbert size", hilbert.hilbertsize)
print ("Wordsize", hilbert.wordsize)
print ("Mask sizes",hilbert.masksize)
print ("Mask shifts",hilbert.shifts)
print ("The masks", [ hilbert.binconf(m) for m in hilbert.sitemasks ])
print ("The negatives", [ bin(m)[2:] for m in hilbert.sitenegativs ])
print(hilbert)
Wave function¶
The wave-function class is the same as before, up to minor changes in methods names. There's nothing new below, just shit-enter the cell.
def cpp_convert(x,p):
if type(x) is complex or type(x) is np.complex64 or type(x) is np.complex128:
return "("+str(round(x.real,p))+","+str(round(x.imag,p))+")"
else:
return str(round(x,p))
class Wavefunction(np.ndarray):
def __new__(cls,hilbert,dtype=float):
obj = np.ndarray.__new__(cls,(hilbert.hilbertsize,),dtype)
obj.hilbert = hilbert
return obj
def __array_finalize__(self, obj):
if obj is None: return
self.hilbert = getattr(obj,'Hilbert',None)
def setvalue(self, conf, value):
"Sets the coefficient in front of conf"
self[self.hilbert.index(conf)] = value
def getvalue(self, conf):
"gets the coefficient in front of conf"
return self[self.hilbert.index(conf)]
def readable(self,precision=4,eps=1e-4):
from itertools import islice
num = min(30,self.hilbert.hilbertsize)
s = ""
for conf,i in islice(self.hilbert.Confs.items(),num):
coef = self[i]
if abs(coef)>eps:
symb = " + "
if type(coef) is np.float64 and np.sign(coef) == -1.0:
symb = " - "
coef = abs(coef)
s += symb + cpp_convert(coef,precision) + self.hilbert.ket(conf)
if s[:3] == " + ": s = s[3:]
return "|psi> = "+s
def norm(self):
return np.linalg.norm(self)
def normalize(self):
norm = self.norm()
try: self /= norm
except: print("can't normalize a vector of norm", norm)
def randomize(self,a=-1.0,b=1.0):
self[:] = (b-a)*np.random.random_sample(self.size).reshape(self.shape)+a
if self.dtype == np.complex128:
self.imag = (b-a)*np.random.random_sample(self.size).reshape(self.shape)+a
self.normalize()
def measure(self,op,hermitic=True):
"computes mean-value of an hermitic operator"
res = scalar(self,op*self)
if abs(res.imag) > 1e-15:
print("Warning, issue with hermiticity")
return res.real
def scalar(psi1,psi2):
"computes overlaps of psi1 and psi2 with complex conjugation on psi1"
return np.vdot(psi1,psi2)
def zeroWf(hilbert,dtype=float):
psi = Wavefunction(hilbert,dtype=dtype)
psi[:] = 0.0
return psi
def basisvectorWf(hilbert,index,dtype=float):
psi = zeroWf(hilbert,dtype=dtype)
psi[index] = 1.0
return psi
def randomWf(hilbert,dtype=float):
psi = Wavefunction(hilbert,dtype=dtype)
psi.randomize()
return psi
def copy(psi):
res = Wavefunction(psi.hilbert,dtype=psi.dtype)
res[:] = psi[:]
return res
psi = randomWf(hilbert)
print(psi.readable())
Creating a coherent state¶
A coherent state $|\alpha\rangle$ (https://en.wikipedia.org/wiki/Coherent_states) is the eigenvector of the $\hat{a}$ operator corresponding to $\alpha \in \mathbb{C}$, ie. $$\hat{a}|\alpha\rangle =\alpha |\alpha\rangle$$ and can be expanded in series in the Fock basis as $$|\alpha\rangle = e^{-|\alpha|^2/2} \sum_{n=0}^{\infty} \frac{\alpha^n}{\sqrt{n!}} |n\rangle$$ we take $\alpha$ real for simplicity in the following. We give below a function that computes normalized coefficients for the coherent state in a truncated local Hilbert space.
def coherentStateCoef(alpha,Nmax):
coef = np.ones(Nmax+1)
for n in range(1,Nmax+1):
coef[n] = alpha/np.sqrt(n)*coef[n-1]
return coef/np.linalg.norm(coef)
print(coherentStateCoef(2.0,10))
hilbert = Hilbert( [Harmonic(10)] )
def coherentState(hilbert,i,alpha):
"create a coherent state alpha for local hilbert space of site i"
nmax = hilbert.shape[i]-1
psi = zeroWf(hilbert,dtype=float)
psi[:] = coherentStateCoef(alpha,nmax)
return psi
print(coherentState(hilbert,0,2.0).readable())
Creating product initial states for time evolution¶
5/ Fill up the following function that must create the following product states: $|g,\alpha\rangle$ and $|e,\alpha\rangle$ depending on the level variable.
def initialState(hilbert,level=0,alpha=2.0):
"create a coherent state alpha for local hilbert space of site i"
nmax = hilbert.shape[1]-1
psi0 = zeroWf(hilbert,dtype=float)
coef = coherentStateCoef(alpha,nmax)
for n in range(nmax+1):
psi0.setvalue(????)
return psi0
hilbert = Hilbert( [TwoLevel(),Harmonic(10)] )
print(initialState(hilbert,1,2.0).readable())
General operators and Hamiltonian¶
We use a simplified version of the GeneralOp seen in classroom that only handles full matrices (non-sparse) as we will treat small Hilbert spaces.
class GeneralOp(list):
"""
A list of 2-element in the format [ [coef, Operator] ]
corresponding to \sum_<i,j> coef_ij * Op_ij
warning: assumes hermitic operators
"""
def __init__(self,iterable=[]):
list.__init__(self,iterable)
self.having_matrix = False
self.storage = 0
def __str__(self):
if not len(self): return ""
s = "("+str(self[0][0])+")*"+str(self[0][1])
for coef, op in self[1:]:
s += "+("+str(coef)+")*"+str(op)
return s
def __iadd__(self,other):
return GeneralOp(self[:]+other[:])
def __add__(self,other):
return GeneralOp(self[:]+other[:])
def __mul__(self,psi):
res = copy(psi)
if self.having_matrix:
res[:] = np.dot(self.matrix,psi)
else:
res[:] = 0.0
for coef, op in self:
res += coef * (op*psi)
return res
def createMatrix(self,hilbert,dtype=float):
if self.having_matrix: return
if hilbert.hilbertsize < 2: return
size = hilbert.hilbertsize
self.matrix = np.zeros(shape=(size,size), dtype = dtype)
for i in np.arange(size):
self.matrix[:,i] = self*basisvectorWf(hilbert,i)
self.having_matrix = True
def diagonalize(self,hilbert,with_eigvec=False):
""""Assumes an hermitic operator by default and that the matrix has been created
stores the energies and eigenstates in self.En and self.Phin"""
if not self.having_matrix:
self.createMatrix(hilbert,dtype)
if self.having_matrix:
if with_eigvec:
self.En, self.Phin = np.linalg.eigh(self.matrix)
else:
self.En = np.linalg.eigvalsh(self.matrix)
def Energies(self,hilbert,dtype=float):
"computes num energies and returns them."
if not self.having_matrix:
self.createMatrix(hilbert,dtype)
self.diagonalize(hilbert)
return self.En
Matrices of the $\hat{a}$ operator¶
6/ Fill up the following code that must display the matrix corresponding to the $\hat{a}$ operator.
H = Hilbert( [Harmonic(10)] )
a = GeneralOp([ [1.0,????] ] )
a.????
np.set_printoptions(precision=3)
print(a.matrix)
The Jayne-Cummings Hamiltonian¶
The atom-field Hamiltonian of the Jaynes-Cummings model reads:
$$ H = \frac{\omega_a}{2} \sigma^Z + \omega_c \hat{N} + g(\sigma^+\hat{a}+\sigma^-\hat{a}^{\dagger})$$
7/ Give a simple physical interpretation of the three terms in the Hamiltonian.
8/ Fill the lines below to create the Hamiltonian.
Hint look back up at the use of the OpTerm class
Nmax = 50
H = Hilbert( [TwoLevel(),Harmonic(Nmax)] )
omega_a = 1.0
omega_c = 1.0
g = 0.1
Ham = GeneralOp()
Ham += [ [ omega_c, Nop(H,1) ] ]
Ham += [ [ g, ???? ] ]
Ham += [ [ g, ????] ]
Ham += [ [ 0.5*omega_a, ???? ] ]
Below are lines that help you compare the energies obtained from your Hamiltonian to the exact ones.
np.set_printoptions(precision=10)
print(Ham.Energies(H)[:40])
delta = omega_c-omega_a
Exact = [ -0.5*omega_c ]
for n in range(1,Nmax+2):
OmegaN = 2*g*np.sqrt(n)
Exact.append( (n-1)*omega_c + 0.5*omega_a +0.5*(delta + np.sqrt(delta**2+OmegaN**2)) )
Exact.append( (n-1)*omega_c + 0.5*omega_a +0.5*(delta - np.sqrt(delta**2+OmegaN**2)) )
En_ex = np.array(sorted(Exact))
print(En_ex[:40])
Time-evolution using exact diagonalization¶
Starting from an initial state, one can follow the time evolution of the system.
9/ Fill the lines below to perform time evolution through the explicit calculation of the evolution operator $U$.
10/ Fill also the lines that allow you to compute the probability that the atom is in the excited state, given by
$$ P_e(t) = \sum_{n=0}^{N_{max}} |\langle{e,n}|{\psi(t)}\rangle|^2$$
Ham.diagonalize(H,with_eigvec=True)
Ham.InvPhin = np.linalg.inv(Ham.Phin)
psi0 = initialState(H,level=1,alpha=5.0)
psi_t = zeroWf(H, dtype=complex)
Pe = []
times = np.linspace(0,40/g,1000)
for t in times:
D = np.diag(np.exp((-1j*t)*Ham.En))
U = np.dot(?????)
psi_t[:] = np.dot(U,psi0)
probas = np.array([ psi_t.????? for n in range(Nmax+1) ])
Pe.append(np.linalg.norm(probas)**2)
The lines below allow you to compare your results with the exact solution
#exact solution
pb0 = np.array([ psi0.getvalue([1,n])**2 for n in range(Nmax+1) ])
def Pe_ex(t):
somme = 0.0
for n in range(Nmax+1):
somme += pb0[n]*(1+np.cos(2*g*t*np.sqrt(1+n)))
return somme/2
from pylab import *
%matplotlib inline
plot(g*times,Pe,'k-')
plot(g*times,Pe_ex(times),'r-')
xlabel("$gt$", fontsize=20)
ylabel("$P_e(t)$", fontsize=20)
show()
plot(g*times,Pe-Pe_ex(times),'k-')
xlabel("$gt$", fontsize=20)
ylabel("Error", fontsize=20)
show()