stochastic quantum dynamics with QuTiP's mcsolve

stochastic

#!/usr/bin/env python
# coding: utf-8

#  # Stochastic Quantum Dynamics with QuTiP
#  
#  
#  We use QuTiP's solvers to study the open dynamics of a quantum system evolving in time. 

# In[12]:


import numpy as np
import matplotlib.pyplot as plt
from qutip import *


# ### Define the operators and Hamiltonian

# In[13]:


# spins
sx_reduced = sigmax()
sy_reduced = sigmay()
sz_reduced = sigmaz()
sp_reduced = sigmap()
sm_reduced = sigmam()

# photons
nph = 4
a_reduced = destroy(nph)

# tensor space
sz = tensor(sz_reduced,qeye(nph))
sx = tensor(sx_reduced,qeye(nph))
sm = tensor(sm_reduced,qeye(nph))
sp = sm.dag()
a = tensor(qeye(2), a_reduced)

# hamiltonians
g = 0.1
Hcav = a.dag()*a
wx = 0.5
Hspin = sz + wx*sx
Hint =  g*sx*(a+a.dag())
HintCheck =  g*tensor(sx_reduced,a_reduced+a_reduced.dag())
H = Hcav + Hspin + Hint
np.testing.assert_(Hint == HintCheck)


# ### Define the initial state

# In[14]:


# initial state
psi0_spin = basis(2,0)
psi0_phot = basis(nph,nph-int(nph/2))
psi0 = tensor(psi0_spin,psi0_phot)
rho0 = ket2dm(psi0)

tlist = np.linspace(0,10,1000)


# In[19]:


kappa = 0.1
gamma = 0.1
# set dynamics options
my_options = Options(average_states = False, store_states = True)
# solve dynamics
results_mc = mcsolve(H, psi0, tlist, 
                  c_ops=[np.sqrt(kappa)*a,np.sqrt(gamma)*sz], 
                  e_ops=[a.dag()*a,sz], 
                  options=my_options, progress_bar=True)

# store time evoluted variables
nph_mc_t = results_mc.expect[0]
sz_mc_t = results_mc.expect[1]


# In[25]:


rho_mc_t = results_mc.states
len(rho_mc_t)
help(expect)


# In[26]:


sz_stoch_t = []
a_stoch_t = []
for i in range(len(rho_mc_t)):
    sz_stoch_t.append(expect(sz,rho_mc_t[i]))
    nph_stoch_t.append(expect(a.dag()*a,rho_mc_t[i]))


# In[16]:


my_options = Options(average_states = True, store_states = True)
results = mesolve(H, psi0, tlist, 
                  c_ops=[np.sqrt(kappa)*a,np.sqrt(gamma)*sz], 
                  e_ops=[a.dag()*a,sz], 
                  options=my_options, progress_bar=True)

# store time evoluted variables
nph_t = results.expect[0]
sz_t = results.expect[1]
rho_t = results.states


# In[6]:


plt.figure()
plt.plot(tlist, sz_t)
plt.plot(tlist, sz_mc_t)
plt.plot(tlist, nph_t/nph/0.5)
plt.plot(tlist, nph_mc_t/nph/0.5)
plt.xlabel(r"$t$")
plt.ylabel(r"$\langle X \rangle$")
plt.show()
plt.close()


# In[11]:


plt.figure()
plt.plot(tlist, sz_t)
plt.plot(tlist, sz_mc_t)
plt.plot(tlist, nph_t/nph/0.5)
plt.plot(tlist, nph_mc_t/nph/0.5)
for i in range(len(sz_stoch_t)):
    plt.plot(tlist, sz_t)
    plt.plot(tlist, sz_mc_t)
    plt.plot(tlist, nph_t/nph/0.5)
    plt.plot(tlist, nph_mc_t/nph/0.5)
plt.xlabel(r"$t$")
plt.ylabel(r"$\langle X \rangle$")
plt.show()
plt.close()


# In[ ]:


results = mesolve(H, rho0, tlist, 
                  c_ops=[np.sqrt(kappa)*a,np.sqrt(gamma)*sz], 
                  e_ops=[a.dag()*a,sz], 
                  options=my_options, progress_bar=None)

# store time evoluted variables
nph_t = results.expect[0]
sz_t = results.expect[1]
rho_t = results.states


# In[ ]:





# In[9]:


qutip.about()