Well-conditioned iterative methods for large open quantum systems

AI summary

Overview

Research area: Quantum physics — specifically open quantum systems and numerical methods for the Lindblad master equation, sitting at the intersection of quantum information, quantum control, and numerical linear algebra.

Technical level: Advanced. The paper assumes familiarity with superoperators, completely positive trace-preserving (CPTP) maps, Krylov subspace methods, and Lyapunov equations.

Scope (one sentence): The paper introduces iterative, preconditioned numerical methods that exploit the cheap invertibility of the no-jump part of the Lindbladian to compute steady states, low-lying spectra, and time evolution of large Markovian open quantum systems at O(n^3) cost per iteration.

What This Paper Is About

For a Markovian open quantum system on a Hilbert space of dimension n, the state is an n × n density matrix and the generator (the Lindbladian ℒ) acts on an n^2-dimensional space. Storing ℒ densely costs n^4, and diagonalizing or inverting it exactly costs O(n^6), which becomes prohibitive beyond n ≳ 100. Yet applying ℒ to a state costs only O(n^3) when the number of jump operators is O(1), which is the typical physical situation. The problem is that ℒ is badly ill-conditioned, so naive iterative methods perform poorly. The goal of this paper is to build well-conditioned iterative methods by exploiting the fact that the no-jump part of the Lindbladian can be inverted efficiently.

Key Contributions

1. A CPTP reformulation of the steady-state problem. The authors show that the steady state can be obtained as the fixed point of an auxiliary completely positive trace-preserving map Φ = -𝒦𝒮⁻¹, where 𝒮 is the no-jump generator and 𝒦 the jump generator. Because Φ is CPTP, its spectrum lies in the closed unit disk and the eigenvalue of interest (1) is of maximal modulus — exactly the regime where Arnoldi/Krylov methods are efficient, in contrast with ℒ whose target eigenvalue is the smallest in magnitude.

2. A structural theorem (Theorem 1) separating two cases. Either the spectrum of G = -iH - ½Σ L_j†L_j touches the imaginary axis, in which case the steady state is pure and is found directly from an eigenstate of H annihilated by all jump operators; or the spectrum lies strictly in the left half-plane, in which case 𝒮 is invertible, Φ is a well-defined CPTP map with a unique invariant state ξ∞, and ρ∞ = -c 𝒮⁻¹ ξ∞ for a positive constant c.

3. A preconditioner built from the no-jump resolvent. The paper introduces ℛ_λ^𝒮 = (λ - 𝒮)⁻¹ for λ > 0, computable at O(n^3) via a continuous Lyapunov equation, and proves (Theorem 2) that 𝒦ℛ_λ^𝒮 is a strict contraction in the trace-induced norm, with norm equal to 1 - λδ_λ < 1. This justifies using ℛ_λ^𝒮 as a right preconditioner for the shifted Lindbladian λ - ℒ.

4. Three applications of that preconditioner: steady-state computation as a linear system, low-lying spectrum computation via shift-invert Arnoldi, and a proof-of-concept implicit time integrator that is competitive on stiff systems in the low-precision regime. All scale as O(n^3) per iteration and are benchmarked on CPU and GPU.

Main Findings

• Steady state found in a few Arnoldi iterations. Because 1 is a simple eigenvalue of maximal modulus of Φ, Arnoldi applied to Φ converges rapidly. The authors retain the Ritz vector whose Ritz value is closest to 1, rather than the dominant one, to guard against peripheral modes.

• Peripheral eigenvalues can exist, unlike for the original Lindbladian. Remark 3 gives a three-level example with L_1 = |0⟩⟨1|, L_2 = |1⟩⟨2|, L_3 = |2⟩⟨0| and no Hamiltonian, whose unique steady state is the maximally mixed state Id/3, yet Φ has three peripheral eigenvalues 1, e^{2πi/3}, e^{-2πi/3}. Power iteration would fail here; Arnoldi is used instead.

• Kernel correspondence holds even without uniqueness. Remark 2 notes that ker ℒ = -𝒮⁻¹ ker(Φ - Id) holds whenever 𝒮 is invertible, and that dim ker(Φ - Id) = dim ker ℒ; uniqueness of the steady state is only used to conclude both kernels are one-dimensional.

• Contraction bound is explicit. Theorem 2 gives ‖𝒦ℛ_λ^𝒮‖ = 1 - λδ_λ in the trace-induced norm, where δ_λ is the minimum of Tr ℛ_λ^𝒮 ρ over density matrices, hence strictly below one for λ > 0.

• Connection to an existing CPTP integrator. Remark 4 shows that a recently proposed first-order CPTP time discretization can be rewritten as e^{δtG} ρ_n e^{δtG†} + Φ(ρ_n - e^{δtG} ρ_n e^{δtG†}), and that in the limit δt → +∞ the scheme reduces to ρ_{n+1} = Φ(ρ_n), whose fixed point is related to the true steady state by Theorem 1.

