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
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]:
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]:
In [13]:
gamma = 0.5 * (W12 + W21)
gamma
Out[13]:
In [14]:
C = d * (g/(2*pi)) ** 2 / (gamma * kappa)
C
Out[14]:
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())
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())
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())
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())
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())
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())
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]: