Efficient quantum thermal simulation

Abstract

Quantum computers promise to tackle quantum simulation problems that are classically intractable¹. Although a lot of quantum algorithms^2,3,4 have been developed for simulating quantum dynamics, a general-purpose method for simulating low-temperature quantum phenomena remains unknown. In classical settings, the analogous task of sampling from thermal distributions has been largely addressed by Markov Chain Monte Carlo (MCMC) methods^5,6. Here we propose an efficient quantum algorithm for thermal simulation that—akin to MCMC methods—exhibits detailed balance, respects locality and serves as a toy model for thermalization in open quantum systems. The enduring impact of MCMC methods suggests that our new construction may play an equally important part in quantum computing and applications in the physical sciences and beyond.

Main

The equations governing quantum many-body physics are too complex to be solved analytically. Therefore, an arsenal of computational methods—including density functional theory⁷, tensor network representations^8,9, and quantum Monte Carlo techniques¹⁰—has been developed for understanding low-temperature phenomena, ranging from phase transitions and reaction times to ground and thermal state equilibrium properties. However, these classical methods often break down in the presence of strong quantum correlations, for example, exhibited by heavy nuclei, high-temperature superconductors and catalysts.

Future quantum computers promise to solve quantum many-body problems that are classically intractable¹. At present, a variety of increasingly sophisticated quantum algorithms have been developed for simulating quantum dynamics^2,11,12. By contrast, for equilibrium properties, the search for particular physical instances in which quantum computers can substantially outperform classical methods remains an area of debate. The efficiency of most existing quantum algorithmic proposals for low-temperature physics^13,14,15,16 relies on additional assumptions, which have been challenged in numerical studies¹⁷.

Classically, the many-body thermal simulation problem has been thoroughly addressed, largely by Markov Chain Monte Carlo methods (MCMC). Across a wide range of systems, Metropolis sampling⁵ and Glauber dynamics⁶ have been successful because of two key features. First, they satisfy detailed balance, which simultaneously lets us (1) prescribe a desired stationary distribution and (2) analyse their convergence (that is, mixing times and spectral gaps). Second, MCMC methods are often local algorithms, that is, the update rule depends only on very few particles at a time; locality gives not only efficient algorithms but also physical connections to how nature thermalizes in a heat bath.

However, quantum computing still lacks a definitive counterpart to MCMC algorithms, as existing proposals^18,19,20,21 are nonlocal, break detailed balance, lack conceptual simplicity or require additional hard-to-verify assumptions. Moreover, the algorithm and its analysis in the promising early approach¹⁸ turns out to be incomplete (see Appendix H of ref. ²²). As a consequence, progress in quantum thermal simulation remained limited for decades, raising the question whether exact quantum detailed balance should be compatible with the time–energy uncertainty principle—the primary source of difficulties rooted in quantum many-body Hamiltonians with noncommuting terms.

Here, we settle the open problem by giving an efficient quantum algorithm for quantum thermal simulation that satisfies quantum detailed balance exactly while also inheriting the locality of the Hamiltonian, up to a range depending on the temperature. Moreover, our algorithm resembles effective descriptions of open quantum system dynamics when coupling a system to a heat bath. Qualitatively, we expect our algorithm to serve as a minimal, self-contained toy model of thermalization: converging rapidly whenever a system thermalizes in nature and capturing what may physically happen in the case of slow mixing. We speculate that our algorithm might be among the first useful quantum simulation algorithms for low-temperature properties in the forthcoming fault-tolerant era of quantum computing. Furthermore, given the key role played by MCMC algorithms in problems across bioinformatics, social sciences, finance, Bayesian inference and machine learning, we expect an efficient and reliable quantum MCMC algorithm to be a cornerstone of quantum algorithmic development and quantum applications in physical sciences and beyond.

Quantum MCMC by master equations for thermalization

Our algorithmic solution to thermal simulation hinges on insights into how physical systems thermalize when coupled to a heat bath. In essence, we introduce a conceptually clear, physically motivated toy model of the thermalization process of nature that is also algorithmically efficient. Although early proposals for simulating thermalization directly considered simulating the global system–bath Hamiltonian evolution²³, more tractable approaches use a master equation, or Lindbladian, to capture the effect of a large bath on the small system using a continuous-time Markovian process^24,25. However, previous approaches inherited issues of the prototypical Davies generator^26,27, which worked well in quantum optics but has unphysical features in noncommuting many-body systems because of the exponentially small level spacing. Our nature-inspired algorithm takes precisely the form of a Lindbladian, inherits the locality from the physical model, exactly satisfies detailed balance and resembles the interactions we expect from weak coupling to a Markovian thermal bath.

Continuous-time Markov chains and detailed balance

We begin by reviewing classical MCMC methods, in particular, continuous-time Markov chains defined using the rate matrix L describing infinitesimal transitions on probability distributions p(t)

$$\frac{{\rm{d}}{\rm{p}}(t)}{{\rm{d}}t}=L{\rm{p}}(t),$$

(1)

where L satisfies that ${L}_{{s}^{{\prime} }s}\ge 0$ for each pair of distinct configurations s ≠ s′ and ${L}_{ss}=-\sum _{{s}^{{\prime} }\ne s}{L}_{{s}^{{\prime} }s}$ for each configuration s. These structural properties ensure that the time evolution can be described by the stochastic matrix p(t) = e^Ltp(0).

The detailed balance condition plays an important part in MCMC methods, requiring that the probability mass transfer between configuration pairs s, s′ is symmetric with respect to a target distribution π

$${L}_{{s}^{{\prime} }s}{{\rm{\pi }}}_{s}={L}_{s{s}^{{\prime} }}{{\rm{\pi }}}_{{s}^{{\prime} }},\,\text{which then implies stationarity}\,L{\rm{\pi }}=0.$$

(2)

If the target distribution π is the unique stationary state, we say the dynamics e^Lt is ergodic, and the algorithmic cost to sample from the stationary state is the simulation cost (per unit time t) multiplied by the convergence time (the mixing time) of the process.

Given an energy function H(s) of the configurations s, an inverse temperature β ≥ 0, and any real nonnegative symmetric transition rate matrix ${A}_{{s}^{{\prime} }s}$, Glauber dynamics⁶ is a specific continuous-time Markov chain satisfying detailed balance such that ${L}_{{s}^{{\prime} }s}={\gamma }_{{s}^{{\prime} }s}{A}_{{s}^{{\prime} }s}\,-$ ${\delta }_{{s}^{{\prime} }s}\sum _{{s}^{{\prime\prime} }}{\gamma }_{{s}^{{\prime\prime} }s}{A}_{{s}^{{\prime\prime} }s}$, where ${\gamma }_{{s}^{{\prime} }s}:= \frac{1}{1+{{\rm{e}}}^{\beta (H({s}^{{\prime} })-H(s))}}$. The process e^Lt can be efficiently simulated as long as the energy differences $H({s}^{{\prime} })-H(s)$ can be computed efficiently for any pair of configurations s′, s. For example, in classical Ising models, the Hamiltonian is a sum over local interactions, and the energy difference H(s′) − H(s) can be computed locally (for transitions that only flip a few neighbouring spins at a time), leading to particularly simple update rules and a deep analytic understanding of mixing times. The robustness and simplicity of Glauber dynamics have led to a comprehensive understanding of phase transitions and thermal properties of classical magnets²⁸.

Lindbladians and quantum detailed balance

The natural quantum counterpart, which arises in the study of open quantum systems, is a master equation, a continuous-time quantum Markov chain, or Lindbladian described by a set of Lindblad operators {L_j} and a coherent term C:

$$\begin{array}{l}\frac{{\rm{d}}\rho (t)}{{\rm{d}}t}={\mathcal{L}}[\,\rho (t)]\,:= \,-\mathop{\underbrace{i[C,\rho (t)]}}\limits_{ \mbox{``} {\rm{coherent}}\mbox{''}}+\sum _{j}\mathop{\underbrace{{L}_{j}\rho (t){L}_{j}^{\dagger }}}\limits_{ \mbox{``} {\rm{transition}}\mbox{''}}\\ \,\,\,\,\,\,\,\,-\mathop{\underbrace{\frac{1}{2}({L}_{j}^{\dagger }{L}_{j}\rho (t)+\rho (t){L}_{j}^{\dagger }{L}_{j})}}\limits_{ \mbox{``} {\rm{decay}}\mbox{''}}.\end{array}$$

(3)

The above generic form is necessary and sufficient for the evolution ${{\rm{e}}}^{{\mathcal{L}}t}$ at every t ≥ 0 to be a valid quantum Markov chain; that is, a completely positive and trace-preserving map. Unlike in the classical case, quantum mechanics allows an additional Hermitian coherent term C = C^† corresponding to the closed-system Schrödinger equation.

Consider an arbitrary noncommuting Hamiltonian with eigen decomposition $H=\sum _{i}{E}_{i}| {\psi }_{i}\rangle \langle {\psi }_{i}|$. The quantum detailed balance condition with respect to the Gibbs state ρ_β ∝ e^−βH requires

$$\langle {\psi }_{2}^{{\prime} }|{\mathcal{L}}[|{\psi }_{2}\rangle \langle {\psi }_{1}|]|{\psi }_{1}^{{\prime} }\rangle =\exp \left(-\beta \frac{{E}_{1}^{{\prime} }-{E}_{1}+{E}_{2}^{{\prime} }-{E}_{2}}{2}\right){(\langle {\psi }_{2}|{\mathcal{L}}[|{\psi }_{2}^{{\prime} }\rangle \langle {\psi }_{1}^{{\prime} }|]|{\psi }_{1}\rangle )}^{\ast }\,\text{for each}\,|{\psi }_{1}\rangle ,|{\psi }_{1}^{{\prime} }\rangle ,|{\psi }_{2}\rangle ,|{\psi }_{2}^{{\prime} }\rangle ,$$

