# QuTiP example: Floquet formalism for dynamics and steadystates¶

J.R. Johansson and P.D. Nation

In [1]:
%matplotlib inline

In [2]:
import matplotlib.pyplot as plt

In [3]:
import numpy as np

In [4]:
from qutip import *
import time


## Floquet modes and quasi energies¶

Find the floquet modes and quasi energies for a driven system and plot the floquet states/quasienergies for one period of the driving.

In [5]:
def hamiltonian_t(t, args):
""" evaluate the hamiltonian at time t. """
H0 = args['H0']
H1 = args['H1']
w  = args['w']

return H0 + sin(w * t) * H1

def H1_coeff_t(t, args):
return sin(args['w'] * t)

In [6]:
def qubit_integrate(delta, eps0_vec, A, omega, gamma1, gamma2, psi0, T, option):

# 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)

quasi_energies = np.zeros((len(eps0_vec), 2))
f_gnd_prob     = np.zeros((len(eps0_vec), 2))
wf_gnd_prob    = np.zeros((len(eps0_vec), 2))

for idx, eps0 in enumerate(eps0_vec):

H0 = - delta/2.0 * sx - eps0/2.0 * sz
H1 = A/2.0 * sz

# H = H0 + H1 * sin(w * t) in the 'list-string' format
H = [H0, [H1, 'sin(w * t)']]
Hargs = {'w': omega}

# find the floquet modes
f_modes,f_energies = floquet_modes(H, T, Hargs)

quasi_energies[idx,:] = f_energies

f_gnd_prob[idx, 0] = expect(sm.dag() * sm, f_modes[0])
f_gnd_prob[idx, 1] = expect(sm.dag() * sm, f_modes[1])

f_states = floquet_states_t(f_modes, f_energies, 0, H, T, Hargs)

wf_gnd_prob[idx, 0] = expect(sm.dag() * sm, f_states[0])
wf_gnd_prob[idx, 1] = expect(sm.dag() * sm, f_states[1])

return quasi_energies, f_gnd_prob, wf_gnd_prob

In [7]:
delta = 0.2 * 2 * np.pi  # qubit sigma_x coefficient
eps0  = 0.5 * 2 * np.pi  # qubit sigma_z coefficient
gamma1 = 0.0        # relaxation rate
gamma2 = 0.0         # dephasing  rate
A      = 2.0 * 2 * np.pi
psi0   = basis(2,0)    # initial state
omega  = 1.0 * 2 * np.pi # driving frequency
T      = (2*np.pi)/omega  # driving period

param  = np.linspace(-5.0, 5.0, 200) * 2 * np.pi

eps0 = param

In [8]:
start_time = time.time()
q_energies, f_gnd_prob, wf_gnd_prob = qubit_integrate(delta, eps0, A, omega, gamma1, gamma2, psi0, T, "dynamics")
print('dynamics: time elapsed = ' + str(time.time() - start_time))

dynamics: time elapsed = 442.35331439971924

In [9]:
fig, ax = plt.subplots()
ax.plot(param, np.real(q_energies[:,0]) / delta, 'b',
param, np.real(q_energies[:,1]) / delta, 'r')
ax.set_xlabel('A or e')
ax.set_ylabel('Quasienergy')
ax.set_title('Floquet quasienergies');

In [10]:
fig, ax = plt.subplots()
ax.plot(param, np.real(f_gnd_prob[:,0]), 'b',
param, np.real(f_gnd_prob[:,1]), 'r')
ax.set_xlabel('A or e')
ax.set_ylabel('Occ. prob.')
ax.set_title('Floquet modes excitation probability');

In [11]:
fig, ax = plt.subplots()
ax.plot(param, np.real(wf_gnd_prob[:,0]), 'b',
param, np.real(wf_gnd_prob[:,1]), 'r')
ax.set_xlabel('A or e')
ax.set_ylabel('Occ. prob.')
ax.set_title('Floquet states excitation probability');


## Floquet modes¶

In [12]:
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

def H1coeff_t(t, args):
w  = args['w']
return sin(w * t)

In [13]:
def qubit_integrate(delta, eps0, A, omega, psi0, tlist):

# 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)
H_args = {'w': omega}
H = [H0, [H1, 'sin(w * t)']]
#H = [H0, [H1, H1coeff_t]]

