In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
import numpy as np
In [3]:
from qutip import *
In [4]:
delta = 0.0 * 2 * np.pi
epsilon = 0.5 * 2 * np.pi
gamma = 0.25
times = np.linspace(0, 10, 100)
In [5]:
H = delta/2 * sigmax() + epsilon/2 * sigmaz()
H
Out[5]:
In [6]:
psi0 = (2 * basis(2, 0) + basis(2, 1)).unit()
In [7]:
c_ops = [np.sqrt(gamma) * sigmam()]
a_ops = [sigmax()]
In [8]:
e_ops = [sigmax(), sigmay(), sigmaz()]
In [9]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [10]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, spectra_cb=[lambda w : gamma * (w > 0)])
In [11]:
plot_expectation_values([result_me, result_brme]);
In [12]:
b = Bloch()
b.add_points(result_me.expect, meth='l')
b.add_points(result_brme.expect, meth='l')
b.make_sphere()
In [13]:
N = 10
w0 = 1.0 * 2 * np.pi
g = 0.05 * w0
kappa = 0.15
times = np.linspace(0, 25, 1000)
In [14]:
a = destroy(N)
In [15]:
H = w0 * a.dag() * a + g * (a + a.dag())
In [16]:
# start in a superposition state
psi0 = ket2dm((basis(N, 4) + basis(N, 2) + basis(N,0)).unit())
In [17]:
c_ops = [np.sqrt(kappa) * a]
a_ops = [a + a.dag()]
In [18]:
e_ops = [a.dag() * a, a + a.dag()]
In [19]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [20]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, spectra_cb=[lambda w : kappa * (w > 0)])
In [21]:
plot_expectation_values([result_me, result_brme]);
In [22]:
times = np.linspace(0, 25, 250)
In [23]:
n_th = 1.5
c_ops = [np.sqrt(kappa * (n_th + 1)) * a, np.sqrt(kappa * n_th) * a.dag()]
In [24]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [25]:
w_th = w0/np.log(1 + 1/n_th)
def S_w(w):
if w >= 0:
return (n_th + 1) * kappa
else:
return (n_th + 1) * kappa * np.exp(w / w_th)
In [26]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, [S_w])
In [27]:
plot_expectation_values([result_me, result_brme]);
In [28]:
result_me = mesolve(H, psi0, times, c_ops, [])
In [29]:
result_brme = brmesolve(H, psi0, times, a_ops, [], [S_w])
In [30]:
n_me = expect(a.dag() * a, result_me.states)
In [31]:
n_brme = expect(a.dag() * a, result_brme.states)
In [32]:
fig, ax = plt.subplots()
ax.plot(times, n_me, label='me')
ax.plot(times, n_brme, label='brme')
ax.legend()
ax.set_xlabel("t");
In [33]:
N = 10
a = tensor(destroy(N), identity(2))
sm = tensor(identity(N), destroy(2))
psi0 = ket2dm(tensor(basis(N, 1), basis(2, 0)))
a_ops = [(a + a.dag())]
e_ops = [a.dag() * a, sm.dag() * sm]
In [34]:
w0 = 1.0 * 2 * np.pi
g = 0.05 * 2 * np.pi
kappa = 0.05
times = np.linspace(0, 5 * 2 * np.pi / g, 1000)
In [35]:
c_ops = [np.sqrt(kappa) * a]
H = w0 * a.dag() * a + w0 * sm.dag() * sm + g * (a + a.dag()) * (sm + sm.dag())
In [36]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [37]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, spectra_cb=[lambda w : kappa*(w > 0)])
In [38]:
plot_expectation_values([result_me, result_brme]);
In the weak coupling regime there is no significant difference between the Lindblad master equation and the Bloch-Redfield master equation.
In [39]:
w0 = 1.0 * 2 * np.pi
g = 0.75 * 2 * np.pi
kappa = 0.05
times = np.linspace(0, 5 * 2 * np.pi / g, 1000)
In [40]:
c_ops = [np.sqrt(kappa) * a]
H = w0 * a.dag() * a + w0 * sm.dag() * sm + g * (a + a.dag()) * (sm + sm.dag())
In [41]:
result_me = mesolve(H, psi0, times, c_ops, e_ops)
In [42]:
result_brme = brmesolve(H, psi0, times, a_ops, e_ops, spectra_cb=[lambda w : kappa*(w > 0)])
In [43]:
plot_expectation_values([result_me, result_brme]);
In the strong coupling regime there are some corrections to the Lindblad master equation that is due to the fact system eigenstates are hybridized states with both atomic and cavity contributions.
In [44]:
from qutip.ipynbtools import version_table
version_table()
Out[44]: