# QuTiP development notebook for floquet and floquet-markov evolution¶

Copyright (C) 2011 and later, Paul D. Nation & Robert J. Johansson

In [1]:
%pylab inline

In [2]:
from qutip import *

In [3]:
import time


### Functions for solving and visualizing the dynamics of driven quantum systems with different solvers¶

In [4]:
def solve_dynamics(H, args, psi0, tlist, T, c_ops, e_ops, fc_ops, J_cbs):
"""
Solve the dynamics for a quantum system with a bunch of different solvers
and return the results. The idea is that the results should be quite similar,
even though some differences can be expected due to the different situation
(approximations) the different solvers are applicable for. Also, some solvers
here solve for unitary evolution and some the dissipative evolution, so
depending on the values c_ops and fc_ops the dynamics can be very different.
"""

sol_list = []

#
# Solve with unitary floquet solver by going through all steps manually
#
sol = Odedata()
sol.times = tlist
sol.solver = "Floquet WF"
sol.expect = zeros((len(e_ops), len(tlist)))

f_modes_0, f_energies = floquet_modes(H, T, args)
f_coeff = floquet_state_decomposition(f_modes_0, f_energies, psi0)
f_modes_table_t = floquet_modes_table(f_modes_0, f_energies, tlist, H, T, args)

for n, t in enumerate(tlist):
f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
psi_t = floquet_wavefunction(f_modes_t, f_energies, f_coeff, t)
for e_idx, e in enumerate(e_ops):
sol.expect[e_idx, n] = expect(e, psi_t)
sol_list.append(sol)

#
# Use fsesolve to do the same as above with a single function call
#
sol = fsesolve(H, psi0, tlist, e_ops, T, args)
sol_list.append(sol)

#
# Solve the floquet-markov master equation by requesting density matrices,
# and then loop through the states and calculate the expectation values.
# Do this with and without asking the solver to transform states back to the
# computational basis.
#
sol = fmmesolve(H, psi0, tlist, fc_ops, [], J_cbs, T, args, floquet_basis=True)
sol.solver = "Floquet ME states #1"
sol.expect = zeros((len(e_ops), len(tlist)))
for n, t in enumerate(sol.times):
f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
for e_idx, e in enumerate(e_ops):
sol.expect[e_idx, n] = expect(e, sol.states[n].transform(f_modes_t, False))
sol_list.append(sol)

sol = fmmesolve(H, psi0, tlist, fc_ops, [], J_cbs, T, args, floquet_basis=False)
sol.solver = "Floquet ME states #2"
sol.expect = zeros((len(e_ops), len(tlist)))
for n, t in enumerate(sol.times):
for e_idx, e in enumerate(e_ops):
sol.expect[e_idx, n] = expect(e, sol.states[n])
sol_list.append(sol)

#
# Solve the floquet-markov master equation by requesting expectation values
#
sol = fmmesolve(H, psi0, tlist, fc_ops, e_ops, J_cbs, T, args)
sol_list.append(sol)

#
# Solve with mesolve
#
sol = mesolve(H, psi0, tlist, c_ops, e_ops, args)
sol_list.append(sol)

return sol_list

In [5]:
def visualize_solutions(sol_list, split=False):

if split:
fig, axes = subplots(len(sol_list), 1, figsize=(12, 3 * len(sol_list)), squeeze=False)
else:
fig, axes = subplots(1, 1, figsize=(12, 6), squeeze=False)

for idx, sol in enumerate(sol_list):
if not split:
idx = 0

for e_idx, e in enumerate(sol.expect):
axes[idx,0].plot(sol.times, real(e), label="%s e_%d" % (sol.solver, e_idx))

axes[idx,0].set_xlabel('Time')
axes[idx,0].set_ylabel('Occupation probability')
axes[idx,0].legend(loc=0)


### Single atom rabi oscillations with classical drive¶

