Reaction coordinate solver for 2-level bosonic system

Author: Neill Lambert (nwlambert@gmail.com) Anubhav Vardhan (anubhavvardhan@gmail.com)

For more information about QuTiP see http://qutip.org

Introduction

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 *









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
from qutip import rcsolve

Defining the 2-level system



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()

Defining the coupling Q such that:

Decoupled Hamiltonian-

$H_0 = \omega a^\dagger a + \frac{1}{2}\omega_q\sigma_z + \frac{1}{2}\delta\sigma_x$

Interaction-

$H_1 = g(a^\dagger + a)\sigma_z$

Total Hamiltonian-

$H = H_0 + H_1$

Solving the RC master equation



In [5]:

    
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 = [tensor(qeye(N), kk) for kk in [Q]]            # List for which to calculate expectation value
eigen_sparse = False
calc_time = True                                             
options = Options(nsteps=15000, store_states=True)        # Options for the solver.

output = rcsolve(Hsys, Q, wc, alpha, N, Temperature, tlist, initial_state,
                 return_vals, eigen_sparse, calc_time, options)









    



Integration required 23.2765 seconds

Plotting the TLS state occupation



In [6]:

    
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[6]:





<matplotlib.legend.Legend at 0x7f8fd23ce3d0>

Plotting the photon distributions at arbitrary times



In [7]:

    
output.states[0]









    Out[7]:





Quantum object: dims = [[10, 2], [10, 2]], shape = [20, 20], type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.654 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.226 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \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 & 3.862\times10^{-04} & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 1.335\times10^{-04} & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 4.616\times10^{-05}\\\end{array}\right)\end{equation*}



In [8]:

    
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()

Plotting the wigner function at arbitrary times



In [9]:

    
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)

Software versions



In [10]:

    
from qutip.ipynbtools import version_table

version_table()









    Out[10]:




Software Version
QuTiP 3.2.0.dev-8f37ac6
Numpy 1.8.0
SciPy 0.14.1
matplotlib 1.3.1
Cython 0.22
IPython 3.0.0
Python 2.7.9 (default, Apr  2 2015, 15:33:21) 
[GCC 4.9.2]
OS posix [linux2]
Fri Jun 05 15:19:05 2015 JST

Software	Version
QuTiP	3.2.0.dev-8f37ac6
Numpy	1.8.0
SciPy	0.14.1
matplotlib	1.3.1
Cython	0.22
IPython	3.0.0
Python	2.7.9 (default, Apr 2 2015, 15:33:21) [GCC 4.9.2]
OS	posix [linux2]
Fri Jun 05 15:19:05 2015 JST