Skip to content

Commit

Permalink
Merge pull request #2271 from Ericgig/doc.dynamics
Browse files Browse the repository at this point in the history 
Improve dynamics section of docs
  • Loading branch information
@Ericgig
Ericgig committed Dec 19, 2023
2 parents a8d041b + 20c6983 commit 7fbb567
Show file tree
Hide file tree
Showing 23 changed files with 1,374 additions and 654 deletions.
1 change: 1 addition & 0 deletions doc/changes/2271.misc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Improve documentation in guide/dynamics
103 changes: 56 additions & 47 deletions doc/guide/dynamics/dynamics-bloch-redfield.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,16 +87,19 @@ This allows us to write master equation in terms of system operators and bath co
g_{\alpha\beta}(-\tau) \left[\rho_S(t)A_\alpha(t-\tau)A_\beta(t) - A_\alpha(t)\rho_S(t)A_\beta(t-\tau)\right]
\right\},
where :math:`g_{\alpha\beta}(\tau) = {\rm Tr}_B\left[B_\alpha(t)B_\beta(t-\tau)\rho_B\right] = \left<B_\alpha(\tau)B_\beta(0)\right>`, since the bath state :math:`\rho_B` is a steady state.
where :math:`g_{\alpha\beta}(\tau) = {\rm Tr}_B\left[B_\alpha(t)B_\beta(t-\tau)\rho_B\right] = \left<B_\alpha(\tau)B_\beta(0)\right>`,
since the bath state :math:`\rho_B` is a steady state.

In the eigenbasis of the system Hamiltonian, where :math:`A_{mn}(t) = A_{mn} e^{i\omega_{mn}t}`, :math:`\omega_{mn} = \omega_m - \omega_n` and :math:`\omega_m` are the eigenfrequencies corresponding the eigenstate :math:`\left|m\right>`, we obtain in matrix form in the Schrödinger picture
In the eigenbasis of the system Hamiltonian, where :math:`A_{mn}(t) = A_{mn} e^{i\omega_{mn}t}`,
:math:`\omega_{mn} = \omega_m - \omega_n` and :math:`\omega_m` are the eigenfrequencies
corresponding the eigenstate :math:`\left|m\right>`, we obtain in matrix form in the Schrödinger picture

.. math::
\frac{d}{dt}\rho_{ab}(t)
=
-i\omega_{ab}\rho_{ab}(t)
-\hbar^{-2}
=&
-i\omega_{ab}\rho_{ab}(t) \nonumber\\
&-\hbar^{-2}
\sum_{\alpha,\beta}
\sum_{c,d}^{\rm sec}
\int_0^\infty d\tau\;
Expand All @@ -107,7 +110,7 @@ In the eigenbasis of the system Hamiltonian, where :math:`A_{mn}(t) = A_{mn} e^{
A^\alpha_{ac} A^\beta_{db} e^{i\omega_{ca}\tau}
\right]
\right. \nonumber\\
+
&+
\left.
g_{\alpha\beta}(-\tau)
\left[\delta_{ac}\sum_n A^\alpha_{dn}A^\beta_{nb} e^{i\omega_{nd}\tau}
Expand Down Expand Up @@ -185,8 +188,9 @@ Bloch-Redfield master equation in QuTiP



In QuTiP, the Bloch-Redfield tensor Eq. :eq:`br-tensor` can be calculated using the function :func:`qutip.bloch_redfield.bloch_redfield_tensor`.
It takes two mandatory arguments: The system Hamiltonian :math:`H`, a nested list of operator :math:`A_\alpha`, spectral density functions :math:`S_\alpha(\omega)` pairs that characterize the coupling between system and bath.
In QuTiP, the Bloch-Redfield tensor Eq. :eq:`br-tensor` can be calculated using the function :func:`.bloch_redfield_tensor`.
It takes two mandatory arguments: The system Hamiltonian :math:`H`, a nested list of operator
:math:`A_\alpha`, spectral density functions :math:`S_\alpha(\omega)` pairs that characterize the coupling between system and bath.
The spectral density functions are Python callback functions that takes the (angular) frequency as a single argument.

To illustrate how to calculate the Bloch-Redfield tensor, let's consider a two-level atom
Expand Down Expand Up @@ -231,13 +235,22 @@ To illustrate how to calculate the Bloch-Redfield tensor, let's consider a two-l
[ 0. +0.j 0. +0.j 0. +0.j
-0.24514517+0.j ]]

Note that it is also possible to add Lindblad dissipation superoperators in the Bloch-Refield tensor by passing the operators via the ``c_ops`` keyword argument like you would in the :func:`qutip.mesolve` or :func:`qutip.mcsolve` functions.
For convenience, the function :func:`qutip.bloch_redfield_tensor` also returns the basis transformation operator, the eigen vector matrix, since they are calculated in the process of calculating the Bloch-Redfield tensor `R`, and the `ekets` are usually needed again later when transforming operators between the laboratory basis and the eigen basis.
The tensor can be obtained in the laboratory basis by setting ``fock_basis=True``, in that case, the transformation operator is not returned.
Note that it is also possible to add Lindblad dissipation superoperators in the
Bloch-Refield tensor by passing the operators via the ``c_ops`` keyword argument
like you would in the :func:`.mesolve` or :func:`.mcsolve` functions.
For convenience, the function :func:`.bloch_redfield_tensor` also returns the basis
transformation operator, the eigen vector matrix, since they are calculated in the
process of calculating the Bloch-Redfield tensor `R`, and the `ekets` are usually
needed again later when transforming operators between the laboratory basis and the eigen basis.
The tensor can be obtained in the laboratory basis by setting ``fock_basis=True``,
in that case, the transformation operator is not returned.


The evolution of a wavefunction or density matrix, according to the Bloch-Redfield master equation :eq:`br-final`, can be calculated using the QuTiP function :func:`qutip.mesolve` using Bloch-Refield tensor in the laboratory basis instead of a liouvillian.
For example, to evaluate the expectation values of the :math:`\sigma_x`, :math:`\sigma_y`, and :math:`\sigma_z` operators for the example above, we can use the following code:
The evolution of a wavefunction or density matrix, according to the Bloch-Redfield
master equation :eq:`br-final`, can be calculated using the QuTiP function :func:`.mesolve`
using Bloch-Refield tensor in the laboratory basis instead of a liouvillian.
For example, to evaluate the expectation values of the :math:`\sigma_x`,
:math:`\sigma_y`, and :math:`\sigma_z` operators for the example above, we can use the following code:

.. plot::
:context:
Expand Down Expand Up @@ -274,24 +287,35 @@ For example, to evaluate the expectation values of the :math:`\sigma_x`, :math:`

sphere.make_sphere()

The two steps of calculating the Bloch-Redfield tensor and evolving according to the corresponding master equation can be combined into one by using the function :func:`qutip.brmesolve`, which takes same arguments as :func:`qutip.mesolve` and :func:`qutip.mcsolve`, save for the additional nested list of operator-spectrum pairs that is called ``a_ops``.
The two steps of calculating the Bloch-Redfield tensor and evolving according to
the corresponding master equation can be combined into one by using the function
:func:`.brmesolve`, which takes same arguments as :func:`.mesolve` and
:func:`.mcsolve`, save for the additional nested list of operator-spectrum
pairs that is called ``a_ops``.

.. plot::
:context: close-figs

output = brmesolve(H, psi0, tlist, a_ops=[[sigmax(),ohmic_spectrum]], e_ops=e_ops)

where the resulting `output` is an instance of the class :class:`qutip.Result`.
where the resulting `output` is an instance of the class :class:`.Result`.

.. note::
While the code example simulates the Bloch-Redfield equation in the secular approximation, QuTiP's implementation allows the user to simulate the non-secular version of the Bloch-Redfield equation by setting ``sec_cutoff=-1``, as well as do a partial secular approximation by setting it to a ``float`` , this float will become the cutoff for the sum in :eq:`br-final` meaning terms with :math:`|\omega_{ab}-\omega_{cd}|` greater than the cutoff will be neglected.
While the code example simulates the Bloch-Redfield equation in the secular
approximation, QuTiP's implementation allows the user to simulate the non-secular
version of the Bloch-Redfield equation by setting ``sec_cutoff=-1``, as well as
do a partial secular approximation by setting it to a ``float`` , this float
will become the cutoff for the sum in :eq:`br-final` meaning terms with
:math:`|\omega_{ab}-\omega_{cd}|` greater than the cutoff will be neglected.
Its default value is 0.1 which corresponds to the secular approximation.
For example the command
::

output = brmesolve(H, psi0, tlist, a_ops=[[sigmax(), ohmic_spectrum]], e_ops=e_ops, sec_cutoff=-1)
output = brmesolve(H, psi0, tlist, a_ops=[[sigmax(), ohmic_spectrum]],
e_ops=e_ops, sec_cutoff=-1)

will simulate the same example as above without the secular approximation. Note that using the non-secular version may lead to negativity issues.
will simulate the same example as above without the secular approximation.
Note that using the non-secular version may lead to negativity issues.

.. _td-bloch-redfield:

Expand All @@ -300,17 +324,23 @@ Time-dependent Bloch-Redfield Dynamics

If you have not done so already, please read the section: :ref:`time`.

As we have already discussed, the Bloch-Redfield master equation requires transforming into the eigenbasis of the system Hamiltonian.
As we have already discussed, the Bloch-Redfield master equation requires transforming
into the eigenbasis of the system Hamiltonian.
For time-independent systems, this transformation need only be done once.
However, for time-dependent systems, one must move to the instantaneous eigenbasis at each time-step in the evolution, thus greatly increasing the computational complexity of the dynamics.
However, for time-dependent systems, one must move to the instantaneous eigenbasis
at each time-step in the evolution, thus greatly increasing the computational complexity of the dynamics.
In addition, the requirement for computing all the eigenvalues severely limits the scalability of the method.
Fortunately, this eigen decomposition occurs at the Hamiltonian level, as opposed to the super-operator level, and thus, with efficient programming, one can tackle many systems that are commonly encountered.
Fortunately, this eigen decomposition occurs at the Hamiltonian level, as opposed to the
super-operator level, and thus, with efficient programming, one can tackle many systems that are commonly encountered.


For time-dependent Hamiltonians, the Hamiltonian itself can be passed into the solver like any other time dependent Hamiltonian, as thus we will not discuss this topic further.
For time-dependent Hamiltonians, the Hamiltonian itself can be passed into the solver
like any other time dependent Hamiltonian, as thus we will not discuss this topic further.
Instead, here the focus is on time-dependent bath coupling terms.
To this end, suppose that we have a dissipative harmonic oscillator, where the white-noise dissipation rate decreases exponentially with time :math:`\kappa(t) = \kappa(0)\exp(-t)`.
In the Lindblad or Monte Carlo solvers, this could be implemented as a time-dependent collapse operator list ``c_ops = [[a, 'sqrt(kappa*exp(-t))']]``.
To this end, suppose that we have a dissipative harmonic oscillator, where the white-noise
dissipation rate decreases exponentially with time :math:`\kappa(t) = \kappa(0)\exp(-t)`.
In the Lindblad or Monte Carlo solvers, this could be implemented as a time-dependent
collapse operator list ``c_ops = [[a, 'sqrt(kappa*exp(-t))']]``.
In the Bloch-Redfield solver, the bath coupling terms must be Hermitian.
As such, in this example, our coupling operator is the position operator ``a+a.dag()``.
The complete example, and comparison to the analytic expression is:
Expand All @@ -320,31 +350,21 @@ The complete example, and comparison to the analytic expression is:
:context: close-figs

N = 10 # number of basis states to consider

a = destroy(N)

H = a.dag() * a

psi0 = basis(N, 9) # initial state

kappa = 0.2 # coupling to oscillator

a_ops = [
([a+a.dag(), f'sqrt({kappa}*exp(-t))'], '(w>=0)')
]

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

