Bloch-Redfield Solver: Two Level System¶
Author: C.Staufenbiel, 2022
with inspirations from the brmesolve notebook by P.D. Nation.
Introduction¶
The Bloch-Redfield solver is another method to solve a master equation. In comparison to the Lindblad Master equation solver qutip.mesolve() the Bloch-Redfield solver qutip.brmesolve() differs in the description of the interaction with the environment. In qutip.mesolve() we described the dissipation by collapse operators, which do not necessarily have a physical interpretation. The qutip.brmesolve() function requires the a dissipation description by the so-called noise-power-spectrum, which gives the intensity of the dissipation depending on the frequency $\omega$.
In this notebook we will introduce the basic usage of qutip.brmesolve() and compare it to qutip.mesolve(). For more information on the Bloch-Redfield solver see the follow-up notebooks and the QuTiP Documentation of the functionality.
Imports¶
import matplotlib.pyplot as plt
import numpy as np
from qutip import (about, basis, bloch_redfield_tensor, brmesolve, expect,
hinton, liouvillian, mesolve, plot_expectation_values,
sigmam, sigmax, sigmay, sigmaz, steadystate)
%matplotlib inline
Two-level system evolution¶
In this example we consider a simple two level system described by the Hamiltonian:
$$ H = \frac{\epsilon}{2} \sigma_z$$Furthermore, we define a constant dissipation rate to the environment $\gamma$.
epsilon = 0.5 * 2 * np.pi
gamma = 0.25
times = np.linspace(0, 10, 100)
Setup the Hamiltonian, initial state and collapse operators for the qutip.mesolve() function. We choose a superposition of states as initial state and want to observe the expectation values of $\sigma_x, \sigma_y, \sigma_z$.
# Setup Hamiltonian and initial state
H = epsilon / 2 * sigmaz()
psi0 = (2 * basis(2, 0) + basis(2, 1)).unit()
# Setup the master equation solver
c_ops = [np.sqrt(gamma) * sigmam()]
e_ops = [sigmax(), sigmay(), sigmaz()]
result_me = mesolve(H, psi0, times, c_ops, e_ops)
For the qutip.brmesolve function we have to give the interaction of the system with the bath as a hermitian operator together with a noise power spectrum, which defines the strength of the interaction per frequency. Here we define a constant interaction whenever the frequency is positive and no dissipation for negative frequencies. This allows us to use sigmax() ( a hermitian operator) instead of the non-hermitian operator sigmam used above.
The usage of hermitian operators simplifies the internal numerical implementation and leads to vanishing cross-correlations between different environment operators (if multiple are given).
a_op = [sigmax(), lambda w: gamma * (w > 0.0)]
Instead of the c_ops we now pass the a_ops to the Bloch-Redfield solver.
result_brme = brmesolve(H, psi0, times, [a_op], e_ops)
We can now compare the expectation values for every operator we passed to the solvers in e_ops. As expected both solvers, mesolve and brmesolve, produce similar results.
plot_expectation_values(
[result_me, result_brme], ylabels=["<X>", "<Y>", "<Z>"], show_legend=True
);
Storing States instead of expectation values¶
As for the other solvers provided in QuTiP, we can obtain the density matrices at each defined time step instead of some expectation values. To do so, we pass an empty list as e_ops argument. If you want to calculate expectation values (i.e. non-empty e_ops) and obtain the states at the same time you can also pass options={"store_states": True} to the solver functions.
# run solvers without e_ops
me_s = mesolve(H, psi0, times, c_ops, e_ops=[])
brme_s = brmesolve(H, psi0, times, [a_op], e_ops=[])
# calculate expecation values
x_me = expect(sigmax(), me_s.states)
x_brme = expect(sigmax(), brme_s.states)
# plot the expectation values
plt.plot(times, x_me, label="ME")
plt.plot(times, x_brme, label="BRME")
plt.legend(), plt.xlabel("time"), plt.ylabel("<X>");
Bloch-Redfield Tensor¶
We described the dynmamics of the system by the Bloch-Redfield master equation, which is constructed from the Bloch-Redfield tensor $R_{abcd}$ (see documentation of Bloch-Redfield master equation). Hence the dynamics are determined by this tensor. We can calculate the tensor in QuTiP using the qutip.bloch_redfield_tensor() function. We have to pass the Hamiltonian of the system and the dissipation description in a_ops to construct $R_{abcd}$. Furthermore, the function gives us the eigenstates of the Hamiltonian, as they are calculated along the way.
R, H_ekets = bloch_redfield_tensor(H, [a_op])
# calculate lindblad liouvillian from H
L = liouvillian(H, c_ops)
We can now use the Bloch-Redfield Tensor and the Lindblad Liouvillain to calculate the steadystate for both approaches. As we saw above the dynamics were the same for using the different solvers, hence we expect the steadystate to be equal too. We use the qutip.hinton() function to plot the steadystate density matrix for both approaches and can see that they are the same.
We have to transform the steadystate density matrix we obtain from the Bloch-Redfield tensor using the eigenstates of the Hamiltonian, as R is expressed in the eigenbasis of H.
# Obtain steadystate from Bloch-Redfield Tensor
rhoss_br_eigenbasis = steadystate(R)
rhoss_br = rhoss_br_eigenbasis.transform(H_ekets, True)
# Steadystate from Lindblad liouvillian
rhoss_me = steadystate(L)
# Plot the density matrices using a hinton plot
hinton(rhoss_br, title="Bloch-Redfield steadystate")
hinton(rhoss_me, title="Lindblad-ME steadystate");
About¶
about()
QuTiP: Quantum Toolbox in Python ================================ Copyright (c) QuTiP team 2011 and later. Current admin team: Alexander Pitchford, Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, Eric Giguère, Boxi Li, Jake Lishman, Simon Cross and Asier Galicia. Board members: Daniel Burgarth, Robert Johansson, Anton F. Kockum, Franco Nori and Will Zeng. Original developers: R. J. Johansson & P. D. Nation. Previous lead developers: Chris Granade & A. Grimsmo. Currently developed through wide collaboration. See https://github.com/qutip for details. QuTiP Version: 4.7.1 Numpy Version: 1.22.4 Scipy Version: 1.8.1 Cython Version: 0.29.33 Matplotlib Version: 3.5.2 Python Version: 3.10.4 Number of CPUs: 2 BLAS Info: Generic OPENMP Installed: False INTEL MKL Ext: False Platform Info: Linux (x86_64) Installation path: /home/runner/work/qutip-tutorials/qutip-tutorials/qutip/qutip ================================================================================ Please cite QuTiP in your publication. ================================================================================ For your convenience a bibtex reference can be easily generated using `qutip.cite()`
Testing¶
# Verify that mesolve and brmesolve generate similar results
assert np.allclose(result_me.expect[0], result_brme.expect[0])
assert np.allclose(result_me.expect[1], result_brme.expect[1])
assert np.allclose(result_me.expect[2], result_brme.expect[2])
assert np.allclose(x_me, x_brme)
# assume steadystate is the same
assert np.allclose(rhoss_br, rhoss_me)