QuTiP example: Landau-Zener-Stuckelberg inteferometry

J.R. Johansson and P.D. Nation

For more information about QuTiP see http://qutip.org

In [1]:
%matplotlib inline
In [2]:
import matplotlib.pyplot as plt
In [3]:
import numpy as np
In [4]:
from qutip import *
from qutip.ui.progressbar import TextProgressBar as ProgressBar

Landau-Zener-Stuckelberg interferometry: Steady state of a strongly driven two-level system, using the one-period propagator.

In [5]:
# set up the parameters and start calculation
delta  = 1.0  * 2 * np.pi  # qubit sigma_x coefficient
w      = 2.0  * 2 * np.pi  # driving frequency
T      = 2 * np.pi / w     # driving period 
gamma1 = 0.00001        # relaxation rate
gamma2 = 0.005          # dephasing  rate

eps_list = np.linspace(-20.0, 20.0, 101) * 2 * np.pi
A_list   = np.linspace(  0.0, 20.0, 101) * 2 * np.pi

# pre-calculate the necessary operators
sx = sigmax(); sz = sigmaz(); sm = destroy(2); sn = num(2)

# collapse operators
c_op_list = [np.sqrt(gamma1) * sm, np.sqrt(gamma2) * sz]  # relaxation and dephasing
In [6]:
# ODE settings (for list-str format)
options = Options()
options.rhs_reuse = True
options.atol = 1e-6 # reduce accuracy to speed
options.rtol = 1e-5 # up the calculation a bit
In [7]:
# for function-callback style time-dependence
def hamiltonian_t(t, args):
    """ evaluate the hamiltonian at time t. """
    H0 = args[0]
    H1 = args[1]
    w  = args[2]
    return H0 + H1 * np.sin(w * t)
In [8]:
# perform the calculation for each combination of eps and A, store the result
# in a matrix
def calculate():

    p_mat = np.zeros((len(eps_list), len(A_list)))

    pbar = ProgressBar(len(eps_list))
    
    for m, eps in enumerate(eps_list):
        H0 = - delta/2.0 * sx - eps/2.0 * sz

        pbar.update(m)

        for n, A in enumerate(A_list):
            H1 = (A/2) * sz

            # function callback format
            #args = (H0, H1, w); H_td = hamiltonian_t

            # list-str format
            #args = {'w': w}; H_td = [H0, [H1, 'sin(w * t)']]

            # list-function format
            args = w; H_td = [H0, [H1, lambda t, w: np.sin(w * t)]]

            U = propagator(H_td, T, c_op_list, args, options)
            rho_ss = propagator_steadystate(U)

            p_mat[m,n] = np.real(expect(sn, rho_ss))

    return p_mat
In [9]:
p_mat = calculate()
10.9%. Run time: 2184.90s. Est. time left: 00:04:57:56
20.8%. Run time: 3892.79s. Est. time left: 00:04:07:09
30.7%. Run time: 5420.77s. Est. time left: 00:03:24:00
40.6%. Run time: 6764.94s. Est. time left: 00:02:44:59
50.5%. Run time: 8029.01s. Est. time left: 00:02:11:11
60.4%. Run time: 9306.51s. Est. time left: 00:01:41:42
70.3%. Run time: 10500.70s. Est. time left: 00:01:13:56
80.2%. Run time: 11884.90s. Est. time left: 00:00:48:54
90.1%. Run time: 13614.38s. Est. time left: 00:00:24:56
In [10]:
fig, ax = plt.subplots(figsize=(8, 8))

A_mat, eps_mat = np.meshgrid(A_list/(2*np.pi), eps_list/(2*np.pi))

ax.pcolor(eps_mat, A_mat, p_mat)
ax.set_xlabel(r'Bias point $\epsilon$')
ax.set_ylabel(r'Amplitude $A$')
ax.set_title("Steadystate excitation probability\n" +
             r'$H = -\frac{1}{2}\Delta\sigma_x -\frac{1}{2}\epsilon\sigma_z - \frac{1}{2}A\sin(\omega t)$' + "\n");

Versions

In [12]:
from qutip.ipynbtools import version_table

version_table()
Out[12]:
SoftwareVersion
SciPy0.14.1
Numpy1.9.1
Python3.4.0 (default, Apr 11 2014, 13:05:11) [GCC 4.8.2]
OSposix [linux]
matplotlib1.4.2
QuTiP3.1.0
IPython2.3.1
Cython0.21.2
Tue Jan 13 21:11:28 2015 JST