out = brmesolve(H, psi0, tlist, a_ops, e_ops=[a.dag() * a])

actual_answer = 9.0 * np.exp(-kappa * (1.0 - np.exp(-tlist)))

plt.figure()

plt.plot(tlist, out.expect[0])

plt.plot(tlist, actual_answer)

plt.show()


Expand All @@ -361,45 +381,34 @@ In this example, the ``a_ops`` list would be:
]


where the first tuple element ``[[a, 'exp(1j*t)'], [a.dag(), 'exp(-1j*t)']]`` tells the solver what is the time-dependent Hermitian coupling operator.
where the first tuple element ``[[a, 'exp(1j*t)'], [a.dag(), 'exp(-1j*t)']]`` tells
the solver what is the time-dependent Hermitian coupling operator.
The second tuple ``f'{kappa} * (w >= 0)'``, gives the noise power spectrum.
A full example is:

.. plot::
:context: close-figs

N = 10

w0 = 1.0 * 2 * np.pi

g = 0.05 * w0

kappa = 0.15

times = np.linspace(0, 25, 1000)

a = destroy(N)

H = w0 * a.dag() * a + g * (a + a.dag())

psi0 = ket2dm((basis(N, 4) + basis(N, 2) + basis(N, 0)).unit())

a_ops = [[
QobjEvo([[a, 'exp(1j*t)'], [a.dag(), 'exp(-1j*t)']]), (f'{kappa} * (w >= 0)')
]]

e_ops = [a.dag() * a, a + a.dag()]

res_brme = brmesolve(H, psi0, times, a_ops, e_ops)

plt.figure()

plt.plot(times, res_brme.expect[0], label=r'$a^{+}a$')

plt.plot(times, res_brme.expect[1], label=r'$a+a^{+}$')

plt.legend()

plt.show()

Further examples on time-dependent Bloch-Redfield simulations can be found in the online tutorials.
Expand Down
166 changes: 166 additions & 0 deletions doc/guide/dynamics/dynamics-class.rst
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,166 @@
.. _solver_class:



*******************************************
Solver Class Interface
*******************************************

In QuTiP version 5 and later, solvers such as :func:`.mesolve`, :func:`.mcsolve` also have
a class interface. The class interface allows reusing the Hamiltonian and fine tuning
many details of how the solver is run.

Examples of some of the solver class features are given below.

Reusing Hamiltonian Data
------------------------

There are many cases where one would like to study multiple evolutions of
the same quantum system, whether by changing the initial state or other parameters.
In order to evolve a given system as fast as possible, the solvers in QuTiP
take the given input operators (Hamiltonian, collapse operators, etc) and prepare
them for use with the selected ODE solver.

These operations are usually reasonably fast, but for some solvers, such as
:func:`.brmesolve` or :func:`.fmmesolve`, the overhead can be significant.
Even for simpler solvers, the time spent organizing data can become appreciable
when repeatedly solving a system.

The class interface allows us to setup the system once and reuse it with various
parameters. Most ``...solve`` function have a paired ``...Solver`` class, with a
``..Solver.run`` method to run the evolution. At class
instance creation, the physics (``H``, ``c_ops``, ``a_ops``, etc.) and options
are passed. The initial state, times and expectation operators are only passed
when calling ``run``:

.. plot::
:context: close-figs

times = np.linspace(0.0, 6.0, 601)
a = tensor(qeye(2), destroy(10))
sm = tensor(destroy(2), qeye(10))
e_ops = [a.dag() * a, sm.dag() * sm]
H = QobjEvo(
[a.dag()*a + sm.dag()*sm, [(sm*a.dag() + sm.dag()*a), lambda t, A: A]],
args={"A": 0.5*np.pi}
)

solver = MESolver(H, c_ops=[np.sqrt(0.1) * a], options={"atol": 1e-8})
solver.options["normalize_output"] = True
psi0 = tensor(fock(2, 0), fock(10, 5))
data1 = solver.run(psi0, times, e_ops=e_ops)
psi1 = tensor(fock(2, 0), coherent(10, 2 - 1j))
data2 = solver.run(psi1, times, e_ops=e_ops)

