Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Log inSign up
Challenges
Datasets
Workspaces
Papers
Discussions
Leaderboard
Blog
Job Board
Q3AS

© 2026 Aqora Quantum S.A.S.

TermsPrivacyLegal Notice

arXiv:2608.30860

Well-conditioned iterative methods for large open quantum systems

Gaspard Beugnot, Paul Gregory, Rémi Robin, Antoine Tilloy

quant-phmath.NASubmitted Aug 31, 2026• Announced Sep 1, 2026• 0 citations
arXivHTMLPDF

Abstract

Markovian open quantum systems are well modeled by the Lindblad Master Equation (ME) $\frac{\mathrm{d}}{\mathrm{d} t} ρ_t = \mathcal{L} ρ_t$, where $\mathcal{L}$ is a linear (super-)operator and $ρ_t$ is the system state, a positive matrix. When designing or characterizing a quantum system, one is usually interested in the steady state $ρ_\infty$ (such that $\mathcal{L} ρ_\infty = 0$), the first few excited states, and trajectories $t\mapsto ρ_t$. In finite dimension, $ρ_t$ is an $n\times n$ matrix, $\mathcal{L}$ thus typically costs $n^4$ to store explicitly as a dense matrix, and $O(n^6)$ to diagonalize or invert exactly, making standard linear algebraic techniques expensive for large systems. However, $\mathcal{L}$ usually costs only $O(n^3)$ to apply. This makes iterative methods appealing, but they do not work without a good preconditioner. In this article, our main observation is that a part of the Lindblad equation, corresponding to the so-called no-jump evolution $\mathcal{S}$, can be inverted efficiently. Using this inverse map, we introduce an auxiliary completely positive trace-preserving (CPTP) map $Φ$ whose fixed point is directly related to $ρ_\infty$, all the other eigenvalues having smaller magnitude. The map $Φ$ is thus well suited to iterative methods, and $ρ_\infty$ can be found in a few Arnoldi iterations. Using the same inverse map $\mathcal{S}^{-1}$ as preconditioner, we compute the low-lying spectrum efficiently via shift-invert Arnoldi, and, as a proof of concept, build an implicit time integrator that is competitive on stiff systems in the low-precision regime. For the steady-state and low excited states problems, our methods scale like $O(n^3)$ per iteration and offer state-of-the-art performance on CPU and GPU.

Related papers

Completely Positive Trace Preserving Methods for the Lindblad Equation

Michael Riesch, Alek Pikl, Christian Jirauschek

The Lindblad master equation is a valuable tool in quantum mechanics, which describes the dynamics of open systems. In the scope of our research, it is combined with the one-dimensional Maxwell's equations to form the generalized Maxwell-Bloch equations. Since analytical solutions are not available in the general case, numerical methods have to be employed to solve the Lindblad equation. In this work, we focus on methods that are completely positive trace preserving (CPTP), i.e., that guarantee to preserve the properties of the density matrix. We review existing approaches and compare the most promising candidates in terms of computational performance.

International Conference on Numerical Simulation of Optoelectronic Devices (NUSOD)•Sep 1, 2020•2 citations

Simulating Non-Markovian Dynamics in Open Quantum Systems

Meng Xu, Vasilii Vadimov, J. T. Stockburger +1 more

Recent advances in quantum technologies and related experiments have created a need for highly accurate, versatile, and computationally efficient simulation techniques for the dynamics of open quantum systems. Long-lived correlation effects (non-Markovianity), system-environment hybridization, and the necessity for accuracy beyond the Born-Markov approximation form particular challenges. Approaches to meet these challenges have been introduced, originating from different fields, such as hierarchical equations of motion, Lindblad-pseudomode formulas, chain-mapping approaches, quantum Brownian motion master equations, stochastic unravelings, and refined quantum master equations. This diversity, while indicative of the field's relevance, has inadvertently led to a fragmentation that hinders cohesive advances and their effective cross-community application to current problems for complex systems. How are different approaches related to each other? What are their strengths and limitations? Here we give a systematic overview and concise discussion addressing these questions. We make use of a unified framework which very conveniently allows to link different schemes and, this way, may also catalyze further progress. In line with the state of the art, this framework is formulated not in a fully reduced space of the system but in an extended state space which in a minimal fashion includes effective reservoir modes. This in turn offers a comprehensive understanding of existing methods, elucidating their physical interpretations, interconnections, and applicability.

quant-pharXiv:2601.02160•Reviews of Modern Physics•Jan 6, 2026•11 citations

