Copyright (C) 2011 and later, Paul D. Nation & Robert J. Johansson
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
from qutip import *
from qutip.ipynbtools import HTMLProgressBar
from qutip.gui.progressbar import TextProgressBar
Consider the mean-field laser master equation on the form (Breuer and Petruccione)
$\displaystyle \dot\rho = -i [H, \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]
where
$\displaystyle H_{\rm eff}(\rho, t) = \omega a^\dagger a + \frac{1}{2}\omega\sigma_z
$
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.
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))
#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
theta = 0.0
H = w * ad * a - 0.5 * w * (cos(2*theta) * sz + sin(2*theta) * sx)
H
e_ops = [ad * a, sm.dag() * sm, sm, sp]
c_ops = [sqrt(2*kappa) * a, sqrt(W21) * sm, sqrt(W12) * sp]
result = mesolve(H, psi0, tlist, c_ops, e_ops, options=Odeoptions(store_final_state=True))
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();
wigner_fock_distribution(ptrace(result.final_state, 0));
d = (W12 - W21) / (W12 + W21)
d
0.6666666666666666
gamma = 0.5 * (W12 + W21)
gamma
0.030000000000000002
C = d * (g/(2*pi)) ** 2 / (gamma * kappa)
C
44.44444444444445
$\displaystyle H_{\rm eff}(\rho, t) = \omega a^\dagger a + \frac{1}{2}\omega\sigma_z
$
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
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
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
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();
wigner_fock_distribution(ptrace(result.final_state, 0));
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)
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
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();
tlist = linspace(0, 50, 150)
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
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
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();
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])
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
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();
tlist = linspace(0, 250, 250)
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
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();
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
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();
from qutip.ipynbtools import version_table
version_table()
Software | Version |
---|---|
matplotlib | 1.4.x |
Python | 3.3.1 (default, Apr 17 2013, 22:30:32) [GCC 4.7.3] |
SciPy | 0.13.0.dev-d74fd00 |
Cython | 0.19.1 |
QuTiP | 2.3.0.dev-b354119 |
Numpy | 1.8.0.dev-75cdf3d |
IPython | 1.0.dev |
OS | posix [linux] |
Thu Jul 04 15:30:06 2013 JST |