J.R. Johansson and P.D. Nation
For more information about QuTiP see http://qutip.org
In [1]:
%pylab inline
In [2]:
from qutip import *
import time
In [3]:
def hamiltonian_t(t, args):
""" evaluate the hamiltonian at time t. """
H0 = args[0]
H1 = args[1]
return H0 + t * H1
In [4]:
def qubit_integrate(delta, eps0, A, gamma1, gamma2, psi0, tlist):
# Hamiltonian
sx = sigmax()
sz = sigmaz()
sm = destroy(2)
H0 = - delta/2.0 * sx - eps0/2.0 * sz
H1 = - A/2.0 * sz
# collapse operators
c_op_list = []
n_th = 0.0 # zero temperature
# relaxation
rate = gamma1 * (1 + n_th)
if rate > 0.0:
c_op_list.append(sqrt(rate) * sm)
# excitation
rate = gamma1 * n_th
if rate > 0.0:
c_op_list.append(sqrt(rate) * sm.dag())
# dephasing
rate = gamma2
if rate > 0.0:
c_op_list.append(sqrt(rate) * sz)
# evolve and calculate expectation values
# method 1: function callback which returns the time-depdent qobj
#H_args = (H0, H1)
#output = mesolve(hamiltonian_t, psi0, tlist, c_op_list, [sm.dag() * sm], H_args)
# method 2: a function callback that returns the coefficient for a qobj
#H = [H0, [H1, lambda x,y: x]]
#output = mesolve(H, psi0, tlist, c_op_list, [sm.dag() * sm], {})
# method 3: a string that defines the coefficient. The solver generates
# and compiles C code using cython. This method is usually the fastest
# for large systems or long time evolutions, but there is fixed-time
# overhead that makes it inefficient for small and short-time evolutions.
H = [H0, [H1, 't']]
output = mesolve(H, psi0, tlist, c_op_list, [sm.dag() * sm], {})
return output.expect[0]
In [5]:
#
# set up the calculation
#
delta = 0.5 * 2 * pi # qubit sigma_x coefficient
eps0 = 0.0 * 2 * pi # qubit sigma_z coefficient
A = 2.0 * 2 * pi # sweep rate
gamma1 = 0.0 # relaxation rate
gamma2 = 0.0 # dephasing rate
psi0 = basis(2,0) # initial state
tlist = linspace(-20.0, 20.0, 5000)
In [6]:
start_time = time.time()
p_ex = qubit_integrate(delta, eps0, A, gamma1, gamma2, psi0, tlist)
print 'time elapsed = ' + str(time.time() - start_time)
In [7]:
figure(figsize=(12,8))
plot(tlist, real(p_ex), 'b', tlist, real(1-p_ex), 'r')
plot(tlist, 1 - exp( - pi * delta **2 / (2 * A)) * ones(shape(tlist)), 'k')
xlabel('Time')
ylabel('Occupation probability')
title('Landau-Zener transition')
legend(("Excited state", "Ground state", "Landau-Zener formula"), loc=0);
In [8]:
def hamiltonian_t(t, args):
""" evaluate the hamiltonian at time t. """
H0 = args[0]
H1 = args[1]
w = args[2]
return H0 + cos(w * t) * H1
In [9]:
def qubit_integrate(delta, eps0, A, omega, gamma1, gamma2, psi0, tlist, option):
# Hamiltonian
sx = sigmax()
sz = sigmaz()
sm = destroy(2)
H0 = - delta/2.0 * sx - eps0/2.0 * sz
H1 = - A/2.0 * sz
H_args = (H0, H1, omega)
# collapse operators
c_op_list = []
n_th = 0.0 # zero temperature
# relaxation
rate = gamma1 * (1 + n_th)
if rate > 0.0:
c_op_list.append(sqrt(rate) * sm)
# excitation
rate = gamma1 * n_th
if rate > 0.0:
c_op_list.append(sqrt(rate) * sm.dag())
# dephasing
rate = gamma2
if rate > 0.0:
c_op_list.append(sqrt(rate) * sz)
if option == "dynamics":
# evolve and calculate expectation values
output = mesolve(hamiltonian_t, psi0, tlist, c_op_list, [sm.dag() * sm], H_args)
return output.expect[0]
else: # option = steadystate
# find the propagator for one driving period
T = 2*pi / omega
U = propagator(hamiltonian_t, T, c_op_list, H_args, options=Odeoptions(nsteps=5000))
# find the steady state of successive application of the propagator
rho_ss = propagator_steadystate(U)
return real(expect(sm.dag() * sm, rho_ss))
In [10]:
#
# set up the calculation: a strongly driven two-level system
# (repeated LZ transitions)
#
delta = 0.05 * 2 * pi # qubit sigma_x coefficient
eps0 = 0.0 * 2 * pi # qubit sigma_z coefficient
A = 2.0 * 2 * pi # sweep rate
gamma1 = 0.0001 # relaxation rate
gamma2 = 0.005 # dephasing rate
psi0 = basis(2,0) # initial state
omega = 0.05 * 2 * pi # driving frequency
T = (2*pi)/omega # driving period
tlist = linspace(0.0, 3 * T, 1500)
In [11]:
start_time = time.time()
p_ex = qubit_integrate(delta, eps0, A, omega, gamma1, gamma2, psi0, tlist, "dynamics")
print 'dynamics: time elapsed = ' + str(time.time() - start_time)
In [12]:
start_time = time.time()
p_ex_ss = qubit_integrate(delta, eps0, A, omega, gamma1, gamma2, psi0, tlist, "steadystate")
print 'steady state: time elapsed = ' + str(time.time() - start_time)
In [13]:
figure(figsize=(12,8))
subplot(211)
plot(tlist, real(p_ex), 'b', tlist, real(1-p_ex), 'r', tlist, ones(shape(tlist)) * p_ex_ss, 'k')
xlabel('Time')
ylabel('Probability')
title('Repeated Landau-Zener-like transitions')
legend(("Excited state", "Ground state", "Excited steady state"), loc=0)
subplot(212)
plot(tlist, -delta/2.0 * ones(shape(tlist)), 'r')
plot(tlist, -(eps0/2.0 + A/2.0 * cos(omega * tlist)), 'b')
legend(("sx coeff", "sz coeff"))
xlabel('Time')
ylabel('Coefficients in the Hamiltonian')
show()
In [14]:
start_time = time.time()
A_vec = 2 * pi * linspace(0.0, 5.0, 100)
p_ex_ss_vec = zeros(len(A_vec))
idx = 0
start_time = time.time()
for A in A_vec:
p_ex_ss_vec[idx] = qubit_integrate(delta, eps0, A, omega, gamma1, gamma2, psi0, tlist, "steadystate")
idx += 1
print 'time elapsed = ' + str(time.time() - start_time)
In [15]:
figure()
plot(A_vec/(2*pi), p_ex_ss_vec, 'b.-')
title("Steady state of repeated LZ transitions")
xlabel("driving amplitude")
ylabel("Occupation probability")
show()
Find the steady state of a strongly driven qubit as a function of driving amplitude and qubit bias.
Note: This calculation can takes a long time.
In [16]:
def hamiltonian_t(t, args):
#
# evaluate the hamiltonian at time t.
#
H0 = args[0]
H1 = args[1]
w = args[2]
return H0 + H1 * sin(w * t)
In [17]:
def sd_qubit_integrate(delta, eps0_vec, A_vec, w, gamma1, gamma2):
# Hamiltonian
sx = sigmax()
sz = sigmaz()
sm = destroy(2)
# collapse operators
c_op_list = []
n_th = 0.0 # zero temperature
# relaxation
rate = gamma1 * (1 + n_th)
if rate > 0.0:
c_op_list.append(sqrt(rate) * sm)
# excitation
rate = gamma1 * n_th
if rate > 0.0:
c_op_list.append(sqrt(rate) * sm.dag())
# dephasing
rate = gamma2
if rate > 0.0:
c_op_list.append(sqrt(rate) * sz)
N = len(A_vec)
M = len(eps0_vec)
p_ex = zeros([N, M]) #, dtype=complex)
T = 2 * pi / w
sn = sm.dag() * sm
# sweep over the driving amplitude and bias point, find the steady state
# for each point and store in a matrix
for n in range(0, N):
for m in range(0, M):
H0 = - delta/2.0 * sx - eps0_vec[m]/2.0 * sz
H1 = - A_vec[n] * sx
H_args = (H0, H1, w)
# find the propagator for one period of the time-dependent
# hamiltonian
U = propagator(hamiltonian_t, T, c_op_list, H_args)
# find the steady state of the driven system
rho_ss = propagator_steadystate(U)
p_ex[n, m] = real(expect(sn, rho_ss))
return p_ex
In [18]:
#
# set up the parameters
#
delta = 0.2 * 2 * pi # qubit sigma_x coefficient
w = 1.0 * 2 * pi # qubit sigma_z coefficient
A_vec = linspace(0.0, 4.0, 100) * 2 * pi # driving amplitude
eps0_vec = linspace(0.0, 4.0, 100) * 2 * pi # qubit sigma-z bias point
gamma1 = 0.05 # relaxation rate
gamma2 = 0.0 # dephasing rate
In [19]:
start_time = time.time()
p_ex = sd_qubit_integrate(delta, eps0_vec, A_vec, w, gamma1, gamma2)
print 'time elapsed = ' + str(time.time() - start_time)
In [22]:
fig, ax = plt.subplots(figsize=(10,10))
p = ax.pcolor(A_vec, eps0_vec, real(p_ex), edgecolors='none')
p.set_cmap('RdYlBu_r')
ax.set_ylabel(r'$A/\omega$', fontsize=20)
ax.set_xlabel(r'$\epsilon_0/\omega$', fontsize=20)
ax.axis('tight')
ax.set_title('Excitation probabilty of qubit, in steady state', fontsize=16);
In [21]:
from qutip.ipynbtools import version_table
version_table()
Out[21]: