Author: Neill Lambert (nwlambert@gmail.com) Anubhav Vardhan (anubhavvardhan@gmail.com)
For more information about QuTiP see http://qutip.org
Reference: link.aps.org/doi/10.1103/PhysRevA.90.032114
In this model, we take a quantum system coupled to a bosonic environment and map to a model in which a collective mode of the environment, known as the reaction coordinate (RC), is incorporated within an effective system Hamiltonian. We then treat the residual environment within a full second-order Born-Markov master equation formalism. Thus, all important system-bath, and indeed intrabath, correlations are incorporated into the system-RC Hamiltonian in the regimes we study.
In [1]:
%pylab inline
from qutip import *
In [2]:
from qutip import rcsolve
In [3]:
Del = 1.0 # The number of qubits in the system.
wq = 0.5 # Energy of the 2-level system.
Hsys = 0.5 * wq * sigmaz() + 0.5 * Del * sigmax()
In [4]:
Q = sigmaz()
In [8]:
wc = 0.05 # Cutoff frequency.
alpha = 2.5/pi # Coupling strength.
N = 10 # Number of cavity fock states.
Temperature = 1/0.95 # Tempertaure.
tlist = np.linspace(0, 40, 600)
initial_state = basis(2,1) * basis(2,1).dag() # Initial state of the system.
return_vals = [Q] # List for which to calculate expectation value
options = Options(nsteps=15000, store_states=True) # Options for the solver.
output = rcsolve(Hsys, initial_state, tlist, return_vals, Q, wc, alpha, N,
Temperature, options=options)
In [9]:
fig, axes = subplots(1, 1, sharex=True, figsize=(8,4))
axes.plot(tlist, real(output.expect[0]), 'b', linewidth=2, label="P12")
axes.legend(loc=0)
Out[9]:
In [10]:
output.states[0]
Out[10]:
In [11]:
t_idx_vec = range(0,len(tlist),200)
fig, axes = subplots(len(t_idx_vec), 1, sharey=True, figsize=(8,2*len(t_idx_vec)))
for idx, t_idx in enumerate(t_idx_vec):
psi_a = ptrace(output.states[t_idx], 0)
cont1 = axes[idx].bar(range(0, N), real(psi_a.diag()))
fig.tight_layout()
In [12]:
xvec = linspace(-5,5,200)
t_idx_vec = range(0,len(tlist),200)
fig, axes = subplots(len(t_idx_vec), 1, sharex=True, sharey=True, figsize=(8,4*len(t_idx_vec)))
for idx, t_idx in enumerate(t_idx_vec):
psi_a = ptrace(output.states[t_idx], 0)
W_a = wigner(psi_a, xvec, xvec)
cont1 = axes[idx].contourf(xvec, xvec, W_a, 100)
In [13]:
from qutip.ipynbtools import version_table
version_table()
Out[13]:
In [ ]: