QuTiPv5 Paper Example: sesolve and mesolve and the new solver class¶

Authors: Maximilian Meyer-Mölleringhof (m.meyermoelleringhof@gmail.com), Neill Lambert (nwlambert@gmail.com), Paul Menczel (paul@menczel.net)

In QuTiP, the Qobj and QobjEvo classes form the very heart of almost all calculations that can be performed. With these classes, open quantum systems with various interactions and structures can be simulated by using the wide variety of solvers provided. Most of the time, these solvers are given an initial state, a Hamiltonian and an environment that is often described using rates or coupling strengths. QuTiP then uses numerical integration to determine the time evolution of the system.

QuTiP v5 [1] introduces a unified interface for interacting with these solvers. This new interface is class-based, allowing users to instantiate a solver object for a specific problem. This can be useful when the same Hamiltonian data is reused with different initial conditions, time steps or other options. A noticable speed-up can thus be achieved if solvers are reused multiple times.

Upon instantiation, one first supplies only the Hamiltonian and the collapse operators (e.g., collapse operators for a Lindbladian master equation). Initial conditions, time steps, etc. are passed to the Solver.run() method which then performs the simulation.

In this notebook we will consider several examples illustrating the usage of the new solver classes.

Part 0: Introduction to the New Solver Class¶

In our first example, we want to look at two interacting qubits that are (for now) decoupled from an environment. Such a system is described by the Hamiltonian

$$H = \dfrac{\epsilon_1}{2} \sigma_z^{(1)} + \dfrac{\epsilon_2}{2} \sigma_z^{(2)} + g \sigma_{x}^{(1)} \sigma_{x}^{(2)}.$$

The Pauli matrices $\sigma_z^{(1/2)}$ describe the respective two-level system of the qubits. The qubits are coupled via $\sigma_{x}^{(1/2)}$ and their interaction strength is given by $g$.

The dynamics of such a system is described by the Schrödinger equation

$$i \hbar \dfrac{d}{dt} \ket{\psi} = H \ket{\psi}.$$

Therefore, we can use SESolver to calculate the dynamics.

Solver and Integrator Options¶

Another change in QuTiP v5 is that the options argument takes a standard Python dictionary. This should increase future flexibility and allow different solvers to provide individual sets of options more easily. The complete list of options can be found in the online documentation for each solver.

As an example of frequently used options, we show store_states, determining whether the output should include the system state at each time step, and store_final_state, determining the same but for the final state. method is another common option, specifying the ODE integration method. Specific options related to the numeric integration depend on the chosen method; here, we show atol controlling the precision (absolute tolerance), nsteps controlling the maximum number of steps between two time steps, and max_step giving the maximum allowed integration step of the default Adams ODE integration method.

Part 1: Lindblad Dynamics and Beyond¶

In general, the Schrödinger equation describes the dynamics of any quantum system. Once systems become large or we consider continuous systems, however, solving it often becomes impossible. Therefore, master equations of various types were developed and have now become the most common way to describe the dynamics of finite (open) quantum systems. Generally, a master equation refers to a first-order linear differential equation for $\rho(t)$ which is the reduced density operator describing the quantum state. In QuTiP, mesolve is the general solver we use to solve such master equations. Although it supports master equations of various forms, the Lindbladian type is implemented by default. The general form of such an equation is given by

$$ \dot{\rho}(t) = - \dfrac{i}{\hbar} [H(t), \rho(t)] + \sum_n \dfrac{1}{2}[ 2 C_n \rho(t) C_n^\dagger - \rho(t) C_n^\dagger C_n - C^\dagger_n C_n \rho(t) ] . $$

Next to the density operator $\rho(t)$ and the Hamiltonian $H(t)$, this equation includes the so-called collapse (or jump) operators $C_n = \sqrt{\gamma_n} A_n$. They define the dissipation due to contact with an environment. The rates $\gamma_n$ describe the frequency of transitions between the states connected by the operator $A_n$.

