2
$\begingroup$

I'm using Qutip to plot some basic two level dynamics using hamiltonians with a temporal envelope defined as the sum of two error functions, designed to make it more representative of experimental conditions.

I'm using the following code.

import numpy as np
from qutip import * 
from scipy.special import erf
import matplotlib.pyplot as plt

def pulseEdge(x, x0, sigma):
    pulseEdge = ((erf((x-x0)/sigma)) +1.)/2.
    return pulseEdge

def pulse(t, tOn, tOff, sigma):
    return pulseEdge(t, tOn, sigma)* (1.-pulseEdge(t, tOff, sigma))

def H(t, args):
    H = Qobj([[0.,1],[1,0.0]])
    envelope = pulse(t, 2,2+np.pi/2,0.2)
    return H*envelope

tlist = np.linspace(0,10,1000)
h = np.array([H(t,0)[0][0] for t in tlist])
H_list = H

psi0 = basis(2,0)
basisStates = [basis(2,i) for i in range(2)]
e_ops = [b*b.dag() for b in basisStates]

me =mesolve(H_list, psi0, tlist,  progress_bar = None, e_ops = e_ops)

plt.plot(me.expect[0])
plt.plot(me.expect[1])

plt.plot(h)
plt.show()

Which works really well and produces the desired output, a $\pi$ Rabi oscillation driven by a $\pi$ pulse. The red is the envelope of the Rabi hamiltonian, and the orange and blue are the (expected) populations of the upper and lower states of a two level system undergoing a pi pulse. This is the expected behaviour.

enter image description here

Except if I want the pulse one time unit later, when no evolution is driven. The code is exactly the same, except the hamiltonian function is modified to the following

def H(t, args):
    H = Qobj([[0.,1],[1,0.0]])
    envelope = pulse(t, 3,3+np.pi/2,0.2)
    H*=envelope
    return H

which produces

enter image description here

Or if I want the edges to be sharper, when again, no evolution occurs

def H(t, args):
    H = Qobj([[0.,1],[1,0.0]])
    envelope = pulse(t, 2,2+np.pi/2,0.01)
    H*=envelope
    return H

strong text

What element of these time dependent Hamiltonians functions is causing mesolve to correctly calculate the dynamics in some cases, and completely fail in others?

The difference between the Hamiltonians in the above plots seems minor and arbitrary to me, and the rest of the code is identical between plots.

EDIT

After Paul Nation's answer, here is updated, working code.

import numpy as np
from qutip import * 
from scipy.special import erf
import matplotlib.pyplot as plt

def pulseEdge(x, x0, sigma):
    pulseEdge = ((erf((x-x0)/sigma)) +1.)/2.
    return pulseEdge

def pulse(t, tOn, tOff, sigma):
    return pulseEdge(t, tOn, sigma)* (1.-pulseEdge(t, tOff, sigma))

def H(t, args):
    H = Qobj([[0.,1],[1,0.0]])
    envelope = pulse(t, 2,2+np.pi/2,0.01)
    H*=envelope
    return H

tlist = np.linspace(0,10,1000)
h = np.array([H(t,0)[0][0] for t in tlist])
H_list = H

psi0 = basis(2,0)
basisStates = [basis(2,i) for i in range(2)]
e_ops = [b*b.dag() for b in basisStates]

options = Options(max_step = 1)
me =mesolve(H_list, psi0, tlist,  progress_bar = None, e_ops = e_ops, options = options)

plt.plot(me.expect[0])
plt.plot(me.expect[1])

plt.plot(h)
plt.show()

enter image description here

$\endgroup$

1 Answer 1

1
$\begingroup$

So, this can happen when the ode solvers automatic step size solver thinks that nothing is going on and takes a big step over the pulse. There is a max_step option that you can set. In practice, setting this to something like half the width of the smallest pulse in a pulse sequence works quite well. Then the solver sees the pulses and adjusts accordingly.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.