J.R. Johansson and P.D. Nation
For more information about QuTiP see http://qutip.org
In [ ]:
%pylab inline
In [ ]:
from qutip import *
import time
Find the floquet modes and quasi energies for a driven system and plot the floquet states/quasienergies for one period of the driving.
In [ ]:
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 [ ]:
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 = zeros((len(A_vec), 2))
#f_gnd_prob = zeros((len(A_vec), 2))
quasi_energies = zeros((len(eps0_vec), 2))
f_gnd_prob = zeros((len(eps0_vec), 2))
wf_gnd_prob = 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 [ ]:
delta = 0.2 * 2 * pi # qubit sigma_x coefficient
eps0 = 0.5 * 2 * pi # qubit sigma_z coefficient
gamma1 = 0.0 # relaxation rate
gamma2 = 0.0 # dephasing rate
A = 2.0 * 2 * pi
psi0 = basis(2,0) # initial state
omega = 1.0 * 2 * pi # driving frequency
T = (2*pi)/omega # driving period
param = linspace(-5.0, 5.0, 200) * 2 * pi
eps0 = param
In [ ]:
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))
In [ ]:
figure(1)
plot(param, real(q_energies[:,0]) / delta, 'b', param, real(q_energies[:,1]) / delta, 'r')
xlabel('A or e')
ylabel('Quasienergy')
title('Floquet quasienergies');
In [ ]:
figure(2)
plot(param, real(f_gnd_prob[:,0]), 'b', param, real(f_gnd_prob[:,1]), 'r')
xlabel('A or e')
ylabel('Occ. prob.')
title('Floquet modes excitation probability');
In [ ]:
figure(3)
plot(param, real(wf_gnd_prob[:,0]), 'b', param, real(wf_gnd_prob[:,1]), 'r')
xlabel('A or e')
ylabel('Occ. prob.')
title('Floquet states excitation probability');
In [ ]:
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 [ ]:
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*pi / omega
f_modes_0,f_energies = floquet_modes(H, T, H_args)
p_ex_0 = zeros(shape(tlist))
p_ex_1 = zeros(shape(tlist))
p_00 = zeros(shape(tlist), dtype=complex)
p_01 = zeros(shape(tlist), dtype=complex)
p_10 = zeros(shape(tlist), dtype=complex)
p_11 = zeros(shape(tlist), dtype=complex)
e_0 = zeros(shape(tlist))
e_1 = 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(real(evals))
e_1[idx] = max(real(evals))
return p_ex_0, p_ex_1, e_0, e_1, f_energies, p_00, p_01, p_10, p_11
In [ ]:
delta = 0.2 * 2 * pi # qubit sigma_x coefficient
eps0 = 1.0 * 2 * pi # qubit sigma_z coefficient
A = 2.5 * 2 * pi # sweep rate
psi0 = basis(2,0) # initial state
omega = 1.0 * 2 * pi # driving frequency
T = (2*pi)/omega # driving period
tlist = linspace(0.0, T, 101)
In [ ]:
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))
In [ ]:
figure(figsize=[8,10])
subplot(2,1,1)
plot(tlist, real(p_ex_0), 'b', tlist, real(p_ex_1), 'r')
xlabel('Time ($T$)')
ylabel('Excitation probabilities')
title('Floquet modes')
legend(("Floquet mode 1", "Floquet mode 2"))
subplot(2,1,2)
plot(tlist, real(e_0), 'c', tlist, real(e_1), 'm')
plot(tlist, ones(shape(tlist)) * f_e[0], 'b', tlist, ones(shape(tlist)) * f_e[1], 'r')
xlabel('Time ($T$)')
ylabel('Energy [GHz]')
title('Eigen- and quasi-energies')
legend(("Ground state", "Excited state", "Quasienergy 1", "Quasienergy 2"));
In [ ]:
figure(figsize=[8,12])
subplot(3,1,1)
plot(tlist, real(p_00), 'b', tlist, real(p_01), 'r')
plot(tlist, real(p_10), 'c', tlist, real(p_11), 'm')
xlabel('Time ($T$)')
ylabel('real')
title('Floquet modes')
legend(("FM1 - gnd", "FM1 - exc", "FM2 - gnd", "FM2 - exc"))
subplot(3,1,2)
plot(tlist, imag(p_00), 'b', tlist, imag(p_01), 'r')
plot(tlist, imag(p_10), 'c', tlist, imag(p_11), 'm')
xlabel('Time ($T$)')
ylabel('imag')
legend(("FM1 - gnd", "FM1 - exc", "FM2 - gnd", "FM2 - exc"))
subplot(3,1,3)
plot(tlist, abs(p_00), 'b', tlist, abs(p_01), 'r.')
plot(tlist, abs(p_10), 'c', tlist, abs(p_11), 'm.')
xlabel('Time ($T$)')
ylabel('abs')
legend(("FM1 - gnd", "FM1 - exc", "FM2 - gnd", "FM2 - exc"));
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 [ ]:
def hamiltonian_t(t, args):
""" evaluate the hamiltonian at time t. """
H0 = args[0]
H1 = args[1]
w = args[2]
return H0 + sin(w * t) * H1
def H1coeff_t(t, args):
w = args['w']
return sin(w * t)
In [ ]:
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 = zeros(N)
tlist = []
# calculate the wavefunctions at times t=nT, for integer n, by using
# the floquet modes and quasienergies
for n in 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 [ ]:
delta = 0.2 * 2 * pi # qubit sigma_x coefficient
eps0 = 0.1 * 2 * pi # qubit sigma_z coefficient
A = 1.0 * 2 * pi # driving amplitude
psi0 = basis(2,0) # initial state
omega = 1.0 * 2 * pi # driving frequency
T = (2*pi)/omega # driving period
tlist = linspace(0.0, 25 * T, 500)
In [ ]:
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))
In [ ]:
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))
In [ ]:
figure(figsize=[6,4])
plot(tlist1, real(p_ex), 'b')
plot(tlist1, real(1-p_ex), 'r')
plot(tlist2, real(f_ex), 'bo', linewidth=2.0)
plot(tlist2, real(1-f_ex), 'ro', linewidth=2.0)
xlabel('Time')
ylabel('Probability')
title('Stroboscopic time-evolution with Floquet states')
legend(("ode $P_1$", "ode $P_0$", "Floquet $P_1$", "Floquet $P_0$"));
In [ ]:
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*pi)
def H1_coeff_t(t, args):
return sin(args['w'] * t)
In [ ]:
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 = [sqrt(gamma1) * sx, 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 = [sqrt(gamma1) * sx, 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
#
import qutip.odeconfig
qutip.odeconfig.tdfunc = None # better way of reseting this?!?
start_time = time.time()
T = 2*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 * pi * 1e-9
f_modes_table_t = floquet_modes_table(f_modes_0, f_energies, 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)
#expt_list4 = fmmesolve(R, f_modes_0, psi0, tlist, [sm.dag() * sm], opt=None) # note: in floquet basis...
rho_list = fmmesolve(R, f_modes_0, psi0, tlist, [])
expt_list3 = zeros(shape(expt_list2), dtype=complex)
expt_list4 = 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, False))
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, False) # in computational basis
expt_list5 = ones(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 [ ]:
delta = 0.0 * 2 * pi # qubit sigma_x coefficient
eps0 = 1.0 * 2 * pi # qubit sigma_z coefficient
A = 0.05 * 2 * pi # driving amplitude (reduce to make the RWA more accurate)
w = 1.0 * 2 * pi # driving frequency
psi0 = (0.3*basis(2,0)+0.7*basis(2,1)).unit() # initial state
tlist = linspace(0, 30.0, 500)
In [ ]:
p_ex1, p_ex2, p_ex3, p_ex4, p_ex5 = qubit_integrate(delta, eps0, A, w, gamma1, gamma2, psi0, tlist)
In [ ]:
figure()
plot(tlist, real(p_ex1), 'b', tlist, real(p_ex2), 'g-') # lindblad
plot(tlist, real(p_ex3), 'r', tlist, real(p_ex4), 'm-', tlist, real(p_ex5), 'c-') # floquet markov
xlabel('Time')
ylabel('Occupation probability')
title('Comparison between time-dependent master equations')
legend(("TD Hamiltonian, Std ME", "RWA Hamiltonian, Std ME",
"Unitary Floquet evol.", "Floquet-Markov ME", "F-M ME steady state"))
In [ ]:
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 + cos(w * t) * H1
In [ ]:
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*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 * 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
rho_ss = floquet_master_equation_steadystate(H0, 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 = zeros(shape(tlist))
p_ex_1 = zeros(shape(tlist))
e_0 = zeros(shape(tlist))
e_1 = 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(real(evals))
e_1[idx] = max(real(evals))
return p_ex_0, p_ex_1, e_0, e_1, f_energies
In [ ]:
delta = 0.2 * 2 * pi # qubit sigma_x coefficient
eps0 = 1.0 * 2 * pi # qubit sigma_z coefficient
A = 2.5 * 2 * pi # sweep rate
psi0 = basis(2,0) # initial state
omega = 1.0 * 2 * pi # driving frequency
T = (2*pi)/omega # driving period
tlist = linspace(0.0, 1 * T, 101)
In [ ]:
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))
In [ ]:
figure(figsize=[8,10])
subplot(2,1,1)
plot(tlist, real(p_ex_0), 'b', tlist, real(p_ex_1), 'r')
xlabel('Time ($T$)')
ylabel('Excitation probabilities')
title('Floquet modes')
legend(("Floquet mode 1", "Floquet mode 2"))
subplot(2,1,2)
plot(tlist, real(e_0), 'c', tlist, real(e_1), 'm')
plot(tlist, ones(shape(tlist)) * f_e[0], 'b', tlist, ones(shape(tlist)) * f_e[1], 'r')
xlabel('Time ($T$)')
ylabel('Energy [GHz]')
title('Eigen- and quasi-energies')
legend(("Ground state", "Excited state", "Quasienergy 1", "Quasienergy 2"));
In [ ]:
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 + sin(w * t) * H1
def H1_coeff_t(t, args):
return sin(args['w'] * t)
In [ ]:
def qubit_integrate(delta, eps0_vec, A, omega, gamma1, gamma2, psi0, T, option):
# Hamiltonian
sx = sigmax()
sz = sigmaz()
sm = destroy(2)
quasi_energies = zeros((len(eps0_vec), 2))
f_gnd_prob = zeros((len(eps0_vec), 2))
wf_gnd_prob = zeros((len(eps0_vec), 2))
ss_prob1 = zeros(len(eps0_vec))
ss_prob2 = 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 * 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
print "="*80
print "rho_ss_fb =\n", rho_ss_fb
print "rho_ss_cb =\n", rho_ss_cb
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 [ ]:
delta = 0.1 * 2 * pi # qubit sigma_x coefficient
eps0 = 1.0 * 2 * pi # qubit sigma_z coefficient
gamma1 = 0.0 # relaxation rate
gamma2 = 0.0 # dephasing rate
A = 2.0 * 2 * pi
psi0 = basis(2,0) # initial state
omega = sqrt(delta**2 + eps0**2) # driving frequency
T = (2*pi)/omega # driving period
param = linspace(-2.0, 2.0, 100) * 2 * pi
eps0 = param
In [ ]:
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))
In [ ]:
figure(1)
plot(param, real(q_energies[:,0]) / delta, 'b', param, real(q_energies[:,1]) / delta, 'r')
xlabel('A or e')
ylabel('Quasienergy')
title('Floquet quasienergies')
In [ ]:
figure(2)
plot(param, real(f_gnd_prob[:,0]), 'b', param, real(f_gnd_prob[:,1]), 'r')
xlabel('A or e')
ylabel('Occ. prob.')
title('Floquet modes excitation probability')
In [ ]:
figure(3)
plot(param, real(wf_gnd_prob[:,0]), 'b', param, real(wf_gnd_prob[:,1]), 'r')
xlabel('A or e')
ylabel('Occ. prob.')
title('Floquet states excitation probability')
In [ ]:
figure(4)
plot(param, real(ss_prob1), 'r')
plot(param, real(ss_prob2), 'b')
xlabel('A or e')
ylabel('Occ. prob. in steady state')
title('Steady state excitation probability')
In [ ]:
from qutip.ipynbtools import version_table
version_table()