QuTiP development notebook for testing new TD function callback signature

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

Note: experimental and under development

In [1]:
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
In [2]:
from qutip import *
In [3]:
from qutip.ipynbtools import HTMLProgressBar
from qutip.gui.progressbar import TextProgressBar

Problem description

Consider the mean-field laser master equation on the form (Breuer and Petruccione)

$\displaystyle \dot\rho = -i [H, \rho]

  • 2\kappa \mathcal{D}[a]
  • W{21} \mathcal{D}[\sigma-]
  • W{12} \mathcal{D}[\sigma+]
  • g [{\rm Tr}[\sigma-\rho] a^\dagger - {\rm Tr}[\sigma+\rho] a, \rho]
  • g [{\rm Tr}[a^\dagger\rho] \sigma- - {\rm Tr}[a\rho] \sigma+, \rho] $

where the dissipator superoperator is

$\displaystyle \mathcal{D}[a] = a\rho a^\dagger - \frac{1}{2}a^\dagger a \rho - \frac{1}{2}\rho a^\dagger a,$

$W_{21}$ and $W_{12}$ are the atomic relaxation rate and pump rate, respectively, and ${\rm Tr}[A\rho]$ is the expectation value of the operator $A$ with respect to the density operator $\rho$.

The Hamiltonian is given by

$\displaystyle H = \omega a^\dagger a + \frac{1}{2}\omega\sigma_z$.

Except for the two last term, the above master equation is on standard Lindblad form and could, if $g = 0$, be solved with standard usage of the QuTiP solver mesolve.

We can write the master equation above on standard Lindblad form using an effective (nonhermitian?) Hamiltonian, corresponding to a nonlinear Schrodinger equation,

$\displaystyle \dot\rho = -i [H_{\rm eff}(\rho, t), \rho]

  • 2\kappa \mathcal{D}[a] + W{21} \mathcal{D}[\sigma-] + W{12} \mathcal{D}[\sigma+]$

where

$\displaystyle H_{\rm eff}(\rho, t) = \omega a^\dagger a + \frac{1}{2}\omega\sigma_z

  • ig\left({\rm Tr}[\sigma-\rho] a^\dagger - {\rm Tr}[\sigma+\rho] a\right)
  • ig\left({\rm Tr}[a^\dagger\rho] \sigma- - {\rm Tr}[a\rho] \sigma+\right) $

In QuTiP (development version required) we implement this effective Hamiltonian using the Python callback function format for time-dependent Hamiltonians, but to calculate $H_{\rm eff}(t)$ we need to have access to $\rho(t)$, and the standard QuTiP time-dependent function callback signature is

def h_t(t, args):
    ...
    return h_eff

However, the solver mesolve that calls the callback function has access to $\rho(t)$ (or $\left|\psi(t)\right>$ for unitary dynamics), so it could change the callback function signature to

def h_t(t, rho, args):
    ...
    return h_eff

To avoid breaking backwards compatibility in the API, the inclusion of rho in the function signature is optionally activated with a keyword argument (rhs_with_state=True) in Odeoptions.

QuTiP implementation

In [4]:
N = 10

w = 1.0 * 2 * pi
g = 0.1 * 2 * pi
kappa = 0.005
W21 = 0.01
W12 = 0.05

tlist = linspace(0, 250, 250)

# cavity operators
a  = tensor(destroy(N), identity(2))
ad = tensor(create(N), identity(2))

# atomic operators
sz = tensor(identity(N), sigmaz())
sx = tensor(identity(N), sigmax())
sm = tensor(identity(N), destroy(2))
sp = tensor(identity(N), create(2))
In [5]:
#psi0 = tensor(fock(N, 0), fock(2, 0)) # start with the vacuum + ground state
psi0 = tensor(coherent(N, 0.5), fock(2, 0)) # start a small coherent sate + ground state
In [6]:
theta = 0.0

H = w * ad * a - 0.5 * w * (cos(2*theta) * sz + sin(2*theta) * sx)