(4)

which ensures exact stationarity ${\mathcal{L}}[{\rho }_{\beta }]=0$ such that ${{\rm{e}}}^{{\mathcal{L}}t}[{\rho }_{\beta }]={\rho }_{\beta }$ for all t ≥ 0 (refs. ^18,29). The above reduces to classical detailed balance (equation (2)) when the states are diagonal in the energy basis (that is, $| {\psi }_{1}^{{\prime} }\rangle =| {\psi }_{2}^{{\prime} }\rangle$ and |ψ₁⟩ = |ψ₂⟩). In general, the input quantum state can be a mixture of superpositions of different energy eigenstates, and quantum detailed balance (equation (4)) (specifically, the Kubo–Martin–Schwinger (KMS) detailed balance) requires the transition amplitudes between matrix elements to be related by a Boltzmann factor corresponding to the average of the energy differences ${E}_{1}^{{\prime} }-{E}_{1}$ and ${E}_{2}^{{\prime} }-{E}_{2}$ (Fig. 1).

Fig. 1: Illustration of detailed balance.

a, Classical: for classical Gibbs distributions, the detailed balance condition is a pairwise relation between heating (red, ν > 0) and cooling (blue, ν < 0) transition rates, depending on the energy difference ν between states. b, Quantum: for quantum Gibbs states, the detailed balance condition refers to matrix elements of the density operator in the energy eigenbasis. Therefore, the relation depends on the respective energy differences ν₁ and ν₂ of the corresponding ‘kets’ and ‘bras’. In classical Markov chains, detailed balance can be achieved by rejecting a jump with a probability proportional to the energy difference. However, two key obstacles arise in implementing quantum detailed balance: (i) measuring energy difference precisely is algorithmically expensive and may disturb the state and (ii) undoing a jump on a quantum state is intricate as we do not have a priori knowledge of the quantum state beyond its energies, and we cannot keep a copy of quantum states before the jump because of quantum no-cloning theroems. Our construction circumvents these obstacles by smoothly manipulating the transition amplitudes using operator Fourier transforms with energy uncertainty of about 1/β and introducing a Lindbladian with a corrective coherent term.

Davies generators, the prototype of quantum detailed balance

In 1974, Davies²⁶ rigorously derived a Lindbladian from first principles by considering a physical system interacting with a large thermal bath in the weak-coupling and infinite-time limit. Suppose that A is a Hermitian jump operator through which the system is perturbed by the bath (for example, a local Pauli operator); then, the corresponding Lindblad operators L_ν in the Davies generator are given by a certain Fourier transform of the Heisenberg dynamics

$$\begin{array}{c}{L}_{\nu }\,=\,\sqrt{\gamma (\nu )}{A}_{\nu },\,\,{\rm{where}}\,\,\,{A}_{\nu }\,:= \,\sum _{i,j:{E}_{i}-{E}_{j}=\nu }\langle {\psi }_{i}|A|{\psi }_{j}\rangle |{\psi }_{i}\rangle \langle {\psi }_{j}|\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,=\,{{\rm{l}}{\rm{i}}{\rm{m}}}_{T\to {\rm{\infty }}}\frac{1}{2T}{\int }_{-T}^{T}{{\rm{e}}}^{{\rm{i}}Ht}A{{\rm{e}}}^{-{\rm{i}}Ht}{{\rm{e}}}^{-{\rm{i}}\nu t}{\rm{d}}t.\end{array}$$

(5)

These Lindblad operators L_ν are labelled by the Bohr frequencies ν = E_i − E_j ∈ B(H), that is, energy differences between any two eigenstates of the Hamiltonian. Colloquially, A_ν contains ‘the transition amplitudes in A that change the energy by ν’, which is precisely what the above operator Fourier transform extracts. The transition weights γ(ν) satisfy the relation γ(ν)e^βν = γ(−ν), whereas their exact form depends on the relaxation dynamics of the bath (see section ‘Comparison with physically derived master equations’, Methods). The transition weights γ(ν) favour cooling (where energy decreases, ν < 0) and exponentially suppress heating (where energy increases, ν > 0). For diagonal input states and jumps A that preserve diagonal states, taking γ(ν) = 1/(1 + e^βν) or $\gamma (\nu )={{\rm{e}}}^{-\max (\beta \nu ,0)}$, the Davies generator recovers the classical Glauber or Metropolis dynamics, respectively.

Issues with the Davies generator and previous works

The Davies generator fails to be the definite quantum analogue of Glauber dynamics: while satisfying quantum detailed balance and thus enjoying the desirable mathematical properties, it cannot be implemented efficiently for general many-body systems. The crucial problem is the infinite time integral (equation (5)), which is necessary, because of energy–time uncertainties, for identifying the exact energy difference between eigenstates. The derivation of the Davies generator intends to capture an infinite-time weak-coupling system–bath dynamics and, therefore, has been considered only physical for models in which the integration time can be effectively reduced because of large structured gaps in the spectrum (such as single-particle or commuting Hamiltonians). Although it is possible to truncate the time integral, all existing schemes^18,19,20,21 fail to satisfy detailed balance, and some also lead to highly nonlocal Lindblad operators.

Our new Lindbladian satisfying quantum detailed balance

Now, we construct our Lindbladian that generates a continuous-time quantum Markov chain featuring the key properties of Glauber and Metropolis dynamics: detailed balance and locality. We can view our proposed Lindblad operators as smoothened versions of those from Davies’s ${L}_{\nu }=\sqrt{\gamma (\nu )}{A}_{\nu }$ (equation (5)).

The key idea is to replace the nonphysical transition amplitudes A_ν, which refer to exact energies, with a smooth quantum operator Fourier transform,

$$\begin{array}{l}\hat{A}(\omega )\,:= \,\sqrt{\frac{\beta }{\sqrt{2{\rm{\pi }}}}}\sum _{\nu \in B(H)}{{\rm{e}}}^{-{\beta }^{2}{(\omega -\nu )}^{2}/4}{A}_{\nu }\\ \,\,\,=\,\sqrt{\frac{\beta }{{\rm{\pi }}\sqrt{2{\rm{\pi }}}}}{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}{{\rm{e}}}^{-{t}^{2}}{{\rm{e}}}^{-{\rm{i}}\beta \omega t}{{\rm{e}}}^{{\rm{i}}\beta Ht}A{{\rm{e}}}^{-{\rm{i}}\beta Ht}{\rm{d}}t,\end{array}$$

(6)

which enjoys the sought properties. First, as $\widehat{A}(\omega )$ is given by a Gaussian-damped integral of Heisenberg evolutions of A, it approximately preserves the locality of the jump A and the Hamiltonian H up to a radius β (by Lieb–Robinson bounds, see equation (21) and ref. ³⁰). This (quasi-)locality of $\widehat{A}(\omega )$ reflects that detailed balance only depends on the energy difference before and after the jump rather than global energies. Second, the operator $\widehat{A}(\omega )$ can be naturally implemented by a phase estimation for operators that selects the transition amplitudes of A that change the energy by roughly ω (Fig. 3). Technically, equation (6) corresponds to a new variant of quantum phase estimation with a guaranteed Gaussian error profile, which has already found independent applications³¹.

Our Lindbladian ${\mathcal{L}}={\sum }_{a\in A}{{\mathcal{L}}}^{a}$ consists of the following carefully constructed Lindbladian terms ${{\mathcal{L}}}^{a}$ associated with the corresponding Hermitian jump operator A^a:

$$\begin{array}{l}{{\mathcal{L}}}^{a}[\rho ]\,:= \,\mathop{\underbrace{-i[{C}^{a},\rho ]}}\limits_{ \mbox{``} {\rm{coherent}}\mbox{''}}+\mathop{\underbrace{{\int }_{-\infty }^{\infty }\gamma (\omega ){\widehat{A}}^{a}(\omega )\rho {\widehat{A}}^{a}{(\omega )}^{\dagger }{\rm{d}}\omega }}\limits_{ \mbox{``} {\rm{transition}}\mbox{''}}-\mathop{\underbrace{\frac{1}{2}({D}^{a}\rho +\rho {D}^{a})}}\limits_{ \mbox{``} {\rm{decay}}\mbox{''}},\\ \text{where}\,\gamma (\omega )\,:= \,\exp \left(-\max (\beta \omega +\frac{1}{2},0)\right),\end{array}$$

(7)

$$\begin{array}{l}{D}^{a}:= {\int }_{-\infty }^{\infty }\gamma (\omega ){\widehat{A}}^{a}{(\omega )}^{\dagger }{\widehat{A}}^{a}(\omega ){\rm{d}}\omega ,\,\text{and}\\ {C}^{a}:= {\int }_{-\infty }^{\infty }\frac{i}{\sinh (2{\rm{\pi }}t)}({e}^{i\beta Ht}{D}^{a}{e}^{-i\beta Ht}-{D}^{a}){\rm{d}}t.\end{array}$$

(8)

A curious feature of Gaussian energy uncertainty is that the above shifted Metropolis transition weights γ(ω) achieve exact detailed balance for the transition part (equation (15)) despite the presence of energy uncertainty—generalizing the analogous classical case³². To fully achieve detailed balance, an extra coherent term C^a is necessary to accompany the decay term D^a, if it does not commute with the system Hamiltonian H. It is intriguing that the physical derivation of the Davies generator (equation (5)) also introduces a Lamb-shift term that somewhat resembles our coherent term. However, that Lamb-shift term does not play a part in the stationarity of the Gibbs state because Davies’s decay term already satisfies detailed balance as it commutes with H (Methods, Lemma 1).