In [6]:
delta = 0.45 * 2*pi
eps0  = 1.0 * 2*pi
A     = 0.3 * 2*pi
omega = sqrt(delta**2 + eps0**2)
T      = (2*pi)/omega
tlist  = linspace(0.0, 10 * T, 501)
psi0   = rand_ket(2)

H0 = - eps0/2.0 * sigmaz() - delta/2.0 * sigmax()
H1 = A/2.0 * sigmax()
args = {'w': omega}
H = [H0, [H1, lambda t,args: sin(args['w'] * t)]]

gamma1 = 0.005
def noise_spectrum(omega):
return 0.25 * gamma1 * omega /(2*pi)

c_ops = [sqrt(gamma1) * sigmam()]
e_ops = [num(2)]
fc_ops = [sigmax()]
J_cbs = [noise_spectrum]

In [7]:
t = time.time()
sol_list = solve_dynamics(H, args, psi0, tlist, T, c_ops, e_ops, fc_ops, J_cbs)
print "elapsed =", time.time() - t

-c:50: ComplexWarning: Casting complex values to real discards the imaginary part
-c:58: ComplexWarning: Casting complex values to real discards the imaginary part

elapsed = 69.7039148808

In [8]:
visualize_solutions(sol_list)


#### Same thing with str-list time-dependent format, and stronger dissipation¶

In [9]:
H = [H0, [H1, "sin(w * t)"]]

gamma1 = 0.075
def noise_spectrum(omega):
return 0.25 * gamma1 * omega /(2*pi)

J_cbs = [noise_spectrum]
c_ops = [sqrt(gamma1) * sigmam()]

In [10]:
t = time.time()
sol_list = solve_dynamics(H, args, psi0, tlist, T, c_ops, e_ops, fc_ops, J_cbs)
print "elapsed =", time.time() - t

elapsed = 72.7279648781

In [11]:
visualize_solutions(sol_list)


### JC off resonant with driving¶

In [12]:
delta = 0.0 * 2*pi
eps0 = 1.0 * 2*pi
A = 0.4 * 2*pi
w0 = 2.0 * 2 * pi
g = 0.15 * 2 * pi

omega = w0 - sqrt(delta**2 + eps0**2)
T      = (2*pi)/omega
tlist  = linspace(0.0, 10 * T, 1001)
N = 2
a  = tensor(destroy(N), qeye(2))
sx = tensor(qeye(N), sigmax())
sz = tensor(qeye(N), sigmaz())
sn = tensor(qeye(N), num(2))
sm = tensor(qeye(N), destroy(2))
n = a.dag() * a

psi0 = tensor(basis(N, 1), basis(2,0))

H0 = w0 * a.dag() * a - eps0/2.0 * sz - delta/2.0 * sx + g * (a + a.dag()) * sx
H1 = A/2.0 * sz
args = {'w': omega}
H = [H0, [H1, lambda t,args: sin(args['w'] * t)]]

gamma1 = 0.001
def noise_spectrum(omega):
return 0.25 * gamma1 * omega /(2*pi)

c_ops = [sqrt(gamma1) * sm]
e_ops = [sn, n]
fc_ops = [sx]
J_cbs = [noise_spectrum]

In [13]:
t = time.time()
sol_list = solve_dynamics(H, args, psi0, tlist, T, c_ops, e_ops, fc_ops, J_cbs)
print "elapsed =", time.time() - t

elapsed = 267.251783848

In [14]:
visualize_solutions(sol_list)


### Software versions¶

In [15]:
from qutip.ipynbtools import version_table
version_table()

Out[15]:
SoftwareVersion
Cython0.16
SciPy0.10.1
QuTiP2.2.0.dev-02eb062
Python2.7.3 (default, Sep 26 2012, 21:51:14) [GCC 4.7.2]
IPython0.13
OSposix [linux2]
Numpy1.7.0.dev-3f45eaa
matplotlib1.3.x
Mon Feb 25 13:51:33 2013