QuTiP 5: The Quantum Toolbox in Python

TL;DR

QuTiP v5 tackles the exponential growth of open-quantum-system simulations by introducing a flexible, multi-format data layer that interplays with CSR, Dense, Dia, and JAX-based formats, enabling GPU acceleration and automatic differentiation. A unified solver class interface couples with a broad suite of solvers (mesolve, mcsolve, nm_mcsolve, brmesolve, HEOM, smesolve) and new time-dependent capabilities, including Diffrax/JAX for GPU performance and automatic differentiation for control and statistics. The release also expands modularity with sub-packages QuTiP-QOC and QuTiP-QIP, enhances visualization, and introduces ENR states and MPI/HPC support, boosting scalability for large-scale quantum circuits and many-body dynamics. Collectively, these innovations position QuTiP as a robust, open-source platform for research, teaching, and industrial quantum computing development, with significant implications for simulating QPUs, open-system dynamics, and quantum control across diverse hardware backends.

Abstract

QuTiP, the Quantum Toolbox in Python, has been at the forefront of open-source quantum software for the past 13 years. It is used as a research, teaching, and industrial tool, and has been downloaded millions of times by users around the world. Here we introduce the latest developments in QuTiP v5, which are set to have a large impact on the future of QuTiP and enable it to be a modern, continuously developed and popular tool for another decade and more. We summarize the code design and fundamental data layer changes as well as efficiency improvements, new solvers, applications to quantum circuits with QuTiP-QIP, and new quantum control tools with QuTiP-QOC. Additional flexibility in the data layer underlying all ``quantum objects'' in QuTiP allows us to harness the power of state-of-the-art data formats and packages like JAX, CuPy, and more. We explain these new features with a series of both well-known and new examples. The code for these examples is available in a static form on GitHub and as continuously updated and documented notebooks in the qutip-tutorials package.

Figures (26)

• Figure 1: A schematic overview of the QuTiP project, describing Qobj and its features/functions, solvers, and QuTiP sub-packages. For a complete list of Qobj methods and attributes see Table \ref{['tab:qobj_methods']}, for a list of libraries which QuTiP uses see Table \ref{['tab:glossary_acronyms']}, for a glossary of terms commonly used in QuTiP see Table \ref{['tab:glossary_terms']}, and for a list of state, operator, superoperator, entanglement measure and metric functions see Tables \ref{['tab:qutip_states_functions']}, \ref{['tab:qutip_operators_functions']}, \ref{['tab:qutip_superoperators_functions']}, \ref{['table:entropy_functions']} and \ref{['table:metrics_functions']}.
• Figure 2: For the example problem of two interacting qubits, we compare the output of two mesolve() simulations. One uses local collapse operators acting on the bare states of each qubit (local Lindblad equation), and the other uses dressed collapse operators acting on the global eigenstates (global Lindblad equation). In addition, we include a simulation performed with the Bloch-Redfield solver (brmesolve()). Panel (a) shows weak coupling between the qubits ($g=0.1\epsilon_1$), where the local Lindblad description is sufficient from the perspective of detailed balance. Panel (b) shows the dynamics with strong coupling ($g=2\epsilon_1$), where local collapse operators predict a steady state that is not compatible with detailed balance. In both cases, the qubits are resonant $\epsilon_1=\epsilon_2$.
• Figure 3: In panel (a), we show the dynamics of a resonantly driven qubit simulated with mesolve, brmesolve and the HEOM solver. We also compare the results to a simulation of the time-independent RWA Hamiltonian, Eq. \ref{['eq:mesolve:H_RWA']}. The bath is chosen to have zero temperature and a flat spectral density $J(\omega) = \gamma$ with a rate $\gamma = 0.005 \Delta /(2\pi)$. The drive is assumed to be weak ($A=0.01 \Delta$) and on resonance ($\omega_{d} = \Delta$). In panel (b), we show the second example, where the energy of a qubit is adiabatically, sinusoidally modulated between positive and negative values with frequency $\omega_d = 0.05 \Delta$ and amplitude $A = \Delta$. The qubit is in contact with a flat zero-temperature bath with rate $\gamma = 0.05 \Delta / (2\pi)$. The naive mesolve implementation with a single collapse operator that is not sensitive to the system's energy predicts just exponential decay. The Bloch-Redfield and HEOM solvers predict reversal of the decay when the energy changes sign, as a true zero-temperature bath can only remove energy from the system. The heomsolve result differs slightly from the brmesolve result since the spectral density used in heomsolve is not flat. When the same spectral density used in the HEOM method is used explicitly in the Bloch-Redfield solver both approaches agree. This is shown with the green dotted line which we called brmesolve non-flat, which uses the same spectral density as HEOM.
• Figure 4: Benchmark of the time required to solve the dynamics of an Ising spin chain as a function of the number of spins $N$ using the JAX data layer. The simulation was performed on an NVIDIA A100 GPU with 80 GB of RAM. The CPU used for comparison was an AMD EPYC 7713 (64 cores). Panel (a) shows a noiseless example for up to $22$ spins ($24$ on CPU) using sesolve. Panel (b) shows the same problem in the presence of noise for up to $11$ spins ($12$ on CPU) using mesolve. The Schrödinger equation solver is able to integrate more spins because of the memory cost of the superoperator constructed by mesolve. In both cases, we see a crossover at a certain system size where the GPU solver becomes more performant, but the GPU memory limit is reached shortly after. For the examples run on CPU we also differentiate between the default SciPy based solver and using JAX in CPU mode. In all JAX examples we used the jaxdia data-layer. The dense Jax format runs out of memory much more quickly (not shown).
• Figure 5: Here we show the same example as Fig. \ref{['fig1']} evaluated with the Monte-Carlo solver mcsolve. The average behavior is compared to the result from mesolve, alongside two example trajectories showing discrete quantum jumps. The yellow shaded region indicates the convergence error based on the standard deviation of the trajectories and the number of trajectories, $\sigma_{\mathrm{err}} = \sigma / \sqrt{N_{\mathrm{traj}}}$. Here, we use a smaller number of trajectories $N_{\mathrm{traj}}=100$ than default, to amplify the error and make it more visible in this plot.