Finally, note that by directly working with Lindbladians, similarly to ref. ²⁰, we circumvent difficulties around implementing the usual MCMC reject step that arose in earlier approaches because of the no-cloning theorem¹⁸. See the Methods for more details about our construction and its implementation.

Simulation costs and guarantees

As inputs to our algorithm, we require access to a set of Hermitian jumps {A^a: a ∈ A}, which is analogous to the symmetric transition rate matrix in the classical case. The jump operators A^a can be as simple as local Pauli operators, but more complex jumps may give faster mixing (for example, cluster updates in Ising models). Our construction can be extended to non-Hermitian jumps, as long as their adjoints are contained in the set³³, but we restrict to Hermitian jumps for simplicity throughout this paper. We quantify the algorithmic cost by the required Hamiltonian simulation time, a subroutine that has been thoroughly studied. The $\widetilde{{\mathcal{O}}}(\cdot )$ notation absorbs logarithmic dependencies on the evolution time t, inverse temperature β, the strength of the Hamiltonian ∥H∥, the number of qubits n, the inverse trace distance error 1/ϵ and the number of jump operators $|A|$.

Theorem 1: efficiency, locality and exact detailed balance

For any β > 0, Hamiltonian H on n-qubits, and self-adjoint jump operators A^a = A^a† normalized by $\parallel \sum _{a}{({A}^{a})}^{2}\parallel \le 1$, the Lindbladian ${\mathcal{L}}=\sum _{a\in A}{{\mathcal{L}}}^{a}$ satisfies the following:

• Each ${{\mathcal{L}}}^{a}$ satisfies quantum detailed balance (equation (4)), and the Gibbs state is stationary ${{\mathcal{L}}}^{a}[\,{\rho }_{\beta }]=0$.

• For any time t ≥ 1, the Lindbladian evolution ${{\rm{e}}}^{{\mathcal{L}}t}$ can be implemented up to ϵ-error (in diamond distance), using

$$\pm \widetilde{{\mathcal{O}}}(\beta \cdot t)\,\,\text{Hamiltonian simulation time},$$

$\widetilde{{\mathcal{O}}}(t)$ applications of the block encoding (and its inverse) for the jumps $\sum _{a\in A}| a\rangle \otimes {A}^{a}$, $\widetilde{{\mathcal{O}}}(t)$ two-qubit gates and (apart from those in the block encoding and the Hamiltonian simulation) $\widetilde{{\mathcal{O}}}(1)$ ancillas.

• If the Hamiltonian H is geometrically local, then each summand ${{\mathcal{L}}}^{a}$ is well-approximated by a local Lindbladian acting on $\widetilde{{\mathcal{O}}}(\beta )$-radius ball surrounding the jump operator A^a.

• If the set of jump operators {A^a} has no nontrivial invariant subspace (for example, all single-site Pauli X and Z operators), then the Gibbs state ρ_β is the unique fixed point.

Detailed balance has been a central feature in the study of Markov chain mixing time. The lack of exact quantum detailed balance, mainly rooted in imprecise energy measurements, has been an obstacle for noncommuting Hamiltonians. Here, achieving exact quantum detailed balance makes a formal connection to the existing framework and enables clean mathematical reasoning, as demonstrated by a series of recent works^34,35,36 that used our Lindbladian to prove rapid mixing for lattice Hamiltonians at high temperatures.

Detailed balance also enables mapping our Lindbladian to a frustration-free parent Hamiltonian ${\mathcal{H}}$ acting on a doubled, purifying system^19,20,37, whose zero-eigenstate is a certain purified Gibbs state $| \sqrt{{\rho }_{\beta }}\rangle$ (known as the thermal field double in high-energy physics) (Fig. 2). This gives a coherent Gibbs sampler, in which we prepare the purified Gibbs state by a natural adiabatic path parametrized by the inverse temperature β. Of course, the purified Gibbs state reduces to the (mixed) Gibbs state once we trace out the purifying replica system, but the purification opens doors to advanced quantum algorithmic tools, such as verification (for example, through the swap test or measuring the energy of the parent Hamiltonian) or faster mean estimation³⁸.

Fig. 2: Quasi-locality of the maps.

For a lattice Hamiltonian, our Lindbladian is a sum over quasi-local terms ${{\mathcal{L}}}^{a}$ surrounding each jump A^a with radius approximately β, with an exponentially decaying tail controlled by the Lieb–Robinson bounds. This locality remains to hold for the parent Hamiltonian ${\mathcal{H}}=\sum _{a}{{\mathcal{H}}}^{a}$ of the purified Gibbs state $| \sqrt{{\rho }_{\beta }}\rangle$ defined on two copies of the system.

The key to the efficiency of our proposal lies in the revelation that detailed balance generally does not require high-precision metrology. The Hamiltonian simulation time ~β (suppressing other dependencies) used in our Lindbladian is not enough for any quantum algorithm to learn energies beyond precision ~1/β. Yet, it is sufficient for coherently modulating the transition amplitudes according to the quantum detailed balance condition (equation (4)), leading to improved implementation methods that are analogous to the exponentially more precise quantum linear equation solvers³⁹, which circumvent the metrology step of the original phase-estimation-based proposal⁴⁰.

We can further exploit that the modulation depends only on the energy differences before and after jumps. Thus, for geometrically local Hamiltonians, our Lindbladian diagnoses energy differences by (quasi-)localized Hamiltonian patches of radius ~β, which, in turn, could further improve the efficiency of implementation by truncating faraway Hamiltonian terms. For classical or commuting Hamiltonians, the energy difference depends only on a strictly local neighbourhood, regardless of β, which simplifies the corresponding Gibbs samplers. This aspect of locality contrasts with earlier quantum methods, which typically involve global energy measurements, for example, through phase estimation.

To implement our synthetic Lindbladian, we also develop high-precision Lindbladian simulation algorithms (Theorem 2), giving an explicit quantum Gibbs sampler with provable accuracy guarantees. However, from a practical perspective, we expect there to be opportunities for heuristic optimization of the actual circuit implementation, for example, by simulating the Lindblad dynamics to low orders (see Fig. 3 for a first-order implementation), truncating the Lindbladian operator depending on the state, or tuning the discretization parameters.

Fig. 3: A circuit implementation of the Lindbladian evolution to first-order with error $\boldsymbol{\mathcal{O}}({{\boldsymbol{\delta }}}^{{\bf{2}}})$.

a, Circuit implementation for the isometry U^a representing the discrete quantum operator Fourier transform ${\widehat{A}}^{a}(\bar{\omega })$, evaluated on input $| \bar{0}\rangle \otimes | \psi \rangle$, in which we assumed for simplicity that A^a is a Hermitian unitary. The subroutines include Gaussian filter state preparation $| \,f\rangle$ (equation (6)), controlled-Hamiltonian simulation and quantum Fourier transform. The label $\bar{\omega }$ is a rough estimate of the energy difference before and after A^a (that is, the Bohr frequency). b, Quantum circuit implementation of an approximate δ-time step corresponding to a single jump operator A^a. This weak-measurement scheme uses the quantum operator Fourier transform isometry U^a and its inverse, a controlled single-qubit rotation ${Y}_{\theta }=\left[\begin{array}{l}\sqrt{1-\theta }\,-\sqrt{\theta }\\ \sqrt{\theta }\quad \quad \sqrt{1-\theta }\end{array}\right]$, Hamiltonian simulation of the coherent part C^a and ancilla resets. To simulate the full Lindbladian, this circuit should be repeatedly applied with the jump operator A^a for a uniformly random a ∈ A. We can implement a block encoding of C^a by a linear combination of unitaries¹¹ according to equation (8). Alternatively, ignoring the C^a term, detailed balance becomes approximate, with error controlled by the precision of energy difference estimation²².

Mixing times

We have shown that given a local Hamiltonian H, our new construction allows for a step-wise efficient quantum MCMC algorithm. However, just as in the classical case, the mixing times can depend crucially on the specific problem at hand. Although large-scale numerical exploration of mixing times for strongly interacting systems seems to be intractable without fault-tolerant quantum computers, we can exactly diagonalize the Lindbladian for small systems, which can, for example, reveal interesting phenomena related to the choice of jump operators⁴¹. For illustration, we choose to look at the transverse field Ising model ${H}_{{\rm{TFI}}}=-\,\sum _{j}{Z}_{j}{Z}_{j+1}+\lambda {X}_{j}$ and the XXZ model ${H}_{{\rm{XXZ}}}=\sum _{j}{X}_{j}{X}_{j+1}+{Y}_{j}{Y}_{j+1}+\gamma {Z}_{j}{Z}_{j+1}$ on a ring. Here, X, Y and Z are the single-site Pauli operators, and λ and γ are the parameters of the models. Neither model exhibits a phase transition in its thermal states, as it is on a one-dimensional lattice⁴². Hence, we expect a system-size-independent gap in the spectrum of the Lindbladian at all temperatures (when each jump is normalized to have unit strength ∥A^a∥ = 1). By contrast, the models show quantum phase transitions in their ground states.

In Fig. 4, we simulate the two models for system size n = 8 by exactly diagonalizing the Lindbladian in its matrix representation (diagonalizing larger systems is prohibitive because the number of matrix elements grows as 2⁴ⁿ). At high temperatures (β ≪ 1), we observe a large spectral gap in all parameter regimes for both models with the Lindbladian constructed from local Pauli X, Y and Z operators only, in agreement with the recent theoretical analysis^34,35. At low temperatures, both models exhibit a critical slowdown at certain values of the parameters λ and γ. However, adding a global jump X^⊗n = ∏_jX_j lifts the critical slowdown in the ferromagnetic phases, while adding two-body terms {X_jX_j+1} lifts the gap in the antiferromagnetic phase of the XXZ model. This is consistent with classical Glauber dynamics, in which slow mixing of the two-dimensional Ising model at low temperatures can be resolved by applying global flips. We take this as positive preliminary evidence that our quantum Gibbs samplers are consistent and can be expected to allow for model-specific heuristic improvements, similar to classical MCMC.