To continue our example of the two qubits, we now connect them to an evironment using the collapse operators $C_1 = \sqrt{\gamma} \sigma_{-}^{(1)}$ and $C_2 = \sqrt{\gamma} \sigma_{-}^{(2)}$ where $\sigma_{-}^{(i)}$ takes qubit (i) from its excited state to its ground state. This time, however, we will be using the mesolve solver.

Global Master Equation - Born-Markov-secular approximation¶

In the previous example, the collapse operators acted locally on each qubit. However, depending on the considered approximation, different kinds of collapse operators are found. One example for this is when the qubits interact more strongly with each other than with the bath. One then arrives at the global master equation under the standard Born-Markox approximation. Although the collapse operators still act like annihilation and creation operators here, they now act on the total coupled eigenstates of the interacting two-qubit system

$$A_{ij} = \ket{\psi_i}\bra{\psi_i}$$

with rates

$$\gamma_{ij} = | \bra{\psi_i} d \ket{\psi_j} |^2 S(\Delta_{ij}).$$

The $\ket{\psi_i}$ are the eigenstates of $H$ and $\Delta_{ij} = E_j - E_i$ are the differences of eigenenergies. $d$ is the coupling operator of the system to the environment. We use the power spectrum

$$S(\omega) = 2 J(\omega) [n_{th} (\omega) + 1] \theta(\omega) + 2J(-\omega)[n_{th}(-\omega)]\theta(-\omega)$$

which depends on details of the environment such as its spectral density $J(\omega)$ and, through the Bose-Einstein distribution $n_{th} (\omega)$, its temperature. $\theta$ denotes the Heaviside function.

By assuming spectral density to be flat, $J(\omega) = \gamma / 2$, and considering zero temperature, we can write

$$S(\omega) = \gamma \theta(\omega).$$

For this example, we manually implement this zero temperature environment for our two-qubit system using mesolve().

It is interesting to note that the long-time evolution leads to a state that is close to the coupled ground state of the two qubit system:

Solver comparison¶

In the following, we compare the results of the local and dressed (global) Lindblad simulations from above with the Bloch-Redfield solver. The Bloch-Redfield solver is explained in more detail in other tutorials, but we use it here to solve the weak-coupling master equation from a given bath power spectrum. For small coupling strengths, the results from the local and global master equations both agree with the Bloch-Redfield solver. However, this changes when considering stronger couplings where the local master equation deviates from the result of the global and the Bloch-Redfield approach.

Manual Liouvillian Superoperator¶

As mentioned before, QuTiP's master equation solver can be used to solve any other master equation. By using spre(), spost() and sprepost(), we can manually construct such equations. These functions specifically convert the operators from the original Hilbert space to operators in the double space which is internally used by QuTiP to optimize computations (more details on this can be found in the QuTiPv5 paper).

For example, the Lindbladian corresponding to the master equation of the previous example can be constructed manually via:

Part 2: Time-Dependent Systems¶

Finally, we compare the results from another example, that of a driven system using, first, mesolve() with and without the rotating-wave approximation (RWA) for the drive, and then without RWA with the Bloch-Redfield and HEOM solvers. For this, we consider a Hamiltonian of the form

$$ H = \frac{\Delta}{2} \sigma_z + \frac{A}{2} \sin (\omega_d t) \sigma_x , $$

where $\Delta$ is the energy splitting, $A$ is the drive amplitude and $\omega_d$ is the drive's frequency.

With mesolve¶

With Rotating-Wave-Approximated mesolve¶

With brmesolve¶

With HEOMSolver¶

Comparison of Solver Results¶

Adiabatic Energy Switching¶

To show where usage of these local-basis collapse operators can lead to errors, we look at a single qubit system. This time however, the energies of the qubit are adiabatically switched between positive and negative values.

$$H = \dfrac{\Delta}{2} \sin{(\omega_d t)} \sigma_z.$$

If the drive is slow enough, we expect the bath to respond to this change and transitions from higher to lower energy levels should be induced. These changes in system energy thus have to be taken into account so that mesolve() still produces a correct result.

References¶

[1] QuTiP 5: The Quantum Toolbox in Python

About¶

Testing¶