plt.figure()
plt.plot(times, data1.expect[0], "b", times, data1.expect[1], "r", lw=2)
plt.plot(times, data2.expect[0], 'b--', times, data2.expect[1], 'r--', lw=2)
plt.title('Master Equation time evolution')
plt.xlabel('Time', fontsize=14)
plt.ylabel('Expectation values', fontsize=14)
plt.legend(("cavity photon number", "atom excitation probability"))
plt.show()


Note that as shown, options can be set at initialization or with the
``options`` property.

The simulation parameters, the ``args`` of the :class:`.QobjEvo` passed as system
operators, can be updated at the start of a run:

.. plot::
:context: close-figs

data1 = solver.run(psi0, times, e_ops=e_ops)
data2 = solver.run(psi0, times, e_ops=e_ops, args={"A": 0.25*np.pi})
data3 = solver.run(psi0, times, e_ops=e_ops, args={"A": 0.125*np.pi})

plt.figure()
plt.plot(times, data1.expect[0], label="A=pi/2")
plt.plot(times, data2.expect[0], label="A=pi/4")
plt.plot(times, data3.expect[0], label="A=pi/8")
plt.title('Master Equation time evolution')
plt.xlabel('Time', fontsize=14)
plt.ylabel('Expectation values', fontsize=14)
plt.legend()
plt.show()


Stepping through the run
------------------------

The solver class also allows to run through a simulation one step at a time, updating
args at each step:


.. plot::
:context: close-figs

data = [5.]
solver.start(state0=psi0, t0=times[0])
for t in times[1:]:
psi_t = solver.step(t, args={"A": np.pi*np.exp(-(t-3)**2)})
data.append(expect(e_ops[0], psi_t))

plt.figure()
plt.plot(times, data)
plt.title('Master Equation time evolution')
plt.xlabel('Time', fontsize=14)
plt.ylabel('Expectation values', fontsize=14)
plt.legend(("cavity photon number"))
plt.show()


.. note::

This is an example only, updating a constant ``args`` parameter between step
should not replace using a function as QobjEvo's coefficient.

.. note::

It is possible to create multiple solvers and to advance them using ``step`` in
parallel. However, many ODE solver, including the default ``adams`` method, only
allow one instance at a time per process. QuTiP supports using multiple solver instances
of these ODE solvers but with a performance cost. In these situations, using
``dop853`` or ``vern9`` integration method is recommended instead.




Feedback: Accessing the solver state from evolution operators
=============================================================

The state of the system during the evolution is accessible via properties of the solver classes.

Each solver has a ``StateFeedback`` and ``ExpectFeedback`` class method that can
be passed as arguments to time dependent systems. For example, ``ExpectFeedback``
can be used to create a system which uncouples when there are 5 or fewer photons in the
cavity.

.. plot::
:context: close-figs

def f(t, e1):
ex = (e1.real - 5)
return (ex > 0) * ex * 10

times = np.linspace(0.0, 1.0, 301)
a = tensor(qeye(2), destroy(10))
sm = tensor(destroy(2), qeye(10))
e_ops = [a.dag() * a, sm.dag() * sm]
psi0 = tensor(fock(2, 0), fock(10, 8))
e_ops = [a.dag() * a, sm.dag() * sm]

H = [a*a.dag(), [sm*a.dag() + sm.dag()*a, f]]
data = mesolve(H, psi0, times, c_ops=[a], e_ops=e_ops,
args={"e1": MESolver.ExpectFeedback(a.dag() * a)}
).expect

plt.figure()
plt.plot(times, data[0])
plt.plot(times, data[1])
plt.title('Master Equation time evolution')
plt.xlabel('Time', fontsize=14)
plt.ylabel('Expectation values', fontsize=14)
plt.legend(("cavity photon number", "atom excitation probability"))
plt.show()

0 comments on commit 7fbb567

Please sign in to comment.