# find the propagator for one driving period
T = 2*np.pi / omega

f_modes_0,f_energies = floquet_modes(H, T, H_args)

p_ex_0 = np.zeros(shape(tlist))
p_ex_1 = np.zeros(shape(tlist))

p_00 = np.zeros(shape(tlist), dtype=complex)
p_01 = np.zeros(shape(tlist), dtype=complex)
p_10 = np.zeros(shape(tlist), dtype=complex)
p_11 = np.zeros(shape(tlist), dtype=complex)

e_0 = np.zeros(shape(tlist))
e_1 = np.zeros(shape(tlist))

f_modes_table_t = floquet_modes_table(f_modes_0, f_energies, tlist, H, T, H_args)

for idx, t in enumerate(tlist):
#f_modes_t = floquet_modes_t(f_modes_0, f_energies, t, H, T, H_args)
f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)

p_ex_0[idx] = expect(sm.dag() * sm, f_modes_t[0])
p_ex_1[idx] = expect(sm.dag() * sm, f_modes_t[1])

p_00[idx] = f_modes_t[0].full()[0][0]
p_01[idx] = f_modes_t[0].full()[1][0]
p_10[idx] = f_modes_t[1].full()[0][0]
p_11[idx] = f_modes_t[1].full()[1][0]

#evals = hamiltonian_t(t, H_args).eigenenergies()
evals = qobj_list_evaluate(H, t, H_args).eigenenergies()
e_0[idx] = min(np.real(evals))
e_1[idx] = max(np.real(evals))

return p_ex_0, p_ex_1, e_0, e_1, f_energies, p_00, p_01, p_10, p_11

In [14]:
delta = 0.2 * 2 * np.pi  # qubit sigma_x coefficient
eps0  = 1.0 * 2 * np.pi  # qubit sigma_z coefficient
A     = 2.5 * 2 * np.pi  # sweep rate
psi0   = basis(2,0)   # initial state
omega  = 1.0 * 2 * np.pi # driving frequency
T      = (2*np.pi)/omega # driving period

tlist = np.linspace(0.0, T, 101)

In [15]:
start_time = time.time()
p_ex_0, p_ex_1, e_0, e_1, f_e, p_00, p_01, p_10, p_11 = qubit_integrate(delta, eps0, A, omega, psi0, tlist)
print('dynamics: time elapsed = ' + str(time.time() - start_time))

dynamics: time elapsed = 4.233538627624512

In [16]:
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=[8,10])

ax1.plot(tlist, np.real(p_ex_0), 'b', tlist, np.real(p_ex_1), 'r')
ax1.set_xlabel('Time ($T$)')
ax1.set_ylabel('Excitation probabilities')
ax1.set_title('Floquet modes')
ax1.legend(("Floquet mode 1", "Floquet mode 2"))

ax2.plot(tlist, np.real(e_0), 'c', tlist, np.real(e_1), 'm')
ax2.plot(tlist, np.ones(shape(tlist)) * f_e[0], 'b', tlist, np.ones(shape(tlist)) * f_e[1], 'r')
ax2.set_xlabel('Time ($T$)')
ax2.set_ylabel('Energy [GHz]')
ax2.set_title('Eigen- and quasi-energies')
ax2.legend(("Ground state", "Excited state", "Quasienergy 1", "Quasienergy 2"));

In [17]:
fig, axes = plt.subplots(3, 1, figsize=[8,12])

axes[0].plot(tlist, np.real(p_00), 'b', tlist, np.real(p_01), 'r')
axes[0].plot(tlist, np.real(p_10), 'c', tlist, np.real(p_11), 'm')
axes[0].set_xlabel('Time ($T$)')
axes[0].set_ylabel('real')
axes[0].set_title('Floquet modes')
axes[0].legend(("FM1 - gnd", "FM1 - exc", "FM2 - gnd", "FM2 - exc"))

axes[1].plot(tlist, np.imag(p_00), 'b', tlist, np.imag(p_01), 'r')
axes[1].plot(tlist, np.imag(p_10), 'c', tlist, np.imag(p_11), 'm')
axes[1].set_xlabel('Time ($T$)')
axes[1].set_ylabel('imag')
axes[1].legend(("FM1 - gnd", "FM1 - exc", "FM2 - gnd", "FM2 - exc"))