Fig. 4: Lindbladian gaps as a function of parameters of two spin chain Hamiltonians.

a, The transverse field Ising model with periodic boundary conditions on n = 8 qubits at various temperatures, with single Pauli jumps. As expected, the Lindbladian is gapped at high temperatures throughout the phase diagram, and the gap closes exponentially in β in the ferromagnetic regime (which is well-known in the λ = 0 case from the classical literature). The inset in the top outer corner shows the same gaps on a logarithmic scale. Vertical dotted lines highlight the parameter values at which a quantum phase transition happens in the ground state of the models, and the phases are labelled above. The dashed line shows that, even for the lowest temperature β = 5, the Lindbladian gap can be opened by including global spin-flip X^⊗n. b, The Heisenberg XXZ model with periodic boundary conditions on n = 8 qubits at a low temperature β = 5, with (i) only local Pauli jumps, (ii) local Pauli jumps and a global Pauli X^⊗n jump and (iii) local Pauli jumps and the nearest-neighbour Pauli XX jumps. The dynamics can freeze with single-site jumps, but adding global and nearest-neighbour Pauli XX jumps speeds up mixing by opening the gap.

Future directions and applications

A central challenge in quantum simulation is the preparation of thermal states. For this purpose, we have constructed a physically inspired, algorithmically efficient and theoretically clean Lindbladian that matches the celebrated features of classical MCMC algorithms. Our new construction opens up several avenues for further studies:

• A provable quantum simulation algorithm: the study of spectral gaps and mixing times of classical Markov chains has been a fruitful and persistent subject^28,43,44. Our construction gives a well-defined, quasi-local map with exact detailed balance and invites a systematic study for noncommuting Hamiltonians. Recently, rapid mixing at perturbative regimes or high temperatures has been proved^34,35,36 using our Lindbladians; through a dedicated mixing time analysis, refs. ^34,45 show that closely related cooling processes at low temperatures are universal for quantum computation, suggesting that simulating low-temperature cooling dynamics can be classically hard. Ultimately, directly showing that physically relevant systems (for example, the Fermi–Hubbard model) mix rapidly at classically intractable regimes gives quantitative evidence of practical quantum advantage in quantum simulation.

• A dynamical angle to quantum Gibbs states: the mixing time of a Lindbladian is naturally connected to the correlation decay and complexity of quantum Gibbs states. It was a celebrated result in classical Gibbs sampling that rapid mixing is deeply tied with the decay of Gibbs correlation²⁸, and only recently has the equivalence been extended to some commuting quantum Hamiltonians^46,47. Our newly defined map gives a dynamic angle for rigorously studying intricate quantum correlations, such as the area-law of entanglement, conditional mutual information^48,49,50, and topological order in noncommuting Gibbs states.

• A self-contained toy model of open-system thermodynamics: as a quantum analogue of Glauber dynamics, the qualitative parallel to physical master equations poses our construction as a succinct toy model of open-system physics. Thus, in cases in which we cannot prepare the Gibbs state due to a mixing time slowdown, we may still uncover properties of new thermal phases. Our Lindbladian provides a theoretical starting point for studying metastability, energy landscapes, quantum spin glasses and self-correcting quantum memories, whose precise formulation for noncommuting Hamiltonians has also been lacking. The explicit form of our Lindbladian also enables direct numerical studies regarding the above notions.

In light of the key theoretical and empirical success of MCMC methods in classical computing, we expect quantum Gibbs sampling to serve similarly important roles in quantum computing. Given the current scepticism about the practical applicability of quantum computers, we anticipate our results will initiate a new set of diverse research directions covering theory, experiment, numerics and applications.

Methods

We provide the explicit form of our construction and explain how exact detailed balance can be achieved through the key algorithmic subroutines and circuits.

Constructing the full Lindbladian

Recall that our main result considers the following Lindbladian in the Schrödinger picture:

$$\begin{array}{c}{\mathcal{L}}[\cdot ]\,:= \,\sum _{a\in A}\,\mathop{\underbrace{-i[{C}^{a},\cdot ]}}\limits_{ \mbox{``} {\rm{c}}{\rm{o}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{e}}{\rm{n}}{\rm{t}}\mbox{''}}\,\,+\mathop{\overbrace{{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}\gamma (\omega )\left(\mathop{\underbrace{{\hat{A}}^{a}(\omega )[\cdot ]{\hat{A}}^{a}{(\omega )}^{\dagger }}}\limits_{ \mbox{``} {\rm{t}}{\rm{r}}{\rm{a}}{\rm{n}}{\rm{s}}{\rm{i}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}\mbox{''}}-\mathop{\underbrace{\frac{1}{2}\{{\hat{A}}^{a}{(\omega )}^{\dagger }{\hat{A}}^{a}(\omega ),\cdot \}}}\limits_{ \mbox{``} {\rm{d}}{\rm{e}}{\rm{c}}{\rm{a}}{\rm{y}}\mbox{''}}\right){\rm{d}}\omega }}\limits^{{\rm{d}}{\rm{i}}{\rm{s}}{\rm{s}}{\rm{i}}{\rm{p}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{v}}{\rm{e}}\,{\rm{p}}{\rm{a}}{\rm{r}}{\rm{t}}}\,.\end{array}$$

(9)

Roughly, it resembles the Davies generator (equation (5)) but is carefully modified to maintain quantum detailed balance while ensuring locality and efficiency. We begin by reviewing the detailed balance condition for the Davies generator to illustrate the key ingredients and differences in our construction. We may regroup the above according to the jumps ${\mathcal{L}}={\sum }_{a\in A}{{\mathcal{L}}}^{a}$, and study each term individually, so in what follows we drop the label a for simplicity, that is, substitute A^a ← A.

Detailed balance of the Davies generator

For a Hermitian jump A, recall that the Davies generator ${{\mathcal{L}}}_{Davies}[\cdot ]\,:= \,\sum _{\nu }\gamma (\nu ){A}_{\nu }[\cdot ]{A}_{\nu }^{\dagger }-$ $\frac{\gamma (\nu )}{2}\{{A}_{\nu }^{\dagger }{A}_{\nu },\cdot \}$ satisfies the KMS-detailed balance condition (equation (4)), that is, satisfies the superoperator equation

$${\mathcal{L}}[\cdot ]={\rho }^{\frac{1}{2}}{{\mathcal{L}}}^{\dagger }[\,{\rho }^{-\frac{1}{2}}\cdot {\rho }^{-\frac{1}{2}}]{\rho }^{\frac{1}{2}},$$

(10)

where ${{\mathcal{L}}}^{\dagger }$ is the superoperator adjoint of ${\mathcal{L}}$ defined by ${\rm{Tr}}[{({\mathcal{L}}[X])}^{\dagger }Y]={\rm{Tr}}[{X}^{\dagger }{{\mathcal{L}}}^{\dagger }[Y]]$ for all X, Y. Here, detailed balance hinges on the following exact operator-valued symmetries: for all Bohr frequency ν ∈ B(H), 

$${\rho }^{-\frac{1}{2}}{A}_{\nu }{\rho }^{\frac{1}{2}}={{\rm{e}}}^{\frac{\beta \nu }{2}}{A}_{\nu }\,(\text{conjugation identity}),$$

$${A}_{-\nu }={({A}_{\nu })}^{\dagger }\,(\text{adjoint property}).$$

The conjugation identity is rooted in that A_ν are eigenoperators [H, A_ν] = νA_ν of the commutator [H, ⋅ ], highlighting a special role played by the Bohr-frequency decomposition A_ν. The adjoint property says that for a Hermitian jump A, the transition amplitudes associated with energy difference ν are paired with the reverse difference –ν, reminiscent of a Fourier transform symmetry of real functions.

As a consequence, the decay part readily satisfies detailed balance by itself because of the adjoint property, because the operator ${({A}_{\nu })}^{\dagger }{A}_{\nu }={A}_{-\nu }{A}_{\nu }$ in the decay part preserves the energies and commutes with the Hamiltonian (and hence with ρ). For the transition part,

$$\begin{array}{c}\sum _{\nu }\gamma (\nu ){\rho }^{\frac{1}{2}}{({A}_{\nu })}^{\dagger }{\rho }^{-\frac{1}{2}}\cdot {\rho }^{-\frac{1}{2}}{A}_{\nu }{\rho }^{\frac{1}{2}}\\ \,\,=\,\sum _{\nu }\gamma (\nu ){e}^{\beta \nu }{({A}_{\nu })}^{\dagger }\cdot {A}_{\nu }\,\,\,\,\,\,\text{(by the conjugation identity)}\,\\ \,\,=\,\sum _{\nu }\gamma (-\nu ){A}_{-\nu }\cdot {({A}_{-\nu })}^{\dagger }\,=\,\sum _{\nu }\gamma (\nu ){A}_{\nu }\cdot {({A}_{\nu })}^{\dagger }.\end{array}$$

(11)

The second equality uses the Kubo–Martin–Schwinger condition (in the frequency domain) for the transition weights

$$\gamma (\nu ){{\rm{e}}}^{\beta \nu }=\gamma (-\nu ).$$

(12)

The main obstacle towards a Lindbladian with exact detailed balance and efficient implementation is the lack of algorithmic access to A_ν in the presence of a dense spectrum, as approximations to A_ν can easily break the exact symmetry conditions in detailed balance^18,20,51.

Operator Fourier transforms

