QuantumToolbox.jl: An efficient Julia framework for simulating open quantum systems

TL;DR

QuantumToolbox.jl delivers a fast, flexible Julia framework for simulating open quantum systems with a QuTiP-like API, leveraging backend-agnostic design, GPU acceleration, and Distributed.jl support. It provides a comprehensive suite of solvers for Schrödinger, Lindblad, Monte Carlo, time-dependent, and stochastic dynamics, along with dynamic Hilbert-space adaptation (DFD/DSF) and experimental automatic differentiation for gradient-based optimization. The paper demonstrates strong performance gains over existing tools across CPU and GPU backends, and showcases advanced capabilities such as steady-state Fourier analysis and quantum trajectory parallelism. The platform aims to accelerate both theoretical investigations and practical quantum technologies, while highlighting current AD limitations and future expansion opportunities.

Abstract

We present QuantumToolbox$.$jl, an open-source Julia package for simulating open quantum systems. Designed with a syntax familiar to users of QuTiP (Quantum Toolbox in Python), it harnesses Julia's high-performance ecosystem to deliver fast and scalable simulations. The package includes a suite of time-evolution solvers supporting distributed computing and GPU acceleration, enabling efficient simulation of large-scale quantum systems. We also show how QuantumToolbox$.$jl can integrate with automatic differentiation tools, making it well-suited for gradient-based optimization tasks such as quantum optimal control. Benchmark comparisons demonstrate substantial performance gains over existing frameworks. With its flexible design and computational efficiency, QuantumToolbox$.$jl serves as a powerful tool for both theoretical studies and practical applications in quantum science.

Figures (8)

• Figure 1: Ecosystem of the QuantumToolbox.jl package.(a) The quantum object structure encapsulates the properties of generic quantum objects through three main fields: , , and . The field specifies the category of the quantum object, the field defines the Hilbert space structure for composite systems, and the field supports arbitrary array types in Julia. (b) Various time evolution solvers are available for studying open quantum system dynamics. (c) All quantum objects support fundamental arithmetic and linear algebra operations. (d)QuantumToolbox.jl also provides a variety of utility functions for analyzing the properties and physical quantities of a given quantum object. (e) Logo of the QuantumToolbox.jl package.
• Figure 2: Examples of time evolution solvers. We observe the average photon number of the cavity $\langle \hat{a}^\dagger \hat{a} \rangle$ with respect to time in the Jaynes-Cummings model ((a), (b) and (c)) and in the optomechanical system ((d)). (a) shows the time evolution following the Schrödinger equation, while (b) and (c) show the time evolution of the Lindblad master equation and the Monte Carlo wave-function, respectively. The dashed light blue curve in (c) represents a single quantum trajectory, where a quantum jump occurs at $t \approx 200$, and the solid blue curve represents the average over 100 quantum trajectories. (d) shows the time evolution of the optomechanical system, where the cavity is driven by a time-dependent drive. The dashed red and solid orange curves represent the value of the cavity population at the time-averaged steady state using $n_\mathrm{max}=1$ and $n_\mathrm{max}=2$ fourier components, respectively.
• Figure 3: Stochastic Schrödinger equation and stochastic master equation.(a) Time evolution under the stochastic Schrödinger equation of the homodyne current $J_x$ (light blue curve) and the expectation value $\langle \hat{X} \rangle$ (dark blue curve) for the Jaynes-Cummings model. (b) Time evolution under the stochastic master equation of the same quantities as in (a). Both the homodyne current and the expectation value are averaged over 500 trajectories.
• Figure 4: Dynamical shifted Fock algorithm.(a) Time evolution of the average photon number $\langle \hat{a}_1^\dagger \hat{a}_1 \rangle$ for two nonlinearly coupled harmonic oscillators with Kerr nonlinearity. In contrast to the second-order quantum cumulant expansion (dash-dotted green curve), the DSF algorithm (dashed orange curve) is able to capture the system’s dynamics while still maintaining a low-dimensional Hilbert space in the shifted Fock basis. (b-c) Wigner function of the first cavity field at the initial time ((b)), and at the final time ((c)). The solid blue circle represents the Hilbert space required to simulate the model using , while the dashed orange circle represents the Hilbert space used by the DSF algorithm. As can be seen, the center of the space is shifted, depending on the coherence of the state. (d) Average photon number of the second mode, which is simulated in the original Fock basis during the DSF algorithm solving process. However, the quantum cumulants expansion is not able to reproduce this behavior. (e-f) Wigner function of the second mode at the initial and final time, respectively. The Wigner function is computed using the function.
• Figure 5: Dynamical shifted Fock algorithm.(a) Time evolution of the average photon number $\langle \hat{a}^\dagger \hat{a} \rangle$ for the driven Jaynes--Cummings model with Kerr nonlinearity. Even though the cavity state is not strictly coherent, its localization in phase space allows the DSF algorithm to reproduce the dynamics accurately, while the quantum cumulant expansion fails to capture the correct behavior. (b--d) Wigner functions of the cavity state at different times, illustrating its localized, non-coherent structure and the efficiency of the DSF representation with $N_\mathrm{dsf}=13$.