H
Out[6]:
\begin{equation}\text{Quantum object: dims = [[10, 2], [10, 2]], shape = [20, 20], type = oper, isHerm = True}\\[1em]\begin{pmatrix}-3.142 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 3.142 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 3.142 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 9.425 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 9.425 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 47.124 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 47.124 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 53.407 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 53.407 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 59.690\\\end{pmatrix}\end{equation}
In [7]:
e_ops = [ad * a, sm.dag() * sm, sm, sp]
In [8]:
c_ops = [sqrt(2*kappa) * a, sqrt(W21) * sm, sqrt(W12) * sp]

Solve without mean-field contribution

In [9]:
result = mesolve(H, psi0, tlist, c_ops, e_ops, options=Odeoptions(store_final_state=True)) 
In [10]:
fig, axes = subplots(2,1)

axes[0].plot(result.times, result.expect[0], label=r'$a^\dagger a$')
axes[0].plot(result.times, result.expect[1], label=r'$\sigma_+\sigma_-$')
axes[0].set_ylim(-0.1, 1.1)
axes[0].legend();

axes[1].plot(result.times, abs(result.expect[2]), label=r'$\sigma_-$')
axes[1].plot(result.times, abs(result.expect[3]), label=r'$\sigma_+$')
axes[1].set_ylim(-0.1, 1.1)
axes[1].legend();

Visualize the cavity final state

In [11]:
wigner_fock_distribution(ptrace(result.final_state, 0));

Threshold

In [12]:
d = (W12 - W21) / (W12 + W21)

d
Out[12]:
0.6666666666666666
In [13]:
gamma = 0.5 * (W12 + W21)

gamma
Out[13]:
0.030000000000000002
In [14]:
C = d * (g/(2*pi)) ** 2 / (gamma * kappa)

C
Out[14]:
44.44444444444445

Solve with mean-field contribution

$\displaystyle H_{\rm eff}(\rho, t) = \omega a^\dagger a + \frac{1}{2}\omega\sigma_z

  • ig\left({\rm Tr}[\sigma-\rho] a^\dagger - {\rm Tr}[\sigma+\rho] a\right)
  • ig\left({\rm Tr}[a^\dagger\rho] \sigma- - {\rm Tr}[a\rho] \sigma+\right) $
In [15]:
def Heff(t, rho, args):
    """ not an optimized implementation ... """
    q_rho = Qobj(vec2mat(rho))
    Heff = 1.0j * g * (expect(sm, q_rho) * ad - expect(sp, q_rho) * a) \
         + 1.0j * g * (expect(ad, q_rho) * sm - expect(a, q_rho) * sp)
    
    return args[0] + liouvillian_fast(Heff, []).data
In [16]:
sop_a_data = liouvillian_fast(a, []).data
sop_ad_data = liouvillian_fast(ad, []).data
sop_sm_data = liouvillian_fast(sm, []).data
sop_sp_data = liouvillian_fast(sp, []).data

sop_a_L_data = spre(a).data
sop_ad_L_data = spre(ad).data
sop_sm_L_data = spre(sm).data
sop_sp_L_data = spre(sp).data

def Heff_fast(t, rho, args):
    """ a somewhat more optimized version ... """
    
    phi_sm = expect_rho_vec1d(sop_sm_L_data, rho)
    phi_a = expect_rho_vec1d(sop_a_L_data, rho)
    
    Heff = 1.0j * g * (phi_sm * sop_ad_data  - conjugate(phi_sm) * sop_a_data) \
         + 1.0j * g * (conjugate(phi_a) * sop_sm_data - phi_a * sop_sp_data)
    
    return args[0] + Heff
In [17]:
result = mesolve(Heff_fast, psi0, tlist, c_ops, e_ops, args=[H],
                 options=Odeoptions(rhs_with_state=True, store_final_state=True),
                 progress_bar=TextProgressBar())