Arnoldi-Lindblad time evolution: Faster-than-the-clock algorithm for the spectrum of time-independent and Floquet open quantum systems

Fabrizio Minganti, Dolf Huybrechts

The characterization of open quantum systems is a central and recurring problem for the development of quantum technologies. For time-independent systems, an (often unique) steady state describes the average physics once all the transient processes have faded out, but interesting quantum properties can emerge at intermediate timescales. Given a Lindblad master equation, these properties are encoded in the spectrum of the Liouvillian whose diagonalization, however, is a challenge even for small-size quantum systems. Here, we propose a new method to efficiently provide the Liouvillian spectral decomposition. We call this method an Arnoldi-Lindblad time evolution, because it exploits the algebraic properties of the Liouvillian superoperator to efficiently construct a basis for the Arnoldi iteration problem. The advantage of our method is double: (i) It provides a faster-than-the-clock method to efficiently obtain the steady state, meaning that it produces the steady state through time evolution shorter than needed for the system to reach stationarity. (ii) It retrieves the low-lying spectral properties of the Liouvillian with a minimal overhead, allowing to determine both which quantum properties emerge and for how long they can be observed in a system. This method is general and model-independent , and lends itself to the study of large systems where the determination of the Liouvillian spectrum can be numerically demanding but the time evolution of the density matrix is still doable. Our results can be extended to time evolution with a time-dependent Liouvillian. In particular, our method works for Floquet (i.e., periodically driven) systems, where it allows not only to construct the Floquet map for the slow-decaying processes, but also to retrieve the stroboscopic steady state and the eigenspectrum of the Floquet map. Although the method can be applied to any Lindbladian evolution (spin, fermions, bosons, …), for the sake of simplicity we demonstrate the efficiency of our method on several examples of coupled bosonic resonators (as a particular example). Our method outperforms other diagonalization techniques and retrieves the Liouvillian low-lying spectrum even for system sizes for which it would be impossible to perform exact diagonalization.

arXiv:2109.01648•Quantum•Feb 10, 2022•15 citations

Succinct Description and Efficient Simulation of Non-Markovian Open Quantum Systems

Xiantao Li, Chunhao Wang

arXiv:2111.03240•Communications in Mathematical Physics•Jan 27, 2023•9 citations

Arnoldi-Lindblad time evolution: Faster-than-the-clock algorithm for the spectrum of (Floquet) open quantum systems

Fabrizio Minganti, Dolf Huybrechts

The characterization of open quantum systems is a central and recurring problem for the development of quantum technologies. For time-independent systems, an (often unique) steady state describes the average physics once all the transient processes have faded out, but interesting quantum properties can emerge at intermediate timescales. Given a Lindblad master equation, these properties are encoded in the spectrum of the Liouvillian whose diagonalization, however, is a challenge even for small-size quantum systems. Here, we propose a new method to efficiently provide the Liouvillian spectral decomposition. We call this method an Arnoldi-Lindblad time evolution, because it exploits the algebraic properties of the Liouvillian superoperator to efficiently construct a basis for the Arnoldi iteration problem. The advantage of our method is double: (i) It provides a faster-than-the-clock method to efficiently obtain the steady state, meaning that it produces the steady state through time evolution shorter than needed for the system to reach stationarity. (ii) It retrieves the low-lying spectral properties of the Liouvillian with a minimal overhead, allowing to determine both which quantum properties emerge and for how long they can be observed in a system. This method is general and model-independent, and lends itself to the study of large systems where the determination of the Liouvillian spectrum can be numerically demanding. Our results can be extended to time evolution with a time-dependent Liouvillian. In particular, our method works for Floquet (i.e., periodically driven) systems, where it allows not only to construct the Floquet map for the slow-decaying processes, but also to retrieve the stroboscopic steady state and the eigenspectrum of the Floquet map. We demonstrate the efficiency of our method on several examples of coupled bosonic systems.

•Sep 3, 2021•0 citations

A Simple GPU-Accelerated Solver for the Schrödinger Operator with Applications to Ground States and Hamiltonian Simulation

Xinyu Liu, Xiangxiong Zhang