axes[2].plot(tlist, abs(p_00), 'b', tlist, abs(p_01), 'r.')
axes[2].plot(tlist, abs(p_10), 'c', tlist, abs(p_11), 'm.')
axes[2].set_xlabel('Time ($T$)')
axes[2].set_ylabel('abs')
axes[2].legend(("FM1 - gnd", "FM1 - exc", "FM2 - gnd", "FM2 - exc"));


## Floquet evolution¶

Find the floquet modes and quasi energies for a driven system and evolve the wavefunction "stroboscopically", i.e., only by evaluating at time mupliples of the driving period t = n * T for integer n.

The system is a strongly driven two-level atom.

In [18]:
def hamiltonian_t(t, args):
""" evaluate the hamiltonian at time t. """
H0 = args[0]
H1 = args[1]
w  = args[2]
return H0 + np.sin(w * t) * H1

def H1coeff_t(t, args):
w  = args['w']
return np.sin(w * t)

In [19]:
def qubit_integrate(delta, eps0, A, omega, psi0, tlist, T, 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)
H_args = {'w': omega}
H = [H0, [H1, 'sin(w * t)']]
#H = [H0, [H1, H1coeff_t]]

if option == "floquet":

# find the floquet modes for the time-dependent hamiltonian
f_modes_0,f_energies = floquet_modes(H, T, H_args)

# decompose the inital state in the floquet modes (=floquet states at t=0)
f_coeff = floquet_state_decomposition(f_modes_0, f_energies, psi0)

# only evaluate the wavefunction at multiples of the driving period
# i.e. a "stroboscopic" evolution
N = max(tlist)/T
p_ex = np.zeros(N)
tlist = []

# calculate the wavefunctions at times t=nT, for integer n, by using
# the floquet modes and quasienergies
for n in np.arange(N):
psi_t = floquet_wavefunction_t(f_modes_0, f_energies, f_coeff, n*T, H, T, H_args)
p_ex[n] = expect(sm.dag() * sm, psi_t)
tlist.append(n*T)

else:

# for reference: evolve and calculate expectation using the full ode solver
#H_args = (H0, H1, omega)
#expt_list = mesolve(hamiltonian_t, psi0, tlist, [], [sm.dag() * sm], H_args)
output = mesolve(H, psi0, tlist, [], [sm.dag() * sm], H_args)
p_ex = output.expect[0]

return tlist, p_ex

In [20]:
delta = 0.2 * 2 * np.pi  # qubit sigma_x coefficient
eps0  = 0.1 * 2 * np.pi  # qubit sigma_z coefficient
A     = 1.0 * 2 * np.pi  # driving amplitude
psi0   = basis(2,0)      # initial state
omega  = 1.0 * 2 * np.pi # driving frequency

T      = (2*np.pi)/omega # driving period
tlist  = np.linspace(0.0, 25 * T, 500)

In [21]:
start_time = time.time()
tlist1, p_ex = qubit_integrate(delta, eps0, A, omega, psi0, tlist, T, "dynamics")
print('dynamics: time elapsed = ' + str(time.time() - start_time))

dynamics: time elapsed = 3.316725730895996

In [22]:
start_time = time.time()
tlist2, f_ex = qubit_integrate(delta, eps0, A, omega, psi0, tlist, T, "floquet")
print('floquet: time elapsed = ' + str(time.time() - start_time))

floquet: time elapsed = 2.8527965545654297

In [23]:
fig, ax = plt.subplots(figsize=[12, 6])
ax.plot(tlist1, np.real(p_ex),   'b')
ax.plot(tlist1, np.real(1-p_ex), 'r')
ax.plot(tlist2, np.real(f_ex),   'bo', linewidth=2.0)
ax.plot(tlist2, np.real(1-f_ex), 'ro', linewidth=2.0)

ax.set_xlabel('Time')
ax.set_ylabel('Probability')
ax.set_title('Stroboscopic time-evolution with Floquet states')
ax.legend(("ode $P_1$", "ode $P_0$", "Floquet $P_1$", "Floquet $P_0$"));


## Floquet-Markov master equation: comparison with other solvers¶

In [24]:
gamma1 = 0.0015          # relaxation rate
gamma2 = 0.0           # dephasing  rate

def J_cb(omega):
""" Noise spectral density """
#print "evaluate J_cb for omega =", omega
return 0.5 * gamma1 * omega/(2*np.pi)

