Rabi Oscillations

Rabi Oscillations#

Author:
Date:
Time spent on this assignment:

import numpy as np 
import pylab as plt
import qutip
import scipy
import scipy.linalg
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib.animation as animation

def resetMe(keepList=[]):
    ll=%who_ls
    keepList=keepList+['FuncAnimation','HTML','resetMe','scipy','np','plt','math','jax','jnp','jit','grad','HTML','animation','qutip','animation','FuncAnimation']
    for iiii in keepList:
        if iiii in ll:
            ll.remove(iiii)
    for iiii in ll:
        jjjj="^"+iiii+"$"
        %reset_selective -f {jjjj}
    ll=%who_ls
    plt.rcParams.update({"font.size": 14})
    return
resetMe()
import datetime;datetime.datetime.now()

def PlotBloch(v, bb, makeVector=False):
    vv = np.asarray(v)
    a = vv[0] * np.conj(vv[0])
    b = vv[1] * np.conj(vv[0])
    x = np.real(2.0 * b.real)
    y = np.real(2.0 * b.imag)
    z = np.real(2.0 * a - 1.0)

    if makeVector:
        old_marker = bb.point_marker
        old_size = bb.point_size
        old_point_color = bb.point_color

        bb.vector_color = [] #clear the vector color list.
        bb.vector_color = ['green'] * (len(bb.vectors) + 1)  # Ensure the list has enough elements

        #bb.add_vectors([x, y, z])
        bb.add_vectors([x, y, z])  # Pass color to add_vectors

        bb.point_marker = ['*']
        bb.point_color = ['blue']

        bb.add_points([x, y, z])
        bb.point_marker = old_marker
        bb.point_size = old_size
        bb.point_color = old_point_color
    else:
        old_point_color = bb.point_color
        bb.add_points([x, y, z])
        bb.point_color = old_point_color

    return



def SetupBloch():
  fig = plt.figure()
  
#  ax = fig.add_subplot(111, projection='3d')
  ax = fig.add_subplot(111,projection="3d")
  ax.view_init(-40, -15) 


  ax.set_aspect('equal')
  ax.set_axis_off()
  bb = qutip.Bloch(axes=ax)
  #bb.point_color = ['red','orange','yellow','green','blue','purple']
  #bb.point_color = ['red']  
  bb.point_marker = ['o']
  bb.vector_color = ['green']
  #bb.point_size = [15]
  bb.point_mode = ['markers']
  bb.frame_color = 'black'
  bb.frame_width = 1
  bb.frame_alpha = 0.5
  bb.sphere_alpha = 0.1
  bb.vector_width = 3
  bb.alpha = 0.1
  return bb



def animate(i, states, sphere,plotPath=False):
  
  sphere.clear()
  #vv = states[i]
  for vv in states[0:i]:
    PlotBloch(vv, sphere, makeVector=False)
  PlotBloch(states[i], sphere, makeVector=True)
  return sphere.render()




sigma_x = qutip.sigmax().full()
sigma_y= qutip.sigmay().full()
sigma_z = qutip.sigmaz().full()
hbar=1

Exercise 1. Time Dependence in the Schrodinger Picture#

a. The Hamiltonian#

In this assignment, we will consider time-dependent Hamiltonians of the form

\[ H = H_0 + V(t) \]

where \(H_0\) has no time-dependence and \(V(t)\) depends on time.

Let us define

\[ H_0 = E_0 |0\rangle \langle 0| + E_1 |1|\rangle \langle 1| \]

where \(E_0 = -0.4\) and \(E_1=0.6\)

and

\[ V(t) = \delta/2 \exp[i \omega t] |0\rangle \langle 1| + \delta/2 \exp[-i \omega t] |1 \rangle \langle 0| \]

Define functions build_H0() and build_V(t,delta,omega) which returns these Hamiltonians. Notice that they are both 2x2 matrices.

###ANSWER HERE

b. Time Evolution#

The most straightforward way to implement time-evolution is to directly use the time-dependent Schrodinger equation

\[ i \hbar \frac{\partial \psi }{\partial t} = H \psi \]

which we can turn into a dynamical equation of motion as

\[ \psi(t + \delta t) = e^{-i\hbar H(t) \delta t} \psi(t) \]

You can assume throughout this assignment that we are working in units where \(\hbar=1\)

Write a function SchrodingerTimeEvolution(T,dt,omega,delta) which takes a total time \(T\) and a time step \(dt\) and the values of \(\delta\) and \(\omega\) for the driving potential. It should then perform Schrodinger evolution starting with a configuration that is in state \(|0\rangle\) (i.e. state = np.array([1,0])). It should return a list of times (ts=np.arange(0,T,dt)) and the state (e.g. 2x1 vectors) at each time-step.

To do this, you want to

Define the Hamiltonian \(H_0\). This should be a 2x2 matrix.
Define the initial state \(\Psi(t) = [1,0]\)
Looping over the time-steps
- Get the matrix \(H(t)=H_0 + V(t)\)
- Define the time-evolution matrix \(M(t) \equiv e^{-i H \delta t}\) - i.e. scipy.linalg.expm(H_t)
- Apply \(M(t)\) to the current state to get the new state.
- Store in a list the current state.

Call your function with

T=17.0
\(dt\) = 0.001
\(\omega = 1.3\)
\(\delta = 0.2\)

You are then going to want to plot two things:

(1) The probability that the state is in one - i.e. np.abs(state[1])**2 as a function of time.

(2) The location of the state on the bloch sphere. To do this you can do

bb=SetupBloch()
for state in states[::100]:
    PlotBloch(state,bb,makeVector=False)
bb.render()
plt.show()

(3) Animate it. To do that you can do