We extend the tensor-product direct solver from the Laplacian to the Schrödinger operator $-Δ+ V$. When the potential $V_1$ is separable, the operator $-Δ+ V_1$ is inverted or exponentiated at cost $O(N^{1+1/d})$ in $d$ dimensions via per-axis eigendecomposition. On a single NVIDIA A100 GPU, this costs less than one second for $10^9$ degrees of freedom in 3D. For non-separable potentials $V = V_1 + V_2$, the same solver provides a preconditioner $(-Δ+ V_1)^{-1}$ for the preconditioned conjugate gradient (PCG) method and a propagator for operator-splitting time integrators. For bounded $V_2$, we prove that the preconditioned operator has a bounded condition number and a clustered spectrum with at most finitely many outlier eigenvalues, independently of the mesh size, and also independently of the domain size when $V_1$ is a confining potential. This explains the mesh- and domain-independent PCG iteration counts observed in practice. We apply this method to ground state computation via inverse iteration for linear problems and via the $a_u$ gradient flow for Gross--Pitaevskii energy in 3D, and also Hamiltonian simulation via the approximated qHOP and Magnus-2 splitting methods from 3D to 9D on a single NVIDIA GH200 GPU.

math.NACross-listedarXiv:2605.20491•May 21, 2026
Markovian open quantum systems are well modeled by the Lindblad Master Equation (ME) ddtρt=Lρt\frac{\mathrm{d}}{\mathrm{d} t} ρ_t = \mathcal{L} ρ_tdtd​ρt​=Lρt​, where L\mathcal{L}L is a linear (super-)operator and ρtρ_tρt​ is the system state, a positive matrix. When designing or characterizing a quantum system, one is usually interested in the steady state ρ∞ρ_\inftyρ∞​ (such that Lρ∞=0\mathcal{L} ρ_\infty = 0Lρ∞​=0), the first few excited states, and trajectories t↦ρtt\mapsto ρ_tt↦ρt​. In finite dimension, ρtρ_tρt​ is an n×nn\times nn×n matrix, L\mathcal{L}L thus typically costs n4n^4n4 to store explicitly as a dense matrix, and O(n6)O(n^6)O(n6) to diagonalize or invert exactly, making standard linear algebraic techniques expensive for large systems. However, L\mathcal{L}L usually costs only O(n3)O(n^3)O(n3) to apply. This makes iterative methods appealing, but they do not work without a good preconditioner. In this article, our main observation is that a part of the Lindblad equation, corresponding to the so-called no-jump evolution S\mathcal{S}S, can be inverted efficiently. Using this inverse map, we introduce an auxiliary completely positive trace-preserving (CPTP) map ΦΦΦ whose fixed point is directly related to ρ∞ρ_\inftyρ∞​, all the other eigenvalues having smaller magnitude. The map ΦΦΦ is thus well suited to iterative methods, and ρ∞ρ_\inftyρ∞​ can be found in a few Arnoldi iterations. Using the same inverse map S−1\mathcal{S}^{-1}S−1 as preconditioner, we compute the low-lying spectrum efficiently via shift-invert Arnoldi, and, as a proof of concept, build an implicit time integrator that is competitive on stiff systems in the low-precision regime. For the steady-state and low excited states problems, our methods scale like O(n3)O(n^3)O(n3) per iteration and offer state-of-the-art performance on CPU and GPU.
The Lindblad master equation is a valuable tool in quantum mechanics, which describes the dynamics of open systems. In the scope of our research, it is combined with the one-dimensional Maxwell's equations to form the generalized Maxwell-Bloch equations. Since analytical solutions are not available in the general case, numerical methods have to be employed to solve the Lindblad equation. In this work, we focus on methods that are completely positive trace preserving (CPTP), i.e., that guarantee to preserve the properties of the density matrix. We review existing approaches and compare the most promising candidates in terms of computational performance.
Recent advances in quantum technologies and related experiments have created a need for highly accurate, versatile, and computationally efficient simulation techniques for the dynamics of open quantum systems. Long-lived correlation effects (non-Markovianity), system-environment hybridization, and the necessity for accuracy beyond the Born-Markov approximation form particular challenges. Approaches to meet these challenges have been introduced, originating from different fields, such as hierarchical equations of motion, Lindblad-pseudomode formulas, chain-mapping approaches, quantum Brownian motion master equations, stochastic unravelings, and refined quantum master equations. This diversity, while indicative of the field's relevance, has inadvertently led to a fragmentation that hinders cohesive advances and their effective cross-community application to current problems for complex systems. How are different approaches related to each other? What are their strengths and limitations? Here we give a systematic overview and concise discussion addressing these questions. We make use of a unified framework which very conveniently allows to link different schemes and, this way, may also catalyze further progress. In line with the state of the art, this framework is formulated not in a fully reduced space of the system but in an extended state space which in a minimal fashion includes effective reservoir modes. This in turn offers a comprehensive understanding of existing methods, elucidating their physical interpretations, interconnections, and applicability.
The characterization of open quantum systems is a central and recurring problem for the development of quantum technologies. For time-independent systems, an (often unique) steady state describes the average physics once all the transient processes have faded out, but interesting quantum properties can emerge at intermediate timescales. Given a Lindblad master equation, these properties are encoded in the spectrum of the Liouvillian whose diagonalization, however, is a challenge even for small-size quantum systems. Here, we propose a new method to efficiently provide the Liouvillian spectral decomposition. We call this method an Arnoldi-Lindblad time evolution, because it exploits the algebraic properties of the Liouvillian superoperator to efficiently construct a basis for the Arnoldi iteration problem. The advantage of our method is double: (i) It provides a faster-than-the-clock method to efficiently obtain the steady state, meaning that it produces the steady state through time evolution shorter than needed for the system to reach stationarity. (ii) It retrieves the low-lying spectral properties of the Liouvillian with a minimal overhead, allowing to determine both which quantum properties emerge and for how long they can be observed in a system. This method is general and model-independent , and lends itself to the study of large systems where the determination of the Liouvillian spectrum can be numerically demanding but the time evolution of the density matrix is still doable. Our results can be extended to time evolution with a time-dependent Liouvillian. In particular, our method works for Floquet (i.e., periodically driven) systems, where it allows not only to construct the Floquet map for the slow-decaying processes, but also to retrieve the stroboscopic steady state and the eigenspectrum of the Floquet map. Although the method can be applied to any Lindbladian evolution (spin, fermions, bosons, …), for the sake of simplicity we demonstrate the efficiency of our method on several examples of coupled bosonic resonators (as a particular example). Our method outperforms other diagonalization techniques and retrieves the Liouvillian low-lying spectrum even for system sizes for which it would be impossible to perform exact diagonalization.
The characterization of open quantum systems is a central and recurring problem for the development of quantum technologies. For time-independent systems, an (often unique) steady state describes the average physics once all the transient processes have faded out, but interesting quantum properties can emerge at intermediate timescales. Given a Lindblad master equation, these properties are encoded in the spectrum of the Liouvillian whose diagonalization, however, is a challenge even for small-size quantum systems. Here, we propose a new method to efficiently provide the Liouvillian spectral decomposition. We call this method an Arnoldi-Lindblad time evolution, because it exploits the algebraic properties of the Liouvillian superoperator to efficiently construct a basis for the Arnoldi iteration problem. The advantage of our method is double: (i) It provides a faster-than-the-clock method to efficiently obtain the steady state, meaning that it produces the steady state through time evolution shorter than needed for the system to reach stationarity. (ii) It retrieves the low-lying spectral properties of the Liouvillian with a minimal overhead, allowing to determine both which quantum properties emerge and for how long they can be observed in a system. This method is general and model-independent, and lends itself to the study of large systems where the determination of the Liouvillian spectrum can be numerically demanding. Our results can be extended to time evolution with a time-dependent Liouvillian. In particular, our method works for Floquet (i.e., periodically driven) systems, where it allows not only to construct the Floquet map for the slow-decaying processes, but also to retrieve the stroboscopic steady state and the eigenspectrum of the Floquet map. We demonstrate the efficiency of our method on several examples of coupled bosonic systems.
We extend the tensor-product direct solver from the Laplacian to the Schrödinger operator −Δ+V-Δ+ V−Δ+V. When the potential V1V_1V1​ is separable, the operator −Δ+V1-Δ+ V_1−Δ+V1​ is inverted or exponentiated at cost O(N1+1/d)O(N^{1+1/d})O(N1+1/d) in ddd dimensions via per-axis eigendecomposition. On a single NVIDIA A100 GPU, this costs less than one second for 10910^9109 degrees of freedom in 3D. For non-separable potentials V=V1+V2V = V_1 + V_2V=V1​+V2​, the same solver provides a preconditioner (−Δ+V1)−1(-Δ+ V_1)^{-1}(−Δ+V1​)−1 for the preconditioned conjugate gradient (PCG) method and a propagator for operator-splitting time integrators. For bounded V2V_2V2​, we prove that the preconditioned operator has a bounded condition number and a clustered spectrum with at most finitely many outlier eigenvalues, independently of the mesh size, and also independently of the domain size when V1V_1V1​ is a confining potential. This explains the mesh- and domain-independent PCG iteration counts observed in practice. We apply this method to ground state computation via inverse iteration for linear problems and via the aua_uau​ gradient flow for Gross--Pitaevskii energy in 3D, and also Hamiltonian simulation via the approximated qHOP and Magnus-2 splitting methods from 3D to 9D on a single NVIDIA GH200 GPU.