def H1_coeff_t(t, args):
return np.sin(args['w'] * t)

In [25]:
def qubit_integrate(delta, eps0, A, w, gamma1, gamma2, psi0, tlist):

# Hamiltonian
sx = sigmax()
sz = sigmaz()
sm = destroy(2)

H0 = - delta/2.0 * sx - eps0/2.0 * sz
H1 = - A * sx

args = {'w': w}
H = [H0, [H1, 'sin(w * t)']]

# --------------------------------------------------------------------------
# 1) time-dependent hamiltonian
#
#
c_op_list = [np.sqrt(gamma1) * sx, np.sqrt(gamma2) * sz]

start_time = time.time()
output = mesolve(H, psi0, tlist, c_op_list, [sm.dag() * sm], args=args)
expt_list1 = output.expect
print('Method 1: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 2) Constant hamiltonian
#
H_rwa = - delta/2.0 * sx - A * sx / 2
c_op_list = [np.sqrt(gamma1) * sx, np.sqrt(gamma2) * sz]

start_time = time.time()
output = mesolve(H_rwa, psi0, tlist, c_op_list, [sm.dag() * sm])
expt_list2 = output.expect
print('Method 2: time elapsed = ' + str(time.time() - start_time))

# --------------------------------------------------------------------------
# 3) Floquet unitary evolution
#
qutip.solver.config.reset()

start_time = time.time()

T = 2*np.pi / w
f_modes_0,f_energies = floquet_modes(H, T, args)

# decompose the initial state vector in terms of the floquet modes (basis
# transformation). used to calculate psi_t below.
f_coeff = floquet_state_decomposition(f_modes_0, f_energies, psi0)

# --------------------------------------------------------------------------
# 4) Floquet markov master equation dynamics
#
kmax = 1
temp = 25e-3
w_th = temp * (1.38e-23 / 6.626e-34) * 2 * np.pi * 1e-9

f_modes_table_t = floquet_modes_table(f_modes_0, f_energies, np.linspace(0, T, 500+1), H, T, args)

# calculate the rate-matrices for the floquet-markov master equation
Delta, X, Gamma, Amat = floquet_master_equation_rates(f_modes_0, f_energies, sx,
H, T, args, J_cb, w_th, kmax, f_modes_table_t)

# the floquet-markov master equation tensor
R = floquet_master_equation_tensor(Amat, f_energies)

rho_list = floquet_markov_mesolve(R, f_modes_0, psi0, tlist, []).states

expt_list3 = np.zeros(shape(expt_list2), dtype=complex)
expt_list4 = np.zeros(shape(expt_list2), dtype=complex)
for idx, t in enumerate(tlist):

# unitary floquet evolution
psi_t = floquet_wavefunction_t(f_modes_0, f_energies, f_coeff, t, H, T, args)
expt_list3[0][idx] = expect(sm.dag() * sm, psi_t)

# the rho_list returned by the floquet master equation is defined in the
# floquet basis, so to transform it back to the computational basis
# before we calculate expectation values.
#f_modes_t = floquet_modes_t(f_modes_0, f_energies, t, H, T, args)
f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
expt_list4[0][idx] = expect((sm.dag() * sm), rho_list[idx].transform(f_modes_t, True))

print('Method 3+4: time elapsed = ' + str(time.time() - start_time))

# calculate the steadystate density matrix according to the floquet-markov
# master equation
rho_ss_fb = floquet_master_equation_steadystate(H0, Amat) # in floquet basis
rho_ss_cb = rho_ss_fb.transform(f_modes_0, True)         # in computational basis
expt_list5 = np.ones(np.shape(expt_list2), dtype=complex) * expect(sm.dag() * sm, rho_ss_cb)

return expt_list1[0], expt_list2[0], expt_list3[0], expt_list4[0], expt_list5[0]

In [26]:
delta = 0.0 * 2 * np.pi   # qubit sigma_x coefficient
eps0  = 1.0 * 2 * np.pi   # qubit sigma_z coefficient
A     = 0.05 * 2 * np.pi  # driving amplitude (reduce to make the RWA more accurate)
w     = 1.0 * 2 * np.pi   # driving frequency
psi0 = (0.3*basis(2,0) + 0.7*basis(2,1)).unit()      # initial state

tlist = np.linspace(0, 30.0, 500)

In [27]:
p_ex1, p_ex2, p_ex3, p_ex4, p_ex5 = qubit_integrate(delta, eps0, A, w, gamma1, gamma2, psi0, tlist)

Method 1: time elapsed = 2.65116810798645
Method 2: time elapsed = 0.08068442344665527
Method 3+4: time elapsed = 1466.1850209236145

In [28]:
fig, ax = plt.subplots()
ax.plot(tlist, np.real(p_ex1), 'b', tlist, np.real(p_ex2), 'g-') # lindblad
ax.plot(tlist, np.real(p_ex3), 'r', tlist, np.real(p_ex4), 'm-', tlist, np.real(p_ex5), 'c-') # floquet markov
ax.set_xlabel('Time')
ax.set_ylabel('Occupation probability')
ax.set_title('Comparison between time-dependent master equations')
ax.legend(("TD Hamiltonian, Std ME", "RWA Hamiltonian, Std ME",
"Unitary Floquet evol.", "Floquet-Markov ME", "F-M ME steady state"));


## Floquet-Markov master equation¶

In [29]:
def J_cb(omega):
""" Noise spectral density """
return omega

def hamiltonian_t(t, args):
""" evaluate the hamiltonian at time t. """
H0 = args[0]
H1 = args[1]
w  = args[2]

return H0 + np.cos(w * t) * H1

In [30]:
def qubit_integrate(delta, eps0, A, omega, psi0, tlist):

# 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)
H_args = {'w': omega}
H = [H0, [H1, 'sin(w * t)']]

# find the propagator for one driving period
T = 2*np.pi / omega

f_modes_0,f_energies = floquet_modes(H, T, H_args)

c_op = sigmax()

kmax = 1

temp = 25e-3
w_th = temp * (1.38e-23 / 6.626e-34) * 2 * np.pi * 1e-9

Delta, X, Gamma, A = floquet_master_equation_rates(f_modes_0, f_energies, c_op, H, T, H_args, J_cb, w_th, kmax)

k_idx = 0
for k in range(-kmax,kmax+1, 1):
print("X[",k,"] =\n", X[:,:,k_idx])
k_idx += 1

k_idx = 0
for k in range(-kmax,kmax+1, 1):
print("Delta[",k,"] =\n", Delta[:,:,k_idx])
k_idx += 1

k_idx = 0
for k in range(-kmax,kmax+1, 1):
print("Gamma[",k,"] =\n", Gamma[:,:,k_idx])
k_idx += 1

print("A =\n", A)

R = floquet_master_equation_tensor(A, f_energies)

print("Floquet-Markov master equation tensor")

print("R =\n", R)

print("Floquet-Markov master equation steady state =\n", rho_ss)

p_ex_0 = np.zeros(shape(tlist))
p_ex_1 = np.zeros(shape(tlist))

e_0 = np.zeros(shape(tlist))
e_1 = np.zeros(shape(tlist))

f_modes_table_t = floquet_modes_table(f_modes_0, f_energies, tlist, H, T, H_args)
for idx, t in enumerate(tlist):
f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)

