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