The key to both efficiency and exact detailed balance is a careful relaxation of A_ν that preserves some aspects of the symmetries using the operator Fourier transform damped by a Gaussian filter $f(t)\,:= \,{{\rm{e}}}^{-{\sigma }^{2}{t}^{2}}\sqrt{\sigma \sqrt{2/{\rm{\pi }}}}$ (normalized by ${\int }_{-\infty }^{\infty }{| f(t)| }^{2}{\rm{d}}t=1$) with a tunable width ∝ 1/σ (setting $\sigma =\frac{1}{\beta }$ recovers (equation (6)):

$$\begin{array}{c}\widehat{A}(\omega )\,:= \,\frac{1}{\sqrt{2{\rm{\pi }}}}{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}{{\rm{e}}}^{{\rm{i}}Ht}A{{\rm{e}}}^{-{\rm{i}}Ht}{{\rm{e}}}^{-{\rm{i}}\omega t}f(t){\rm{d}}t\\ \,\,\,=\,\frac{1}{\sqrt{{\sigma }\sqrt{2{\rm{\pi }}}}}\sum _{\nu \in B(H)}{A}_{\nu }\,{{\rm{e}}}^{-\frac{{(\omega -\nu )}^{2}}{4{{\sigma }}^{2}}}.\end{array}$$

(13)

Like A_ν in the Davies generator (equation (5)), our operator Fourier transforms $\widehat{A}(\omega )$ are labelled by energy differences, but here the parameter ω ∈ [−∞, ∞] takes continuous values, without referring to the true Bohr frequencies. When the uncertainty vanishes σ → 0, we recover $\mathop{\mathrm{lim}}\limits_{\sigma \to 0}\sqrt{\sigma \sqrt{2{\rm{\pi }}}}\widehat{A}(\omega )=\sum _{\nu \in B(H)}{\mathbb{1}}(\nu =\omega ){A}_{\nu }$; when the energy uncertainty is finite σ ≠ 0, the operator Fourier transforms still maintain exact algebraic properties crucial for detailed balance. First, conjugating with the Gibbs state preserves the form of operator Fourier transforms, albeit causing a shift and rescaling

$$\begin{array}{c}{\rho }^{-\frac{1}{2}}\widehat{A}(\omega ){\rho }^{\frac{1}{2}}\,=\,\frac{1}{\sqrt{{\sigma }\sqrt{2{\rm{\pi }}}}}\sum _{\nu \in B(H)}{A}_{\nu }{{\rm{e}}}^{\frac{\beta \nu }{2}}\exp \left(-\frac{{(\omega -\nu )}^{2}}{4{{\sigma }}^{2}}\right)\,\,\,\,({\rm{b}}{\rm{y}}\,(13)\,\,\text{and the conjugation identity)}\,\\ \,\,\,\,=\,\frac{1}{\sqrt{{\sigma }\sqrt{2{\rm{\pi }}}}}\sum _{\nu \in B(H)}{{\rm{e}}}^{\frac{\beta \omega }{2}+\frac{{{\sigma }}^{2}{\beta }^{2}}{4}}{A}_{\nu }\exp \left(-\frac{{(\omega -\nu +{{\sigma }}^{2}\beta )}^{2}}{4{{\sigma }}^{2}}\right)\,\,\,\,\,\text{(by completing the square)}\,\\ \,\,\,\,=\,{{\rm{e}}}^{\frac{\beta \omega }{2}+\frac{{{\sigma }}^{2}{\beta }^{2}}{4}}\widehat{A}(\omega +{{\sigma }}^{2}\beta ),\,\,\,\,\,\text{(by shift-rescale symmetry)}\end{array}$$

which reflects the fact that multiplying a Gaussian distribution by an exponential weight must shift the mean ${{\rm{e}}}^{{(x-a)}^{2}-2bx}={{\rm{e}}}^{{(x-a-b)}^{2}-2ab-{b}^{2}}$. Second, even though the operator Fourier transform is a linear combination of different Bohr frequencies, the adjoint property holds exactly when A = A^† as

$$\begin{array}{l}\widehat{A}(\,-\,\omega )\,=\,\frac{1}{\sqrt{2{\rm{\pi }}}}{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}{{\rm{e}}}^{{\rm{i}}Ht}A{{\rm{e}}}^{-{\rm{i}}Ht}{{\rm{e}}}^{{\rm{i}}\omega t}f(t){\rm{d}}t\\ \,\,\,\,=\,\frac{1}{\sqrt{2{\rm{\pi }}}}{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}{{\rm{e}}}^{{\rm{i}}Ht}{A}^{\dagger }{{\rm{e}}}^{-{\rm{i}}Ht}{({{\rm{e}}}^{-{\rm{i}}\omega t}f(t))}^{\ast }{\rm{d}}t=\widehat{A}{(\omega )}^{\dagger }\,\,\,\,({\rm{s}}{\rm{i}}{\rm{n}}{\rm{c}}{\rm{e}}\,\,{f}^{\ast }(t)=\,f(t)).\end{array}$$

The above two exact symmetries appear to be absent in previous approaches that attempted to directly measure energy differences. Now we show that quantum detailed balance is a consequence of these exact symmetries and related algebraic properties of the Gaussian uncertainty in operator Fourier transforms, which hold exactly despite not measuring energies to high precision.

Exact detailed balance with finite uncertainty

Although the operator Fourier transform does not yield an exact representation of A_ν, the uncertainty has a specific structure with distinctive symmetries—arising from the interplay between the Gaussian weighing and the exponential form of the Boltzmann factors—and consequently can be exactly compensated by an appropriate shift in the transition weights. In the following, we prove that the transition part (equation (9)) satisfies detailed balance (equation (10)): ${\mathcal{T}}\,[\cdot ]={\rho }^{\frac{1}{2}}{{\mathcal{T}}}^{\dagger }[\,{\rho }^{-\frac{1}{2}}\cdot {\rho }^{-\frac{1}{2}}]{\rho }^{\frac{1}{2}}$. The above shift-rescale and adjoint symmetries yield

$$\begin{array}{c}{\rho }^{\frac{1}{2}}{{\mathcal{T}}}^{\dagger }[{\rho }^{-\frac{1}{2}}\cdot {\rho }^{-\frac{1}{2}}]{\rho }^{\frac{1}{2}}\\ \,\,=\,{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}\gamma (\omega ){\left({\rho }^{-\frac{1}{2}}\widehat{A}(\omega ){\rho }^{\frac{1}{2}}\right)}^{\dagger }[\cdot ]{\rho }^{-\frac{1}{2}}\widehat{A}(\omega ){\rho }^{\frac{1}{2}}{\rm{d}}\omega \,\,\,\,\text{(by definition)}\\ \,\,=\,{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}\gamma (\omega ){{\rm{e}}}^{\beta \omega +\frac{{{\sigma }}^{2}{\beta }^{2}}{2}}\widehat{A}(-\omega -{{\sigma }}^{2}\beta )[\cdot ]\widehat{A}{(-\omega -{{\sigma }}^{2}\beta )}^{\dagger }{\rm{d}}\omega ,\end{array}$$

(14)

showing that we may compensate for the uncertainty by imposing a shifted KMS condition (equation (12))

$$\begin{array}{c}\gamma (\omega ){{\rm{e}}}^{\beta \omega +\frac{{{\sigma }}^{2}{\beta }^{2}}{2}}\,=\,\gamma (\,-\,\omega -{{\sigma }}^{2}{\beta }^{2}),\,\,\,\text{such that}\\ \,\,(14)\,=\,{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}\gamma (\,-\,\omega -{{\sigma }}^{2}\beta )\widehat{A}(\,-\,\omega -{{\sigma }}^{2}\beta )[\cdot ]\widehat{A}{(-\omega -{{\sigma }}^{2}\beta )}^{\dagger }{\rm{d}}\omega ={\mathcal{T}}\,[\cdot ].\end{array}$$

(15)

Therefore, we can take any transition weight function satisfying the KMS condition γ₀(ν)e^βν = γ₀(−ν) and pretend we underestimated the Bohr frequency by $\frac{{\sigma }^{2}\beta }{2}$, that is, substitute $\nu \leftarrow {\omega }_{+}:= \omega +\frac{{\sigma }^{2}\beta }{2}$ to obtain our shifted γ(ω). The canonical examples include the Metropolis and Glauber weights

$${\gamma }^{M}(\omega )\,:= \,\exp (-\beta \,\max ({\omega }_{+},0)).\,(\text{shifted Metropolis})$$

(16)

$${\gamma }^{G}(\omega )\,:= \,\frac{1}{1+{{\rm{e}}}^{\beta {\omega }_{+}}}.\,(\text{shifted Glauber})$$

(17)

See ref. ³³ for the original step-by-step detailed derivation; the above streamlined derivation draws partly from subsequent simplification and development (Lemma 7.1 of ref. ⁵² and Lemma IX.2 of ref. ⁴⁹).

Under the natural normalization γ(ω) ∈ [0, 1], the shifted symmetry (equation (15)) implies a low transition rate $\gamma (0)\le \exp (-{\sigma }^{2}{\beta }^{2}/2)$ around ω = 0. Therefore, to avoid unnecessary idling of the process, it is imperative to choose an uncertainty σ ≤ 1/β not exceeding the temperature; in the main text, we have simply set $\sigma =\frac{1}{\beta }$. This is why implementing a single step uses ~β Hamiltonian simulation time in our construction. By contrast, for classical or commuting systems with gapped periodic spectrum, the uncertainty can be discretized, and there is no need to scale the Hamiltonian simulation time with β.

Achieving full detailed balance by tuning the coherent term

Even if the transition part (equation (9)) satisfies detailed balance exactly, the decay part D of the Lindbladian may still break detailed balance when it does not commute with the Hamiltonian H.

A second insight of our construction is to prescribe and efficiently implement a dedicated coherent term C that perfectly cancels out the deviation from quantum detailed balance in a uniquely quantum way, as shown by the following lemma.

Lemma 1: prescribing the coherent term

(Lemma II.1 of ref. ³³) For any full-rank density operator $0\,\prec \,\rho \in {{\mathbb{C}}}^{d\times d}$ and Hermitian operator $D\in {{\mathbb{C}}}^{d\times d}$, there is a unique Hermitian operator $C\in {{\mathbb{C}}}^{d\times d}$ (up to adding any scalar multiples of the identity I) such that the superoperator

$$-\frac{1}{2}\{D,\cdot \}-{\rm{i}}[C,\cdot ]$$

(18)

satisfies ρ-detailed balance. For a Gibbs state $\rho \propto \exp (-\beta H)$, we can express C as

$$\begin{array}{c}\,C=\sum _{\nu \in B(H)}\frac{{\rm{i}}}{2}\,\tanh \left(\frac{\beta \nu }{4}\right){D}_{\nu }\\ \text{where}\,\,{D}_{\nu }:= \sum _{{E}_{i}-{E}_{j}=\nu }|{\psi }_{i}\rangle \langle {\psi }_{i}|D|{\psi }_{j}\rangle \langle {\psi }_{j}|.\end{array}$$

(19)

Proof

The detailed balance condition is equivalent to the following matrix being self-adjoint:

$$\begin{array}{c}{\rho }^{-\frac{1}{4}}\left(-\frac{D}{2}-{\rm{i}}C\right){\rho }^{\frac{1}{4}}\,=\,\frac{1}{2}\,{\rho }^{-\frac{1}{4}}\left(\sum _{\nu \in B(H)}{D}_{\nu }\left(\tanh \left(\frac{\beta \nu }{4}\right)-1\right)\right){\rho }^{\frac{1}{4}}\\ \,\,\,\,\,\,\,=\,-\frac{1}{2}\sum _{\nu \in B(H)}\frac{{D}_{\nu }}{\cosh \left(\frac{\beta \nu }{4}\right)}.\end{array}$$

This is self-adjoint because ${({D}_{\nu })}^{\dagger }={D}_{-\nu }$, $\cosh (x)=\cosh (-x)$, and B(H) = −B(H). □

See also ref. ⁵³ on prescribing fixed points in general, without necessarily assuming detailed balance. The coherent term C is a Hermitian matrix obtained by reweighing the given D operator with the profile ${\rm{i}}\,\tanh (\beta \nu /4)/2$ on each Bohr frequency component D_ν. The coherent term is completely determined by ρ and D (ref. ⁵⁴), and we found a particularly simple and useful closed-form representation in the time domain (recall ${\rm{sinhc}}(x)\,:= \,\sinh (x)/x$ for real x ≠ 0)

$$C={\int }_{-\infty }^{\infty }\frac{{\rm{i}}}{\sinh (2{\rm{\pi }}t)}({{\rm{e}}}^{{\rm{i}}\beta Ht}D{{\rm{e}}}^{-{\rm{i}}\beta Ht}-D){\rm{d}}t=\frac{\beta }{2{\rm{\pi }}}{\int }_{-\infty }^{\infty }\frac{-1}{{\rm{sinhc}}(2{\rm{\pi }}t)}{\int }_{0}^{1}{{\rm{e}}}^{{\rm{i}}s\beta Ht}[H,D]{{\rm{e}}}^{-{\rm{i}}s\beta Ht}{\rm{d}}s{\rm{d}}t,$$

(20)

which in turn implies the algorithmic efficiency and quasi-locality of the coherent term. The above integral form and the exponential decay in t allow using the linear combination of unitaries technique to implement a block encoding of C (refs. ^33,55) by truncating and discretizing the time integral. For the Davies generator (equation (5)), we have that [H, D] = 0. Therefore, the dissipative part readily satisfies detailed balance, and D_ν = 0 for all ν ≠ 0, implying that the above coherent term prescription simply vanishes.

Quasi-locality

A salient feature of our construction is the quasi-locality of the Lindbladian terms, inherited from the operator Fourier transform in systems that feature a Lieb–Robinson bound. For all $\omega \in {\mathbb{R}}$ and geometrically local jump A, truncating the Hamiltonian to a local Hamiltonian patch H_ℓ within distance ℓ from the jump yields an error

$$\left\Vert {\widehat{A}}_{(H)}(\omega )-{\widehat{A}}_{({H}_{{\ell }})}(\omega )\right\Vert \le \frac{1}{\sqrt{2{\rm{\pi }}}}{\int }_{-\infty }^{\infty }\left\Vert {{\rm{e}}}^{{\rm{i}}Ht}A{{\rm{e}}}^{-{\rm{i}}Ht}-{{\rm{e}}}^{{\rm{i}}{H}_{{\ell }}t}A{{\rm{e}}}^{-{\rm{i}}{H}_{{\ell }}t}\right\Vert |\,f(t)|{\rm{d}}t.$$

(21)

For a wide variety of local Hamiltonian systems, off-the-shelf Lieb–Robinson bounds³⁰ for the Heisenberg dynamics e^iHtAe^−iHt state that the integrand in equation (21) exponentially reduces with the distance ℓ but degrades with the evolution time t. As the Gaussian weight function f(t) of equation (13) effectively dampens the integral to small values of t ~ 1/σ, the operator ${\widehat{A}}_{(H)}(\omega )$ is well-approximated by a Hamiltonian patch with radius scaling with the energy uncertainty σ, which can be independent of the system size. See, for example, Appendix A of ref. ⁴⁹ for a quantitative estimate for both the transition and coherent parts. Of course, after truncation, the resulting strictly local Lindbladian may no longer satisfy exact quantum detailed balance.

Fixed point and its uniqueness

The detailed balance condition (equation (10)) directly implies that the Gibbs state ρ is a fixed point for the Lindbladian ${\mathcal{L}}$.

$${\mathcal{L}}[\,\rho ]={\rho }^{\frac{1}{2}}{{\mathcal{L}}}^{\dagger }\,[{\rho }^{-\frac{1}{2}}({\rho }^{\frac{1}{2}}){\rho }^{-\frac{1}{2}}]\,{\rho }^{\frac{1}{2}}={\rho }^{\frac{1}{2}}{{\mathcal{L}}}^{\dagger }[I]{\rho }^{\frac{1}{2}}=0,$$

where the last equality used the trace-preservation property ${{\mathcal{L}}}^{\dagger }[I]=0$ of Lindbladians.

We now explain why the Gibbs state is the unique fixed point whenever the set of jump operators has no (nontrivial) invariant subspaces, which holds, for example, when the jumps include all single-site Pauli X and Z operators. Decomposing each jump operator by the operator Fourier transform ${A}^{a}\propto {\int }_{-\infty }^{\infty }{\widehat{A}}^{a}(\omega ){\rm{d}}\omega$ cannot create new invariant subspaces; likewise, multiplying by strictly positive transition weights γ(ω) cannot. Thus, the resulting set of Lindblad operators $\sqrt{\gamma (\omega )}{\widehat{A}}^{a}(\omega )$ has no invariant subspaces, which is known⁵⁶ to imply the uniqueness of the fixed point. But this uniqueness argument says little about the quantitative convergence rate, and the mixing times can depend on the particular Hamiltonian.

A single-qubit example

Let us illustrate the details of our Lindbladian construction through a pedagogical example. Consider a single-qubit Hamiltonian with a single jump:

$$\begin{array}{l}H=Z=\left[\begin{array}{ll}1 & 0\\ 0 & -1\end{array}\right]=| 0\rangle \langle 0| -| 1\rangle \langle 1| \quad \text{and}\\ A=X=\left[\begin{array}{ll}0 & 1\\ 1 & 0\end{array}\right]=| 1\rangle \langle 0| +| 0\rangle \langle 1| .\end{array}$$

(22)

The eigenvalues of the Hamiltonian are ±1, and the Bohr frequencies, their differences, are B(H) = {2, 0, −2}. We can decompose the jump by the Bohr frequencies, into the ν = 2 and ν = −2 components as follows:

$${A}_{2}=| 0\rangle \langle 1| \,\,\text{and}\,\,{A}_{-2}=| 1\rangle \langle 0| .$$

(23)

We can then directly obtain the Davies generator for any transition weights satisfying γ(2)e^2β = γ(−2). By contrast, our operator Fourier transform yields

$$\widehat{A}(\omega )=\frac{1}{\sqrt{{\sigma }\sqrt{2{\rm{\pi }}}}}\left({{\rm{e}}}^{-\frac{{(\omega -2)}^{2}}{4{{\sigma }}^{2}}}{A}_{2}+{{\rm{e}}}^{-\frac{{(\omega +2)}^{2}}{4{{\sigma }}^{2}}}{A}_{-2}\right).$$

(24)

Thus, for all $\omega \in {\mathbb{R}}$ the resulting $\widehat{A}(\omega )$ has contributions coming from both the A₂ and A₋₂ components; these components only get separated in the σ → 0 limit. Nevertheless, detailed balance still holds at every finite uncertainty σ with suitable choices of γ(ω) (equation (15)) and the additional coherent term.

Recovering the Davies generator

As we show here, our Lindbladian exactly recovers the Davies generator in the σ → 0 limit. Moreover, the Davies generator reduces to Glauber dynamics when the Hamiltonian H is a classical function of bitstrings and the jumps A map a bitstring to another: in this case, inputs that are a probabilistic mixture of bitstrings undergo a classical Glauber dynamics (see section ‘Quantum MCMC by master equations for thermalization’).

For any bounded, continuous function γ₀ satisfying the KMS condition γ₀(ν)e^βν = γ₀(−ν), consider the shift ${\gamma }_{\beta \sigma }(\omega )\,:= \,{\gamma }_{0}\left(\omega +\frac{{\beta }^{2}{\sigma }^{2}}{2}\right),$ then the transition part (similarly for the decay term) converges to

$$\mathop{\mathrm{lim}}\limits_{\sigma \to 0}{\int }_{-\infty }^{\infty }{\gamma }_{\beta \sigma }(\omega )\widehat{A}(\omega )[\cdot ]\widehat{A}{(\omega )}^{\dagger }{\rm{d}}\omega =\sum _{\nu \in B(H)}{\gamma }_{0}(\nu ){A}_{\nu }[\cdot ]{({A}_{\nu })}^{\dagger }.$$

When we expand the left-hand side as a sum over ${A}_{{\nu }_{1}}[\cdot ]{({A}_{{\nu }_{2}})}^{\dagger }$, the coefficients for each ν₁, ν₂

$${\mathrm{lim}}_{\sigma \to 0}{\int }_{-\infty }^{\infty }\frac{{\gamma }_{\beta \sigma }(\omega )}{\sigma \sqrt{2{\rm{\pi }}}}{{\rm{e}}}^{-\frac{{(\omega -{\nu }_{1})}^{2}}{4{\sigma }^{2}}}{{\rm{e}}}^{-\frac{{(\omega -{\nu }_{2})}^{2}}{4{\sigma }^{2}}}{\rm{d}}\omega ={\mathbb{1}}({\nu }_{1}={\nu }_{2}){\gamma }_{0}({\nu }_{1}),$$

approach the corresponding value in the Davies generator. By contrast, the coherent term vanishes, because the decay term commutes with the Hamiltonian in the σ → 0 limit and tanh(0) = 0 (see Lemma 1).

Implementation

In this section, we give low-level quantum algorithmic implementations of our Lindbladian dynamics.

A phase estimation for operators

The most general form of the discrete operator Fourier transform is a subroutine similar to quantum phase estimation: it combines controlled Hamiltonian simulation with the quantum Fourier transform, as shown in Fig. 3a, acting on the time–frequency and system registers. To approximately implement the (continuous) operator Fourier transform, we need to discretize the time integral by introducing a time mesh $\bar{t}\in {S}_{{t}_{0}}$ and a corresponding frequency mesh $\bar{\omega }\in {S}_{{\omega }_{0}}$ for the QFT, each having M points, resulting in the discretized operation

$$\begin{array}{l}\mathop{\underbrace{{\sum }_{\bar{t}\in {S}_{{t}_{0}}}f(\bar{t})| \bar{t}\rangle \otimes A\to {\sum }_{\bar{\omega }\in {S}_{{\omega }_{0}}}| \bar{\omega }\rangle \otimes \widehat{A}(\bar{\omega })}}\limits_{{\rm{discrete}}\,{\rm{operator}}\,{\rm{Fourier}}\,{\rm{transform}}}\\ \text{where}\,\widehat{A}(\bar{\omega })\,:= \,\frac{1}{\sqrt{M}}\sum _{\bar{t}\in {S}_{{t}_{0}}}{e}^{-i\bar{\omega }\bar{t}}f(\bar{t}){e}^{iHt}A{e}^{-iHt}.\end{array}$$

(25)

As given in Corollary C.2 of ref. ²², choosing

$$\begin{array}{l}\,{\omega }_{0}=2\sigma \sqrt{\frac{2{\rm{\pi }}}{M}},\,{t}_{0}=\frac{1}{2\sigma }\sqrt{\frac{2{\rm{\pi }}}{M}},\\ {S}^{\lceil M\rfloor }:= \,\{-\lceil (M-1)/2\rceil ,\ldots ,-1,0,1,\ldots ,\lfloor (M-1)/2\rfloor \},\end{array}$$

(26)

$${\rm{and}}\quad {S}_{{\omega }_{0}}:= {\omega }_{0}\cdot {S}^{\lceil M\rfloor },\quad {S}_{{t}_{0}}:= {t}_{0}\cdot {S}^{\lceil M\rfloor },$$

(27)

the above equation (25) recovers the advertised continuous operator Fourier transform (equation (13)) in the M → ∞ limit. As f is a smooth Gaussian function, we can achieve any finite precision ϵ-approximation of the dissipative part of ${\mathcal{L}}$ with a moderately scaling dimension $M \sim {\rm{Poly}}(\parallel H\parallel ,\beta ,1/{\epsilon },\sigma +1/\sigma ,| A| )$, which requires only $\log (M)=\widetilde{{\mathcal{O}}}(1)$-many ancilla qubits, and no more than ${\mathcal{O}}(\sqrt{\log (1/{\epsilon })}\,/\sigma )$ (controlled) Hamiltonian simulation time, because of the truncation of the Gaussian tail.

A general-purpose Lindbladian simulation algorithm

Another key algorithmic component is an improved black-box Lindbladian simulation subroutine. It achieves the sought nearly linear dependence on t even for Lindbladians with high Kraus rank, as long as the Lindblad operators are given in a block-encoded form, which represents an improvement on previous Lindbladian simulation algorithms^51,57.

Definition 1: block encoding for Lindblad operators

(Definition I.2 of ref. ²²) We say that a unitary U is a block encoding for Lindblad operators {L_j: j ∈ J}, if

$$(\langle {0}^{b}|\otimes I)\cdot U\cdot (|{0}^{q}\rangle \otimes I)=\sum _{j\in J}|\,j\rangle \,\otimes {L}_{j}\,\text{for some}\,b\le q\in {\mathbb{N}}.$$

In particular, discretized operator Fourier transforms (equation (25)) naturally give block encodings of this form, where J refers to the discretized set of frequency labels ${S}_{{\omega }_{0}}$. Given a block encoding U of the Lindblad operators as above, we can directly obtain a block encoding of the decay term D by a single use of U and U^† by standard multiplication of block-encoded matrices¹¹. As a consequence of the exponentially decaying tails in the integral representation of equation (20), by using an additional ~ β-time controlled Hamiltonian simulation and the linear combination of unitaries (LCU) technique⁵⁵, we can also obtain a sufficiently good approximate block encoding of the coherent term C.

We begin with a brief analysis of the first-order simulation circuit shown in Fig. 3b. Assuming the system register is in the pure state |ψ⟩, the first three gates act as follows (for simplicity, we drop the superscripts A^a ← A):

$$\begin{array}{cc}|0\rangle |{0}^{q}\rangle |\psi \rangle & \mathop{\to }\limits^{(1)}|0\rangle \otimes \sum _{\bar{\omega }\in {S}_{{\omega }_{0}}}|\bar{\omega }\rangle \otimes \hat{A}(\bar{\omega })|\psi \rangle \\ & \mathop{\to }\limits^{(2)}\sum _{\bar{\omega }\in {S}_{{\omega }_{0}}}(\mathop{\underbrace{\sqrt{1-\delta \gamma (\bar{\omega })}}}\limits_{1-\frac{\delta \gamma (\bar{\omega })}{2}+{\mathcal{O}}({\delta }^{2})}|0\rangle +\sqrt{\delta \gamma (\bar{\omega })}|1\rangle )|\bar{\omega }\rangle \hat{A}(\bar{\omega })|\psi \rangle \\ & \mathop{\to }\limits^{(3)}|0\rangle |{0}^{q}\rangle \left({I}-\frac{\delta }{2}D\right)|\psi \rangle +|1\rangle \sum _{\bar{\omega }\in {S}_{{\omega }_{0}}}\sqrt{\delta \gamma (\bar{\omega })}|\bar{\omega }\rangle \hat{A}(\bar{\omega })|\psi \rangle -\frac{\delta }{2}|0\rangle \otimes |{0}^{q}\perp \rangle +{\mathcal{O}}({\delta }^{2}),\end{array}$$

(28)

where $| {0}^{q}\perp \rangle$ is a quantum state such that $\parallel | {0}^{q}\perp \rangle \parallel \le 1$ and $(\langle {0}^{q}| \otimes I)\cdot | {0}^{q}\perp \rangle =0$, see Theorem III.1 of ref. ²², for details.

Let $| {\psi }^{{\prime} }\rangle$ be the resulting state in equation (28), and let us denote the dissipative part ${{\mathcal{L}}}^{{\prime} }:= {{\mathcal{L}}}^{a}+{\rm{i}}[{C}^{a},\cdot ]$. Tracing out the first q + 1 qubits, we get that $| {\psi }^{{\prime} }\rangle$ is ${\mathcal{O}}({\delta }^{2})$-close to the desired state, ignoring the coherent term. We now show that

$${\parallel ({\mathcal{I}}+\delta {{\mathcal{L}}}^{{\prime} })[| \psi \rangle \langle \psi | ]-{{\rm{Tr}}}_{q+1}[| {\psi }^{{\prime} }\rangle \langle {\psi }^{{\prime} }| ]\parallel }_{1}={\mathcal{O}}({\delta }^{2})$$

(29)

by observing that

$$\begin{array}{c}{{\rm{T}}{\rm{r}}}_{q+1}[|{\psi }^{{\prime} }\rangle \langle {\psi }^{{\prime} }|]\,=\,{{\rm{T}}{\rm{r}}}_{q}[(\langle 0|\otimes I)\cdot |{\psi }^{{\prime} }\rangle \langle {\psi }^{{\prime} }|\cdot (|0\rangle \otimes I)]+{{\rm{T}}{\rm{r}}}_{q}[(\langle 1|\otimes I)\cdot |{\psi }^{{\prime} }\rangle \langle {\psi }^{{\prime} }|\cdot (|1\rangle \otimes I)]\\ \,\,\,\,\,\,=\,\left(I-\frac{\delta }{2}D\right)|\psi \rangle \langle \psi |\left(I-\frac{\delta }{2}D\right)\\ \,\,\,\,\,\,+\,\delta \sum _{\overline{\omega }\in {S}_{{\omega }_{0}}}\gamma (\overline{\omega })\hat{A}(\overline{\omega })|\psi \rangle \langle \psi |\hat{A}{(\overline{\omega })}^{\dagger }+{\mathcal{O}}({\delta }^{2})\\ \,\,\,\,\,\,=\,({\mathcal{I}}+\delta {{\mathcal{L}}}^{{\prime} })[|\psi \rangle \langle \psi |]+{\mathcal{O}}({\delta }^{2}).\end{array}$$

Convexity and the triangle inequality imply that equation (29) also holds for mixed input states. By including the δ-time Hamiltonian simulation for C^a, we get an ${\mathcal{O}}({\delta }^{2})$ approximation of δ-time evolution by ${{\mathcal{L}}}^{a}$. Once again, owing to linearity and the triangle inequality, this also implies that performing the circuit Fig. 3b for a uniformly random jump A^a we get a δ-time evolution by $\frac{1}{|A|}\sum _{a\in A}{{\mathcal{L}}}_{a}$ up to ${\mathcal{O}}({\delta }^{2})$ error in trace distance, see Corollary III.1 of ref. ²². Repeating the entire argument after replacing the jump operators A^a by A^a ⊗ I, we can see that actually the above results in a ${\mathcal{O}}({\delta }^{2})$-precise implementation of the quantum channel $\exp \left(\frac{\delta }{|A|}\sum _{a\in A}{{\mathcal{L}}}_{a}\right)$ in the completely bounded 1-1 superoperator norm, that is, the diamond-norm.

Choosing $\delta =\varTheta \left(\frac{{\epsilon }}{t}\right)$ ensures that the error in a single time step is bounded by ${\mathcal{O}}\left(\frac{{{\epsilon }}^{2}}{{t}^{2}}\right)$, and repeating the process $\varTheta \left(\frac{{t}^{2}}{{\epsilon }}\right)$-times induces an error that is bounded by ϵ for the entire time-t evolution. The complexity is then $\varTheta \left(\frac{{t}^{2}}{{\epsilon }}\right)$ times the cost of implementing the circuit in Fig. 3b.

Building on the compression techniques in ref. ⁵⁷, we can bootstrap the first-order weak-measurement circuit of Fig. 3b by observing that for very small time steps, the circuit is very close to identity. Exploiting this by effectively only using the circuit Fig. 3b in parts of the trajectories that are nontrivial, we can achieve almost linear scaling in t and polylogarithmic scaling in the desired diamond-norm accuracy of the simulation. However, we need to apply some modifications, as it turns out the original compressed measurement scheme described in ref. ⁵⁷ does not work as intended. Thus, we use a variant of the analogous measurement scheme described in ref. ⁵⁸ (see Appendix F of ref. ²², for more details). This amounts to our ultimate near-linear-time simulation result.

Theorem 2: almost linear-time Lindbladian simulation

(Theorem III.2 of ref. ²²) Suppose U is a unitary block encoding of the Lindblad operators of ${\mathcal{L}}$ as in Definition 1, and V is a block encoding of the coherent term C. Let t ≥ 1 and ϵ ≤ 1/2, then we can simulate the map ${e}^{{\mathcal{L}}t}$ to error ϵ in diamond norm using

$${\mathcal{O}}((q+\log (t/{\epsilon }))\log (t/{\epsilon }))\,\,({\rm{resettable}})\,{\rm{ancilla}}\,{\rm{qubits}},$$

(30)

$${\mathcal{O}}\left(t\frac{\log (t/{\epsilon })}{\log \,\log (t/{\epsilon })}\right)\,\,({\rm{controlled}})\,{\rm{uses}}\,{\rm{of}}\,U,V,{U}^{\dagger }\,{\rm{and}}\,{V}^{\dagger },$$

(31)

$$\,\text{and}\,\,{\mathcal{O}}(t(q+1){\rm{polylog}}(t/{\epsilon }))\,\,\,\text{other two-qubit gates},$$

(32)

where q is the number of ancilla qubits used for the block encodings.

To place the above general result in context, we give some simple bounds on the resources required for implementing our Lindbladian. For instance, if the jump operators are K different Pauli strings, then $q={\log }_{2}(K)+1$ ancilla qubits suffice for block encoding them. The overall gate complexity should be dominated by the controlled Hamiltonian simulation subroutine. Thus, we focus on estimating the required Hamiltonian simulation time per call to U and V (assuming an operator Fourier transform width $\sigma =\frac{1}{\beta }$). Under the normalization $\parallel \sum _{a\in A}{A}^{a\dagger }{A}^{a}\parallel \le 1$, the Lindblad operators from our construction (equation (7)) can be ϵ-accurately block-encoded by discretized operator Fourier transform using ${\mathcal{O}}(\beta \sqrt{\log (1/{\epsilon })})$ (controlled) Hamiltonian simulation time by truncating the Gaussian integrand in equation (6). Meanwhile, as implied by Corollary III.2 of ref. ³³, a slightly subnormalized coherent term C/α for $\alpha ={\mathcal{O}}\left(\log \left(\frac{\beta \,\parallel \,H\,\parallel }{{\epsilon }}\right)\right)$ can be ϵ-accurately block-encoded by LCU by using a mere ${\mathcal{O}}(\beta \,\log (1/{\epsilon }))$ (controlled) Hamiltonian simulation time. The extra subnormalization factor α can be absorbed into the number of uses of the block encoding V using state-of-the-art block-encoded Hamiltonian simulation techniques^4,11.

To approximate ${e}^{{\mathcal{L}}t}$ to ϵ-diamond distance, we can control the accumulation of errors by setting an increased accuracy goal of ${\rm{poly}}({\epsilon }/(t\beta \parallel H\parallel ))$ for the approximate block encodings U, V. This results in an overall ${\mathcal{O}}(\beta \sqrt{\log (t\beta \parallel H\parallel /{\epsilon })})$ and ${\mathcal{O}}(\beta {\log }^{2}(t\beta \parallel H\parallel /{\epsilon }))$ (controlled) Hamiltonian simulation time overhead for implementing sufficiently accurate block encodings U and V, respectively.

Comparison with physically derived master equations

Our synthetic Lindbladian is mainly presented as an algorithmic construction. In this section, we compare it with physically motivated master equations derived from first-principles open-system calculations. Overall, we believe that our Lindbladian can serve as a self-contained toy model of thermalization. Nevertheless, for quantitative modelling of particular system–bath interactions, the role of strict detailed balance is less clear, and it may be preferable to use a physically derived master equation.

The textbook setup in open-system thermalization^24,59 considers a system governed by a Hamiltonian H_S that couples weakly to a large thermal bath with Hamiltonian H_B, which are together governed by the total Hamiltonian ${H}_{{\rm{tot}}}={H}_{{\rm{S}}}\otimes {I}_{{\rm{B}}}+{I}_{{\rm{S}}}\otimes {H}_{{\rm{B}}}+\lambda \sum _{a\in A}{A}^{a}\otimes {B}^{a}$. The Hermitian operators {A^a: a ∈ A} act on the system (mirroring our jump operators, hence the same notation) and {B^a: a ∈ A} acts on the bath, whereas λ represents the coupling strength. Tracing out the bath, we can obtain an effective master equation governing the system dynamics under the assumptions that the thermal bath is Markovian and the coupling λ is sufficiently weak.

The aforementioned Davies generator was originally derived in the weak-coupling limit λ → 0 (relative to all other energy scales). In this limit, the rotating-wave approximation (or the secular approximation) removes the cross terms ${A}_{{\nu }_{1}}[\cdot ]{A}_{{\nu }_{2}}^{\dagger }$; this perfect isolation of Bohr frequencies causes the Davies generator to satisfy detailed balance, but at the same time, makes it implausible for many-body systems with exponentially small level spacings. To focus on issues related to detailed balance, here we studied only the simplest form of the Davies generator. In principle, different jumps A^a may have different transition weights, as they correspond to the Fourier transform of the bath correlation function $\langle {{\rm{e}}}^{{\rm{i}}{H}_{B}t}{B}^{a}{{\rm{e}}}^{-{\rm{i}}{H}_{B}t}{B}^{a}\rangle$, but for algorithmic purposes, it is natural to use universal transition weights. Also, we have studied only the dissipative part of the Davies generator (equation (5)). However, the full Davies generator also includes a Lamb-shift term that somewhat resembles our coherent term, but depends on additional details of the bath correlation functions. As the Lamb-shift term commutes with the Hamiltonian, adding this term keeps the Gibbs state stationary.

Recently, there have been several attempts to derive more realistic master equations for many-body systems by avoiding the λ → 0 limit needed for the rotating-wave approximation. This essentially translates to introducing Lindblad operators with a finite energy uncertainty, or, equivalently, in the time domain, integrals over Heisenberg-evolved jump operators weighted by a finite-width window function. In particular, the coarse-grained master equation⁵⁹ takes a very similar form as ours (equation (7)), but its operator Fourier transform features a uniform integral $\frac{1}{\sqrt{2{\rm{\pi }}T}}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}{{\rm{e}}}^{{\rm{i}}Ht}{A}^{a}{{\rm{e}}}^{-{\rm{i}}Ht}{{\rm{e}}}^{-{\rm{i}}\omega t}{\rm{d}}t$, and its coherent term also resembles ours (equation (19)) but does not seem to strictly enforce detailed balance. The universal Lindblad equation⁶⁰ combines (square root of) the transition weights into the operator Fourier transforms and is closer to subsequent constructions^61,62. Very recently⁶³ derived a master equation with exact KMS-detailed balance using slightly different operator Fourier transform weights compared with ours. All the above recent master equations feature some finite energy uncertainty, derived from various system–bath parameters, that parallels our tunable Gaussian width σ.