p_ex_0[idx] = expect(sm.dag() * sm, f_modes_t[0])
p_ex_1[idx] = expect(sm.dag() * sm, f_modes_t[1])

#evals = hamiltonian_t(t, H_args).eigenenergies()
evals = qobj_list_evaluate(H, t, H_args).eigenenergies()
e_0[idx] = min(np.real(evals))
e_1[idx] = max(np.real(evals))

return p_ex_0, p_ex_1, e_0, e_1, f_energies

In [31]:
delta = 0.2 * 2 * np.pi  # qubit sigma_x coefficient
eps0  = 1.0 * 2 * np.pi  # qubit sigma_z coefficient
A     = 2.5 * 2 * np.pi  # sweep rate
psi0   = basis(2,0)      # initial state
omega  = 1.0 * 2 * np.pi # driving frequency
T      = (2*np.pi)/omega # driving period

tlist = np.linspace(0.0, 1 * T, 101)

In [32]:
start_time = time.time()
p_ex_0, p_ex_1, e_0, e_1, f_e = qubit_integrate(delta, eps0, A, omega, psi0, tlist)
print('dynamics: time elapsed = ' + str(time.time() - start_time))

X[ -1 ] =
[[ -9.03801955e-07-0.19524686j   2.02821511e-01+0.31089132j]
[ -1.22745970e-03+0.00188092j   9.03801953e-07+0.19524686j]]
X[ 0 ] =
[[-0.48721741+0.j          0.26340704-0.17184237j]
[ 0.26340704+0.17184237j  0.48721741+0.j        ]]
X[ 1 ] =
[[ -9.03801955e-07+0.19524686j  -1.22745970e-03-0.00188092j]
[  2.02821511e-01-0.31089132j   9.03801953e-07-0.19524686j]]
Delta[ -1 ] =
[[ -6.28318531 -11.94580954]
[ -0.62056108  -6.28318531]]
Delta[ 0 ] =
[[ 0.         -5.66262423]
[ 5.66262423  0.        ]]
Delta[ 1 ] =
[[  6.28318531   0.62056108]
[ 11.94580954   6.28318531]]
Gamma[ -1 ] =
[[-0. -0.]
[-0. -0.]]
Gamma[ 0 ] =
[[ 0.         -0.        ]
[ 3.51925944  0.        ]]
Gamma[ 1 ] =
[[  1.50497000e+00   1.96690829e-05]
[  1.03422039e+01   1.50497000e+00]]
A =
[[  2.02170555   1.03320708]
[ 14.89465076   2.02170555]]
Floquet-Markov master equation tensor
R =
Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isherm = False
Qobj data =
[[-14.89465076+0.j           0.00000000+0.j           0.00000000+0.j
1.03320708+0.j        ]
[  0.00000000+0.j          -9.98563447-5.66262423j   0.00000000+0.j
0.00000000+0.j        ]
[  0.00000000+0.j           0.00000000+0.j          -9.98563447+5.66262423j
0.00000000+0.j        ]
[ 14.89465076+0.j           0.00000000+0.j           0.00000000+0.j
-1.03320708+0.j        ]]
Floquet-Markov master equation steady state =
Quantum object: dims = [[2], [2]], shape = [2, 2], type = oper, isherm = True
Qobj data =
[[ 0.06818126+0.j         -0.03313332+0.04199644j]
[-0.03313332-0.04199644j  0.93181874+0.j        ]]
dynamics: time elapsed = 8.292602300643921

