Optical Bloch equations: Lindblad dynamics using the matvec function

This example uses the tools.misc.matvec() function to define a Lindblad equation and solve the ODE using the tools.evolution.evolve() function.

Consider the the two-level system:

\[\begin{split}H &= \delta\sigma^z + \Omega_0\sigma^x, \\ L &= \sigma^+.\end{split}\]

where \(H\) is the Hamiltonian of the two-level system, and \(L\) is the Lindblad (or jump) operator. The Lindblad equation is a non-unitary extension of the Liouville-von Neumann equation for the density matrix \(\rho(t)\):

\[\partial_t\rho(t) = -i[H,\rho(t)] + 2\gamma\left(L\rho(t)L^\dagger - \frac{1}{2}\{L^\dagger L,\rho(t) \} \right),\]

where \([\cdot,\cdot]\) is the commutator, and \(\{\cdot,\cdot\}\) is the anti-commutator. The system of equations for this specific problem is also known as the optical Bloch equations.

Below, we provide three ways to define the function for the Lindblad ODE:

  • the first version is very intuitive as it follows the definition of the Lindblad equation above; but it is rather slow;

  • the second version uses the hamiltonian.dot() and hamiltonian.rdot() functions and is a bit more sophisticated and a bit faster;

  • the third version uses the matvec() function and is faster than the previous two (but may not be memory efficient for large systems). Notice that there are different kinds of matvec functions for sparse matrices, depending on the format (csc, csr, etc.). In using them, one should keep in mind the format of the matrix, which can be determined using the get_matvec_function.

Note that this way of simulating the Lindblad equation has severe limitations for many-body systems. An alternative, parallelizable way to effectively simulate a subset of Lindblad dynamics using unitary evolution is described in arXiv:1608.01317.

Script

download script

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
from __future__ import print_function, division
#
import sys,os
os.environ['KMP_DUPLICATE_LIB_OK']='True' # uncomment this line if omp error occurs on OSX for python 3
os.environ['OMP_NUM_THREADS']='1' # set number of OpenMP threads to run in parallel
os.environ['MKL_NUM_THREADS']='1' # set number of MKL threads to run in parallel
#
quspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,quspin_path)
###########################################################################
#                            example 17                                   #
#  In this script we demonstrate how to apply the matvec function         #
#  to define the Lundblad equation for a two-leel system, and solve it    #
#  using the evolve funcion.                                              #
###########################################################################
from quspin.operators import hamiltonian, commutator, anti_commutator
from quspin.basis import spin_basis_1d # Hilbert space spin basis_1d
from quspin.tools.evolution import evolve
from quspin.tools.misc import get_matvec_function
import numpy as np
from six import iteritems # loop over elements of dictionary
import matplotlib.pyplot as plt # plotting library
#
###### model parameters
#
L=1 # one qubit
delta=1.0 # detuning
Omega_0=np.sqrt(2) # bare Rabi frequency
Omega_Rabi=np.sqrt(Omega_0**2+delta**2) # Rabi frequency
gamma = 0.5*np.sqrt(3) # decay rate
#
##### create Hamiltonian to evolve unitarily
# basis
basis=spin_basis_1d(L,pauli=-1) # uses convention "+" = [[0,1],[0,0]]
# site-coupling lists
hx_list=[[Omega_0,i] for i in range(L)]
hz_list=[[delta,i] for i in range(L)]
# static opstr list 
static_H=[['x',hx_list],['z',hz_list]]
# hamiltonian
H=hamiltonian(static_H,[],basis=basis,dtype=np.float64)
print('Hamiltonian:\n',H.toarray())
#
##### create Lindbladian
# site-coupling lists
L_list=[[1.0j,i] for i in range(L)]
# static opstr list 
static_L=[['+',L_list]]
# Lindblad operator
L=hamiltonian(static_L,[],basis=basis,dtype=np.complex128,check_herm=False)
print('Lindblad operator:\n',L.toarray())
# pre-compute operators for efficiency
L_dagger=L.getH()
L_daggerL=L_dagger*L
#
#### determine the corresponding matvec routines ####
#
# different matvec functions are required since we use both matrices and their transposed, cf. Lindblad_EOM_v3
matvec_csr=get_matvec_function(H.static) # csr matvec function
matvec_csc=get_matvec_function(H.static.T) # csc matvec function
#
##### define Lindblad equation in diagonal form
#
# slow, straightforward function
#
def Lindblad_EOM_v1(time,rho):
	"""
	This function solves the complex-valued time-dependent GPE:
	$$ \dot\rho(t) = -i[H,\rho(t)] + 2\gamma\left( L\rho L^\dagger - \frac{1}{2}\{L^\dagger L, \rho \} \right) $$
	"""
	# solve static part of Lindblad equation
	rho = rho.reshape((H.Ns,H.Ns))
	rho_dot = -1j*commutator(H,rho).static + 2.0*gamma*(L*rho*L_dagger).static - gamma*anti_commutator(L_daggerL,rho).static 
	# solve dynamic part of Lindblad equation (no time-dependence in Lindbladian for this example)
	for f,Hd in iteritems(H.dynamic):
		rho_dot += -1j*f(time)*commutator(Hd,rho)
	return rho_dot.ravel()