Completed:  0.0%. Elapsed time:   0.00s. Est. remaining time: 00:00:00:00.
Completed: 10.0%. Elapsed time:  12.88s. Est. remaining time: 00:00:01:55.
Completed: 20.0%. Elapsed time:  26.35s. Est. remaining time: 00:00:01:45.
Completed: 30.0%. Elapsed time:  39.72s. Est. remaining time: 00:00:01:32.
Completed: 40.0%. Elapsed time:  54.83s. Est. remaining time: 00:00:01:22.
Completed: 50.0%. Elapsed time:  72.08s. Est. remaining time: 00:00:01:12.
Completed: 60.0%. Elapsed time:  89.05s. Est. remaining time: 00:00:00:59.
Completed: 70.0%. Elapsed time: 106.85s. Est. remaining time: 00:00:00:45.
Completed: 80.0%. Elapsed time: 124.90s. Est. remaining time: 00:00:00:31.
Completed: 90.0%. Elapsed time: 143.07s. Est. remaining time: 00:00:00:15.
Elapsed time: 159.90s
In [18]:
fig, ax = subplots()

ax.plot(result.times, result.expect[0], label=r'$a^\dagger a$')
ax.plot(result.times, result.expect[1], label=r'$\sigma_+\sigma_-$')
ax.set_ylim(-0.1, 3)
ax.legend();

Visualize the cavity final state

In [19]:
wigner_fock_distribution(ptrace(result.final_state, 0));

Solve the same problem using the list-function format

In [20]:
def a_coeff(t, rho, args):
    return - 1.0j * g * expect_rho_vec1d(sop_sp_L_data, rho)

def ad_coeff(t, rho, args):
    return + 1.0j * g * expect_rho_vec1d(sop_sm_L_data, rho)

def sm_coeff(t, rho, args):
    return + 1.0j * g * expect_rho_vec1d(sop_ad_L_data, rho)

def sp_coeff(t, rho, args):
    return - 1.0j * g * expect_rho_vec1d(sop_a_L_data, rho)
In [21]:
result = mesolve([H, [a, a_coeff], [ad, ad_coeff], [sm, sm_coeff], [sp, sp_coeff]],
                 psi0, tlist, c_ops, e_ops, args={}, 
                 options=Odeoptions(rhs_with_state=True, store_final_state=True),
                 progress_bar=TextProgressBar())
Completed:  0.0%. Elapsed time:   0.00s. Est. remaining time: 00:00:00:00.
Completed: 10.0%. Elapsed time:  16.56s. Est. remaining time: 00:00:02:29.
Completed: 20.0%. Elapsed time:  33.92s. Est. remaining time: 00:00:02:15.
Completed: 30.0%. Elapsed time:  50.92s. Est. remaining time: 00:00:01:58.
Completed: 40.0%. Elapsed time:  70.16s. Est. remaining time: 00:00:01:45.
Completed: 50.0%. Elapsed time:  92.30s. Est. remaining time: 00:00:01:32.
Completed: 60.0%. Elapsed time: 114.58s. Est. remaining time: 00:00:01:16.
Completed: 70.0%. Elapsed time: 136.75s. Est. remaining time: 00:00:00:58.
Completed: 80.0%. Elapsed time: 158.77s. Est. remaining time: 00:00:00:39.
Completed: 90.0%. Elapsed time: 180.88s. Est. remaining time: 00:00:00:20.
Elapsed time: 202.96s
In [22]:
fig, ax = subplots()

ax.plot(result.times, result.expect[0], label=r'$a^\dagger a$')
ax.plot(result.times, result.expect[1], label=r'$\sigma_+\sigma_-$')
ax.set_ylim(-0.1, 3)
ax.legend();

Solve with mean-field contribution, but without dissipation (using sesolve)

In [23]:
tlist = linspace(0, 50, 150)
In [24]:
H_data = H.data
a_data = a.data
ad_data = ad.data
sm_data = sm.data
sp_data = sp.data

def Heff_func(t, psi, args):
    
    phi_sm = expect_psi(sm_data, psi[:, newaxis])
    phi_a  = expect_psi(a_data, psi[:, newaxis])
    phi_sp = expect_psi(sp_data, psi[:, newaxis])
    phi_ad = expect_psi(ad_data, psi[:, newaxis])
    
    Heff = 1.0j * g * (phi_sm * ad_data - phi_sp * a_data) \
         + 1.0j * g * (phi_ad * sm_data - phi_a * sp_data)
    
    return H_data + Heff