In [33]:
fig, axes = plt.subplots(2, 1, figsize=[8,10])

axes[0].plot(tlist, np.real(p_ex_0), 'b', tlist, np.real(p_ex_1), 'r')
axes[0].set_xlabel('Time ($T$)')
axes[0].set_ylabel('Excitation probabilities')
axes[0].set_title('Floquet modes')
axes[0].legend(("Floquet mode 1", "Floquet mode 2"))

axes[1].plot(tlist, np.real(e_0), 'c', tlist, np.real(e_1), 'm')
axes[1].plot(tlist, np.ones_like(tlist) * f_e[0], 'b', tlist, np.ones_like(tlist) * f_e[1], 'r')
axes[1].set_xlabel('Time ($T$)')
axes[1].set_ylabel('Energy [GHz]')
axes[1].set_title('Eigen- and quasi-energies')
axes[1].legend(("Ground state", "Excited state", "Quasienergy 1", "Quasienergy 2"));


In [34]:
def J_cb(omega):
""" Noise spectral density """
return omega

def hamiltonian_t(t, args):
""" evaluate the hamiltonian at time t. """
H0 = args['H0']
H1 = args['H1']
w  = args['w']

return H0 + np.sin(w * t) * H1

def H1_coeff_t(t, args):
return np.sin(args['w'] * t)

In [35]:
def qubit_integrate(delta, eps0_vec, A, omega, gamma1, gamma2, psi0, T, option):

# Hamiltonian
sx = sigmax()
sz = sigmaz()
sm = destroy(2)

quasi_energies = np.zeros((len(eps0_vec), 2))
f_gnd_prob     = np.zeros((len(eps0_vec), 2))
wf_gnd_prob    = np.zeros((len(eps0_vec), 2))
ss_prob1       = np.zeros(len(eps0_vec))
ss_prob2       = np.zeros(len(eps0_vec))

Hargs = {'w': omega}

for idx, eps0 in enumerate(eps0_vec):

H0 = - delta/2.0 * sx - eps0/2.0 * sz
H1 = A/2.0 * sz
H = [H0, [H1, 'sin(w * t)']]