sphere=SetupBloch()
anim = FuncAnimation(plt.gcf(), animate, frames=len(states[::100]), fargs=(states[::100], sphere,True), blit=False,interval=50, repeat=False)
display(HTML(anim.to_jshtml()))
plt.close()

###ANSWER HERE

###ANSWER HERE

###ANSWER HERE

b. The Interaction Picture#

Given the states represented in the Schrodinger picture, it’s possible to convert them into the interaction picture. The interaction picture essentially takes away the rotation due to the non-interacting Hamiltonian. The way to make this conversion is to do

\[ \Psi_\textrm{interaction} = \exp\left(\frac{i}{\hbar} t H_0 \right)\Psi_\textrm{schrodinger} \]

Write a function S2I(psi,t,H0) which takes the non-interacting Hamiltonian and the state in the Schrodinger picture and gives back the state in the interaction picture. Use this function to convert all the Schrodinger picture states into Interaction picture states. Then plot (and animate) the interaction picture states on the Bloch sphere. What is different about these states?

###ANSWER HERE

###ANSWER HERE

Exercise 2. Time evolving in the Interaction Picture#

a. Time evolving in the interaction picture#

In the previous problem, we time-evolved in the Schrodinger picture. In this problem we are going to time-evolve directly in the interaction picture. First let us derive this.

Consider a Hamiltonian

\[ H = H_0 + V(t) \]

where

\[ |\psi(t) \rangle_I = e^{iH_0 t/\hbar} |\psi(t)\rangle_S \]

where \(|\psi(t)\rangle_S\) is the time-dependent wave-function in the Schrodinger picture and \(|\psi(0)\rangle_S = \psi(0)\rangle_I\)

Then we find that

\[ i\hbar \partial_t |\psi(t)\rangle_I = V_I(t) |\psi(t)\rangle_I \]

where

\[ V_I(t) = e^{iH_0 t /\hbar}V e^{-i H_0 t/\hbar} \]

Solving, we get a series of first order differential equations

\[ i\hbar \dot{c_m}(t) = \sum_n V_{mn}(t) e^{i\omega_{mn}t}c_n(t) \]

where \(\omega_{mn} \equiv (E_m-E_n)/\hbar\)

which we can now write as

\[ c_m(t+dt)= c_m(t) + dt \sum_n V_{mn}(t) e^{i\omega_{mn}t}c_n(t) \]

Now write a function InteractionTimeEvolution(T,dt,omega,delta) which uses this equation of motion to time-evolve in the interaction picture returning both the list of times and the list of states. Using the PlotBloch function go ahead and plot this evolution on the Bloch sphere.

###ANSWER HERE

###ANSWER HERE

|

b. Transforming to the Schrodinger Picture#

Now go ahead and write a function I2S(psi,t,H0) to transform an interaction-picture state to a schrodinger-picture state. Then use your function to convert all the states generated in part (a) to the Schrodinger picture. After you do this conversion,

plot the probability of the state being in \(|1\rangle\).
the state on the Bloch sphere

###ANSWER HERE

c. Probability to reach \(|1\rangle\)#

Let’s measure the maximum probability that you get to \(|1\rangle\) as a function of \(\Delta \omega \equiv \omega - \omega_{10}\) - i.e. omega - (E[1]-E[0])/hbar. To do this, you want to let the maximum \(T=50\).

Change \(\omega\) from 0.1 to 1.9 in small incremements and then plot the maximum probability to get \(|1\rangle\) versus \(\Delta \omega\).

Also plot the probability of being at \(|1\rangle\) at \(T=8.3\) as a function of \(\Delta \omega\)

###ANSWER HERE

###ANSWER HERE

Exercise 3. Time dependent perturbation theory#

a. Probability (at first order) to jump to state \(|1\rangle\)#

We are now going to evaluate these probabilities using time-dependent perturbation theory.

Time-dependent perturbation theory tells us that, to first order, the probability of going from state \(m \rightarrow n\) is

\[ P_{m\leftarrow n }^{(1)}(t) = \left|\int_0^{t}e^{i\omega_{mn}t'}\frac{V_{mn}(t')}{i\hbar}dt'\right|^2 \]

Write a function Prob_10(T,dt,omega,delta) which gives the first order probability for going from \(0 \rightarrow 1\) returning both the times and these probabilities. You may find it useful to use np.cumsum to do the integral out to different \(t\).

Using \(T=40\), graph that probability as a function of time and plot it on top of the actual probabilities.

Print out \(T=2\pi/\Delta \omega\) and \( | |V|/\hbar \Delta \omega|^2 \) and see that this corresponds to the wave-length and height of the first-order perturbation terms.

###ANSWER HERE

b. Probability (at first order) to jump to state \(|1\rangle\) using the formula#

In exercise (a), we computed the first order perturbative correction using a general formula where we needed to take an integral. For this specific problem, Griffith’s has done the integral analytically and worked out a general formula,

\[ P_{10} = -i \frac{\delta}{\hbar} \frac{\sin(\Delta \omega t/2)}{\Delta \omega } e^{i\Delta \omega t/2} \]

Write a function Prob_10_formula(T,dt,omega,delta) that generates these probabilities using the formula and make sure it matches your results from exercise (a).

###ANSWER HERE

c. Maximum Probabability and Probability at \(T=8.3\)#

At first order, let’s measure the maximum probability that you get to \(|1\rangle\) as a function of \(\Delta \omega \equiv \omega - \omega_{10}\) - i.e. omega - (E[1]-E[0])/hbar.

Change \(\omega\) from 0.1 to 5.0 in small incremements and then plot the probability versus \(\Delta \omega\).

Also plot the probability of being at \(|1\rangle\) at \(T=8.3\) as a function of \(\Delta \omega\) Compare it against the exact answer by plotting it on the same plot.

###ANSWER HERE

###ANSWER HERE