#
# intermediate, straightforward function using dot and rdot
#
def Lindblad_EOM_v2(time,rho,rho_out,rho_aux):
	"""
	This function solves the complex-valued time-dependent GPE:
	$$ \dot\rho(t) = -i[H,\rho(t)] + 2\gamma\left( L\rho L^\dagger - \frac{1}{2}\{L^\dagger L, \rho \} \right) $$
	"""
	rho = rho.reshape((H.Ns,H.Ns)) # reshape vector from ODE solver input
	### Hamiltonian part
	# commutator term (unitary)
	# rho_out = H.static.dot(rho))
	H.dot( rho,out=rho_out  ,a=+1.0,overwrite_out=True)
	# rho_out -= (H.static.T.dot(rho.T)).T // RHS~rho.dot(H) 
	H.rdot(rho,out=rho_out,a=-1.0,overwrite_out=False)
	# multiply by -i
	rho_out *= -1.0j
	#
	### Lindbladian part (static only)
	## 1st Lindblad term (nonunitary)
	# rho_aux = 2\gamma*L.dot(rho)
	L.dot(rho       ,out=rho_aux, a=+2.0*gamma,overwrite_out=True)
	# rho_out += (L.static.conj().dot(rho_aux.T)).T // RHS~rho_aux.dot(L_dagger) 
	L.H.rdot(rho_aux,out=rho_out, a=+1.0,overwrite_out=False) # L.H = L^\dagger
	#
	## anticommutator (2nd Lindblad) term (nonunitary)
	# rho_out += gamma*L_daggerL._static.dot(rho)
	L_daggerL.dot(rho ,out=rho_out, a=-gamma,overwrite_out=False)
	# # rho_out += gamma*(L_daggerL._static.T.dot(rho.T)).T // RHS~rho.dot(L_daggerL) 
	L_daggerL.rdot(rho,out=rho_out, a=-gamma,overwrite_out=False) 
	
	return rho_out.ravel() # ODE solver accepts vectors only