f_modes,f_energies = floquet_modes(H, T, Hargs)

quasi_energies[idx,:] = f_energies

f_gnd_prob[idx, 0] = expect(sm.dag() * sm, f_modes[0])
f_gnd_prob[idx, 1] = expect(sm.dag() * sm, f_modes[1])

f_states = floquet_states_t(f_modes, f_energies, 0, H, T, Hargs)

wf_gnd_prob[idx, 0] = expect(sm.dag() * sm, f_states[0])
wf_gnd_prob[idx, 1] = expect(sm.dag() * sm, f_states[1])

c_op = sigmax()
kmax = 5
temp = 0e-3
w_th = temp * (1.38e-23 / 6.626e-34) * 2 * np.pi * 1e-9
Delta, X, Gamma, Amat = floquet_master_equation_rates(f_modes, f_energies, c_op, H, T, Hargs, J_cb, w_th, kmax)

rho_ss_fb = floquet_master_equation_steadystate(H0, Amat) # floquet basis
rho_ss_cb = rho_ss_fb.transform(f_modes, True) #False           # computational basis

ss_prob1[idx] = expect(sm.dag() * sm, rho_ss_fb)
ss_prob2[idx] = expect(sm.dag() * sm, rho_ss_cb)

return quasi_energies, f_gnd_prob, wf_gnd_prob, ss_prob1, ss_prob2

In [36]:
delta = 0.1 * 2 * np.pi  # qubit sigma_x coefficient
eps0  = 1.0 * 2 * np.pi  # qubit sigma_z coefficient
gamma1 = 0.0             # relaxation rate
gamma2 = 0.0             # dephasing  rate
A      = 2.0 * 2 * np.pi
psi0   = basis(2,0)    # initial state
omega  = np.sqrt(delta**2 + eps0**2) # driving frequency
T      = (2*np.pi)/omega  # driving period

param  = np.linspace(-2.0, 2.0, 100) * 2 * np.pi

eps0 = param

In [37]:
start_time = time.time()
q_energies, f_gnd_prob, wf_gnd_prob, ss_prob1, ss_prob2 = qubit_integrate(delta, eps0, A, omega,
gamma1, gamma2, psi0, T, "dynamics")
print('dynamics: time elapsed = ' + str(time.time() - start_time))

dynamics: time elapsed = 1151.9405281543732

In [39]:
fig, ax = plt.subplots()
ax.plot(param, np.real(q_energies[:,0]) / delta, 'b',
param, np.real(q_energies[:,1]) / delta, 'r')
ax.set_xlabel('A or e')
ax.set_ylabel('Quasienergy')
ax.set_title('Floquet quasienergies')

Out[39]:
<matplotlib.text.Text at 0x7f57f755e3c8>
In [40]:
fig, ax = plt.subplots()
ax.plot(param, np.real(f_gnd_prob[:,0]), 'b',
param, np.real(f_gnd_prob[:,1]), 'r')
ax.set_xlabel('A or e')
ax.set_ylabel('Occ. prob.')
ax.set_title('Floquet modes excitation probability')

Out[40]:
<matplotlib.text.Text at 0x7f57f6ad5630>
In [41]:
fig, ax = plt.subplots()
ax.plot(param, np.real(wf_gnd_prob[:,0]), 'b',
param, np.real(wf_gnd_prob[:,1]), 'r')
ax.set_xlabel('A or e')
ax.set_ylabel('Occ. prob.')
ax.set_title('Floquet states excitation probability')

Out[41]:
<matplotlib.text.Text at 0x7f57bd767828>
In [42]:
fig, ax = plt.subplots()
ax.plot(param, np.real(ss_prob1), 'r')
ax.plot(param, np.real(ss_prob2), 'b')
ax.set_xlabel('A or e')

Out[42]:
<matplotlib.text.Text at 0x7f57d6ae0518>

## Versions¶

In [43]:
from qutip.ipynbtools import version_table

version_table()

Out[43]:
SoftwareVersion
Numpy1.9.1
IPython2.3.1
Cython0.21.2
Python3.4.0 (default, Apr 11 2014, 13:05:11) [GCC 4.8.2]
matplotlib1.4.2
QuTiP3.1.0
SciPy0.14.1
OSposix [linux]
Wed Jan 14 13:18:15 2015 JST