• Benchmark results. On a dense random Lindbladian (random Hermitian Hamiltonian plus three dense jump operators with i.i.d. complex-normal entries) and on a sparse memory-buffer cat-qubit Lindbladian (four sparse jump operators, small spectral gap), the Arnoldi-on-Φ method and the preconditioned GMRES solver outperform the competitors on CPU and gain a further ×50 speedup on an Nvidia H100 GPU. Competitors were QuTiP SVD, QuTiP ILU-GMRES, MUMPS, Krylov.jl GMRES, and a time-evolution baseline; the dense SVD solver shows the expected O(n^6) scaling, and the unpreconditioned and ILU-preconditioned GMRES baselines scale poorly and hit memory limits at larger dimensions. MUMPS performs well on CPU but is unavailable on GPU and cannot handle the dense operators of the random benchmark. The ODE integration baseline is much slower in absolute runtime on the cat qubit due to strong metastability from exponential suppression of bit-flips.

• Spectral illustration at small size. Figure 1 plots the spectra of ℒ and Φ for a random dense Lindbladian at n = 6, showing the wide eigenvalue spread of ℒ versus the unit-disk confinement of Φ.

Methodology in Plain English

The starting point is a split of the Lindblad equation into two pieces: the "no-jump" part, which describes how the state evolves when no quantum jump is detected, and the "jump" part, which describes the effect of the jumps themselves. The key practical fact is that inverting the no-jump part amounts to solving a continuous Lyapunov equation, a classic problem from control theory that can be solved in O(n^3) rather than the O(n^6) one would fear from inverting an operator on an n^2-dimensional space.

The authors use this cheap inverse in two ways. First, they compose it with the jump part to build a new map whose fixed point is the steady state, and whose eigenvalue structure is favorable for Krylov iteration. Second, for problems where the eigenvalue structure of that map is not directly useful — such as finding excited states — they use the same inverse as a preconditioner for the shifted Lindbladian, so that the iterative solver only has to resolve the part of the problem that the Lyapunov solve does not already handle.

Numerically, the methods are tested against standard alternatives on a single node with an Intel Xeon Platinum 8481C CPU and 70 GB of RAM, and on an Nvidia H100 GPU with 80 GB, in double precision. Convergence is declared when every entry of ℒ(ρ) is below 10^-8 in absolute value after enforcing that ρ is Hermitian with trace 1, and runs are cut off at 100 s or when memory is exhausted. The reported quantity is wall-clock time as a function of Hilbert space dimension n.

Why This Matters

Impact on research. Steady states of the GKSL equation determine long-time observables, are used to calibrate experimental parameters, and appear in the evaluation of expectation values of continuous matrix product states (CMPS) in the thermodynamic limit. Faster steady-state and spectral solvers directly enlarge the system sizes that can be simulated, and the paper's observation that only the no-jump part needs to be inverted exactly gives a general recipe that other practitioners can reuse.

Real-world applications:

• Quantum state preparation by dissipation: designing a Lindbladian whose unique stationary state is the desired state, which requires solving for that state repeatedly during design.
• Device calibration: extracting physical parameters from measured steady-state observables.
• Bosonic and cat-qubit hardware: the paper benchmarks a memory-buffer cat-state model representative of dynamics in Alice & Bob quantum processors, where metastability makes long-time integration impractical.
• Tensor network methods: computing CMPS expectation values in the thermodynamic limit, which reduces to a Lindblad steady-state problem.

Industry relevance. Quantum hardware companies need to characterize and design dissipative dynamics for large devices; the demonstrated ×50 GPU speedup and O(n^3) per-iteration scaling mean larger, more realistic device models become tractable on commodity accelerated hardware. The methods are also directly compatible with existing sparse and GPU linear algebra stacks.

Future Directions

• Optimizing the no-jump/jump split. The decomposition ℒ = 𝒮 + 𝒦 is not unique: shifting each jump operator by a multiple of the identity, L_j → L_j + c_j Id, with a compensating Hamiltonian change, yields other admissible splits. The authors state that they stick to the c = 0 split throughout and defer discussion of optimizing over c to the conclusion, so choosing c to improve conditioning is an open avenue.
• Scaling up tensor network calculations. Since efficient steady-state computation would allow working with far larger CMPS, applying these solvers inside CMPS pipelines is a natural next step.
• Production-grade implicit time integrators. The implicit integrator is presented explicitly as a proof of concept that is competitive only on stiff systems in the low-precision regime; extending it to broader regimes and higher accuracy is left open.
• Truncation error control for bosonic systems. For the infinite-dimensional examples (bosonic systems) the Hilbert space is truncated to finite dimension and truncation errors are explicitly not considered, leaving a gap between the linear-algebra benchmarks and physically converged results.

Target Audience

This paper is aimed at researchers and practitioners comfortable with quantum master equations and numerical linear algebra: quantum information and open-quantum-systems theorists, quantum control engineers, and scientific computing specialists working on Krylov methods and preconditioners. It is also relevant to computational physicists simulating bosonic or superconducting-qubit hardware who need steady states or low-lying Lindblad spectra at sizes beyond n ≈ 100. Readers without a background in superoperators, CPTP maps, and GMRES will find the theorems and benchmark discussion demanding.

Authors’ 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.