#
# fast function using matvec (not as memory efficient)
#
def Lindblad_EOM_v3(time,rho,rho_out,rho_aux):
	"""
	This function solves the complex-valued time-dependent GPE:
	$$ \dot\rho(t) = -i[H,\rho(t)] + 2\gamma\left( L\rho L^\dagger - \frac{1}{2}\{L^\dagger L, \rho \} \right) $$
	"""
	rho = rho.reshape((H.Ns,H.Ns)) # reshape vector from ODE solver input
	### Hamiltonian part
	# commutator term (unitary
	# rho_out = H.static.dot(rho))
	matvec_csr(H.static  ,rho  ,out=rho_out  ,a=+1.0,overwrite_out=True)
	# rho_out -= (H.static.T.dot(rho.T)).T // RHS~rho.dot(H) 
	matvec_csc(H.static.T,rho.T,out=rho_out.T,a=-1.0,overwrite_out=False)
	# 
	for func,Hd in iteritems(H._dynamic):
		ft = func(time)
		# rho_out += ft*Hd.dot(rho)
		matvec_csr(Hd  ,rho  ,out=rho_out  ,a=+ft,overwrite_out=False)
		# rho_out -= ft*(Hd.T.dot(rho.T)).T 
		matvec_csc(Hd.T,rho.T,out=rho_out.T,a=-ft,overwrite_out=False)
	# multiply by -i
	rho_out *= -1.0j
	#
	### Lindbladian part (static only)
	## 1st Lindblad term (nonunitary): 2\gamma * L*rho*L^\dagger
	# rho_aux = 2\gamma*L.dot(rho)
	matvec_csr(L.static  ,rho             ,out=rho_aux  ,a=+2.0*gamma,overwrite_out=True)
	# rho_out += (L.static.conj().dot(rho_aux.T)).T // RHS ~ rho_aux.dot(L_dagger)
	matvec_csr(L.static.conj(),rho_aux.T,out=rho_out.T,a=+1.0,overwrite_out=False)
	# 
	## anticommutator (2nd Lindblad) term (nonunitary): -\gamma \{L^\dagger * L, \eho}
	# rho_out += gamma*L_daggerL.static.dot(rho)
	matvec_csr(L_daggerL.static,rho  ,out=rho_out  ,a=-gamma,overwrite_out=False)
	# rho_out += gamma*(L_daggerL.static.T.dot(rho.T)).T // RHS~rho.dot(L_daggerL) 
	matvec_csc(L_daggerL.static.T,rho.T,out=rho_out.T,a=-gamma,overwrite_out=False) 
	#
	return rho_out.ravel() # ODE solver accepts vectors only
#
# define auxiliary arguments
EOM_args=(np.zeros((H.Ns,H.Ns),dtype=np.complex128,order="C"),    # auxiliary variable rho_out
		  np.zeros((H.Ns,H.Ns),dtype=np.complex128,order="C"),  ) # auxiliary variable rho_aux
#
##### time-evolve state according to Lindlad equation
# define real time vector
t_max=6.0
time=np.linspace(0.0,t_max,101)
# define initial state
rho0=np.array([[0.5,0.5j],[-0.5j,0.5]],dtype=np.complex128)
# slow solution, uses Lindblad_EOM_v1
#rho_t = evolve(rho0,time[0],time,Lindblad_EOM_v1,iterate=True,atol=1E-12,rtol=1E-12)
# intermediate function, uses Lindblad_EOM_v2
rho_t = evolve(rho0,time[0],time,Lindblad_EOM_v2,f_params=EOM_args,iterate=True,atol=1E-12,rtol=1E-12) 
# fast solution (but 3 times as memory intensive), uses Lindblad_EOM_v3
#rho_t = evolve(rho0,time[0],time,Lindblad_EOM_v3,f_params=EOM_args,iterate=True,atol=1E-12,rtol=1E-12) 
#
# compute state evolution
population_down=np.zeros(time.shape,dtype=np.float64)
for i,rho_flattened in enumerate(rho_t):
	rho=rho_flattened.reshape(H.Ns,H.Ns)
	population_down[i]=rho_flattened[1,1].real
	print("time={0:.2f}, population of down state = {1:0.8f}".format(time[i],population_down[i]) )
#
##### plot population dynamics of down state
#
plt.plot(Omega_Rabi*time, population_down)
plt.xlabel('$\\Omega_R t$')
plt.ylabel('$\\rho_{\\downarrow\\downarrow}$')
plt.xlim(time[0],time[-1])
plt.ylim(0.0,1.0)
plt.grid()
plt.tight_layout()
plt.savefig('example17.pdf', bbox_inches='tight')
plt.close()