In [25]:
result = sesolve(Heff_func, psi0, tlist, e_ops,
                 options=Odeoptions(rhs_with_state=True, store_final_state=True),
                 progress_bar=TextProgressBar())
Completed:  0.0%. Elapsed time:   0.00s. Est. remaining time: 00:00:00:00.
Completed: 10.0%. Elapsed time:   1.82s. Est. remaining time: 00:00:00:16.
Completed: 20.0%. Elapsed time:   3.62s. Est. remaining time: 00:00:00:14.
Completed: 30.0%. Elapsed time:   5.42s. Est. remaining time: 00:00:00:12.
Completed: 40.0%. Elapsed time:   7.22s. Est. remaining time: 00:00:00:10.
Completed: 50.0%. Elapsed time:   9.03s. Est. remaining time: 00:00:00:09.
Completed: 60.0%. Elapsed time:  10.84s. Est. remaining time: 00:00:00:07.
Completed: 70.0%. Elapsed time:  12.67s. Est. remaining time: 00:00:00:05.
Completed: 80.0%. Elapsed time:  14.48s. Est. remaining time: 00:00:00:03.
Completed: 90.0%. Elapsed time:  16.32s. Est. remaining time: 00:00:00:01.
Elapsed time:  18.13s
In [26]:
fig, ax = subplots()

ax.plot(result.times, real(result.expect[0]), label=r'$a^\dagger a$')
ax.plot(result.times, real(result.expect[1]), label=r'$\sigma_+\sigma_-$')
ax.set_ylim(-0.1, 1)
ax.legend();

Solve the same problem with list function format

In [27]:
def a_coeff(t, psi, args):
    return - 1.0j * g * expect_psi(sp_data, psi[:, newaxis])

def ad_coeff(t, psi, args):
    return + 1.0j * g * expect_psi(sm_data, psi[:, newaxis])

def sm_coeff(t, psi, args):
    return + 1.0j * g * expect_psi(ad_data, psi[:, newaxis])

def sp_coeff(t, psi, args):
    return - 1.0j * g * expect_psi(a_data, psi[:, newaxis])
In [28]:
result = sesolve([H, [a, a_coeff], [a.dag(), ad_coeff], [sm, sm_coeff], [sp, sp_coeff]],
                 psi0, tlist, e_ops, 
                 options=Odeoptions(rhs_with_state=True, store_final_state=True),
                 progress_bar=TextProgressBar())
Completed:  0.0%. Elapsed time:   0.00s. Est. remaining time: 00:00:00:00.
Completed: 10.0%. Elapsed time:   1.89s. Est. remaining time: 00:00:00:16.
Completed: 20.0%. Elapsed time:   3.76s. Est. remaining time: 00:00:00:15.
Completed: 30.0%. Elapsed time:   5.63s. Est. remaining time: 00:00:00:13.
Completed: 40.0%. Elapsed time:   7.52s. Est. remaining time: 00:00:00:11.
Completed: 50.0%. Elapsed time:   9.39s. Est. remaining time: 00:00:00:09.
Completed: 60.0%. Elapsed time:  11.25s. Est. remaining time: 00:00:00:07.
Completed: 70.0%. Elapsed time:  13.10s. Est. remaining time: 00:00:00:05.
Completed: 80.0%. Elapsed time:  14.96s. Est. remaining time: 00:00:00:03.
Completed: 90.0%. Elapsed time:  16.85s. Est. remaining time: 00:00:00:01.
Elapsed time:  18.73s
In [29]:
fig, ax = subplots()

ax.plot(result.times, real(result.expect[0]), label=r'$a^\dagger a$')
ax.plot(result.times, real(result.expect[1]), label=r'$\sigma_+\sigma_-$')
ax.set_ylim(-0.1, 1)
ax.legend();

Solve with mean-field contribution and dissipation using mcsolve

In [30]:
tlist = linspace(0, 250, 250)
In [31]:
result = mcsolve(Heff_func, psi0, tlist, c_ops, e_ops, ntraj=64,
                 options=Odeoptions(rhs_with_state=True, store_final_state=True, gui=False, num_cpus=4),
                 progress_bar=TextProgressBar())
Completed:  1.6%. Elapsed time: 141.29s. Est. remaining time: 00:02:28:21.
Completed: 10.9%. Elapsed time: 505.69s. Est. remaining time: 00:01:08:37.
Completed: 20.3%. Elapsed time: 884.28s. Est. remaining time: 00:00:57:49.
Completed: 31.2%. Elapsed time: 1122.07s. Est. remaining time: 00:00:41:08.
Completed: 40.6%. Elapsed time: 1404.55s. Est. remaining time: 00:00:34:12.
Completed: 50.0%. Elapsed time: 1624.73s. Est. remaining time: 00:00:27:04.
Completed: 60.9%. Elapsed time: 1885.82s. Est. remaining time: 00:00:20:08.
Completed: 70.3%. Elapsed time: 2173.70s. Est. remaining time: 00:00:15:17.
Completed: 81.2%. Elapsed time: 2467.79s. Est. remaining time: 00:00:09:29.
Completed: 90.6%. Elapsed time: 2683.16s. Est. remaining time: 00:00:04:37.
Completed: 100.0%. Elapsed time: 2915.00s. Est. remaining time: 00:00:00:00.
Elapsed time: 2915.10s
In [32]:
fig, ax = subplots()

ax.plot(result.times, real(result.expect[0]), label=r'$a^\dagger a$')
ax.plot(result.times, real(result.expect[1]), label=r'$\sigma_+\sigma_-$')
ax.set_ylim(-0.1, 3)
ax.legend();

Using the list function format and mcsolve

In [33]:
result = mcsolve([H, [a, a_coeff], [ad, ad_coeff], [sm, sm_coeff], [sp, sp_coeff]],
                 psi0, tlist, c_ops, e_ops, args={}, ntraj=512,
                 options=Odeoptions(rhs_with_state=True, store_final_state=True, gui=False, num_cpus=4),
                 progress_bar=TextProgressBar())
Completed:  0.2%. Elapsed time:  37.77s. Est. remaining time: 00:05:21:39.
Completed: 10.2%. Elapsed time: 518.60s. Est. remaining time: 00:01:16:27.
Completed: 20.1%. Elapsed time: 1019.87s. Est. remaining time: 00:01:07:29.
Completed: 30.1%. Elapsed time: 1520.25s. Est. remaining time: 00:00:58:54.
Completed: 40.0%. Elapsed time: 2034.98s. Est. remaining time: 00:00:50:47.
Completed: 50.0%. Elapsed time: 2513.66s. Est. remaining time: 00:00:41:53.
Completed: 60.2%. Elapsed time: 3036.83s. Est. remaining time: 00:00:33:31.
Completed: 70.1%. Elapsed time: 3589.60s. Est. remaining time: 00:00:25:29.
Completed: 80.1%. Elapsed time: 4684.49s. Est. remaining time: 00:00:19:25.
Completed: 90.0%. Elapsed time: 5647.46s. Est. remaining time: 00:00:10:24.
Completed: 100.0%. Elapsed time: 6778.12s. Est. remaining time: 00:00:00:00.
Elapsed time: 6778.18s
In [34]:
fig, ax = subplots()

ax.plot(result.times, real(result.expect[0]), label=r'$a^\dagger a$')
ax.plot(result.times, real(result.expect[1]), label=r'$\sigma_+\sigma_-$')
ax.set_ylim(-0.1, 3)
ax.legend();

Software versions

In [35]:
from qutip.ipynbtools import version_table

version_table()
Out[35]:
SoftwareVersion
matplotlib1.4.x
Python3.3.1 (default, Apr 17 2013, 22:30:32) [GCC 4.7.3]
SciPy0.13.0.dev-d74fd00
Cython0.19.1
QuTiP2.3.0.dev-b354119
Numpy1.8.0.dev-75cdf3d
IPython1.0.dev
OSposix [linux]
Thu Jul 04 15:30:06 2013 JST