Thermal State Preparation via Rounding Promises

A promising avenue for the preparation of Gibbs states on a quantum computer is to simulate the physical thermalization process. The Davies generator describes the dynamics of an open quantum system that is in contact with a heat bath. Crucially, it does not require simulation of the heat bath itself, only the system we hope to thermalize. Using the state-of-the-art techniques for quantum simulation of the Lindblad equation, we devise a technique for the preparation of Gibbs states via thermalization as specified by the Davies generator. In doing so, we encounter a severe technical challenge: implementation of the Davies generator demands the ability to estimate the energy of the system unambiguously. That is, each energy of the system must be deterministically mapped to a unique estimate. Previous work showed that this is only possible if the system satisfies an unphysical 'rounding promise' assumption. We solve this problem by engineering a random ensemble of rounding promises that simultaneously solves three problems: First, each rounding promise admits preparation of a 'promised' thermal state via a Davies generator. Second, these Davies generators have a similar mixing time as the ideal Davies generator. Third, the average of these promised thermal states approximates the ideal thermal state.


Introduction
Motivation.Preparing Gibbs states is a major task for quantum computers.There are several reasons for this.First, the Gibbs state is one of the most important states of matter.For quantum models comprised of many locally interacting particles, it describes a wide range of physical situations, relevant to condensed matter physics, high energy physics, quantum chemistry [Alh22].Therefore, to use the quantum computer as a universal simulator of quantum systems, it is desirable to be able to prepare Gibbs states.Second, for general Hamiltonians, the Gibbs state is a crucial ingredient in some quantum algorithms such as those for solving semidefinite programs [VAGGdW20, BKL + 19, vAG19, GSLW19] and for training quantum Boltzman machines [KW17].Third, the problem of estimating the quantum partition function, which is connected to the problem of approximately preparing the Gibbs state, plays an important role in quantum complexity theory [BCGW21].

Background and overview of prior work.
There are three main approaches to preparing Gibbs states.
The first one is a Grover-based approach in which an initial state is mapped onto a certain purification of the Gibbs state at inverse temperature β (see [PW09,CS17]).The resulting running time is dictated by the overlap between the initial and target states, which is exponentially small in the system size.The Gibbs state can be approximately prepared in time on the order of D/Z β , where D denotes the dimension of the quantum system and Z β the partition function at inverse temperature β.In [HMS + 22] a quantum algorithm is presented for preparing a purification of the Gibbs state for the Hamiltonian H 1 = H 0 + V at inverse temperature β starting from a purification of the thermal state of H 0 .
The second approach is based on Davies generators, which is a differential equation that describes how nature thermalizes a quantum system to its thermal equilibrium [Dav76,Dav79].Davies generators are special cases of Lindbladians [Lin76,BP02], which describe the most general continuous-time Markovian dynamics of an open quantum system, i.e., a quantum system that is weakly coupled to the environment and the dynamics of the environment are fast enough so that the information only flows from the system to the environment while no information is flowing back to the system.The method in [CB21] simulates a Davies generator by attaching a heat bath and simulating time evolution on the joint system while repeatedly refreshing the bath.Their algorithm in some sense follows the original derivation of the Davies generator as the limit of such joint system-bath Hamiltonian time evolution.When the system Hamiltonian satisfies the Eigenstate Thermalization Hypothesis (ETH), they show that the implemented quantum map not only converges to the desired Gibbs state, but also does so in polynomial time.
The third approach is based on the quantum Metropolis algorithm [TOV + 11].This approach also avoids the exponential scaling when the system Hamiltonian satisfies ETH [SC22,CB21].
The advantages of the second approach are two fold.First, for every quantum system that thermalizes fast in nature, it is expected that our algorithms can also prepare its corresponding Gibbs state efficiently without suffering from the exponential dependence on the number of qubits.Second, this approach fits well for many physics-motivated applications.For example, we can use our algorithm to prepare a "partially thermalized" (in the natural thermalization process) thermal state, which might be of interest in some scenarios.
Main result.The present work examines how to approximately prepare Gibbs states for arbitrary system Hamiltonians by simulating time evolutions according to carefully engineered Lindbladians.These Lindbladians are derived from an ideal Davies generator having the Gibbs state as unique fixed point.
There are some similarities but also important differences to the work [CB21].Obviously, both are based on Davies generators.In contrast, we do not approximate Lindblad evolution with the help of the Hamiltonian evolution of the system and a bath, a method with which it is provably impossible to achieve linear scaling in evolution time (see [CW17]).Instead, we rely on a method for directly simulating Lindbladian time evolution specified by any jump operators.Using this method can lead to a reduction in resources: specifically, we achieve linear scaling in the mixing time t mix .In addition, a direct Lindblad simulation approach avoids the complication of dealing with the dynamics of the bath and the interaction Hamiltonians.We also seek to prove that our quantum map approximates the Gibbs state for any system Hamiltonian, that is, we do not need to make an assumption such as ETH.
The quantum algorithms for simulating Lindbladian time evolution in [CW17,LW22] assume that the jump operators have been suitably encoded.Unfortunately, it is not possible to construct jump operators of a Davies generator due to inherent imperfections of energy estimation of general Hamiltonians.However, we show how to construct a family of Lindbladians from the given Davies generators such that simulating them with the help of the simulation algorithms in [CW17,LW22] and taking the average of the resulting quantum states provides a good approximation of the Gibbs state.This is formulated in more detail in the theorem statement below.
Theorem 1 (Main result -informal statement).Assume we are given block encodings of the Hamiltonian H, coupling operators S α , and a filter function G.With appropriately chosen S α , these give rise to a Davies generator L that has the Gibbs state ρ β for inverse temperature β as a unique fixed point.Assume that after time t mix the time evolved state exp(t mix L)(σ 0 ) is ε-close to the Gibbs state ρ β for any initial state σ 0 .
We engineer a certain family of 2 r many Lindbladians L(j) from the above Davies generator L. These Lindbladians have a new 'attenuated' mixing time t mix , and their jump operators can be encoded efficiently with imperfect energy estimation.This makes it possible to simulate their time evolutions exp t L(j) .We prove that the average is ε-close to the Gibbs state ρ β for any initial state σ 0 .Furthermore, we show how to implement the time evolution according to L(j) to arbitrarily small failure probability δ L , and that the total number of invocations of the block encoding of the Hamiltonian is bounded by: where γ is an attenuation coefficient that affects the attenuated mixing time tmix .
Remark.After the first version of this manuscript was made public, recent work [CKBG23] resolved an open question raised in this manuscript (see Section 7) on utilizing the block encoding of the form j |j⟩ ⊗ L j .Using the Theorem III.1 of [CKBG23] together with slight adaptation of the simulation algorithm in [LW22], the complexity of our algorithm can be improved to Õ( tmix • γ −1 • βε −2 ).We also note that the additional factor of log 2 (β/ε) can be removed via existing techniques from [Ral21].
In Section 5 we give some numerical experiments indicating that the attenuated mixing time tmix is on the order of the original mixing time t mix for suitably chosen γ −1 .
We formulate in Section 7 some open research questions whose solution could lead to improvements of our current methods.
Technical overview.As mentioned above, the implementation of jump operators of a Davies generator requires perfect energy estimation.More precisely, what we mean by perfect is that energy estimation would have to be unambiguous: for each energy of the Hamiltonian, it must yield a unique energy estimate.Unfortunately, energy estimation unavoidably produces superpositions of different energy estimates for general Hamiltonians: if |ψ⟩ is an eigenstate of the Hamiltonian with eigenvalue λ, energy estimation implements a map: where λ1 and λ2 are two different estimates of λ.This superposition of estimates cause constant-size errors in quantum algorithms implementing Davies generators, and must be eliminated.With perfect estimation, the superposition on the second register is not present and there is always a unique estimate.It is possible to construct an 'approximate Davies generator' from an energy estimation algorithm that produces superpositions of estimates.However, the resulting dynamics no longer correspond to the Davies generator of any particular Hamiltonian, making it challenging to analyze.To our knowledge, no method exists for rigorously proving the accuracy of an algorithm based on such an approximate Davies generator.We highlight some of the challenges of this task in Appendix C.
It was shown in [Ral21] that perfect energy estimation is possible when the Hamiltonian satisfies a 'rounding promise' assumption.The rounding promise prohibits the Hamiltonian from having any eigenvalues that induce a superposition of estimates.However, such an assumption on the Hamiltonian is extremely unphysical and will not be satisfiable in practice.The main technical idea of the present work is to shift the notion of a rounding promise away from the Hamiltonian itself, but rather to a family of states.The basic idea is very simple: if a quantum state has no support on any eigenstates whose eigenvalues have a superposition of estimates, then it is as if the Hamiltonian did not have such eigenvalues.This family of states is defined by a 'promised subspace'.
Since perfect energy estimation is possible on a promised subspace, implementation of Davies generators is possible as well.While the analysis involves a wide variety of error parameters, we find that all of them admit a mathematically rigorous treatment.
Our construction relies on just one assumption: that we can construct coupling operators for each promised subspace that ensure that the Davies generator converges in a reasonable amount of time.We give numerical evidence that projecting a coupling operator into a promised subspace does not significantly reduce the Davies generator's convergence time.Notice how this assumption has nothing to do with the protocol's accuracy, only its convergence time.
We now summarize the different types of errors that occur in our implementation, and the rough ideas behind keeping them under control: • Lindblad simulation error: There exists a quantum algorithm, e.g.[CW17], for simulating Lindblad evolution given a description of the jump operators of a Lindbladian L.
For evolution time t and accuracy δ sim , this algorithm implements a quantum channel with accuracy δ sim in terms of the diamond norm to the ideal channel e Lt .In our case, we simulate the dynamics that approximates a special Lindbladian called the Davies generator.
• Knowledge of mixing time: Given a Hamiltonian H, coupling operators S α , and a filter function G, the Davies generator is a certain Lindbladian whose fixed point is the desired thermal state ρ β = e −βH / Tr e −βH .We assume we are given the Hamiltonian and suitable coupling operators of the Davies generator together with a bound t mix after which, the state is sufficiently close the desired thermal state.
• Precision and failure probability of energy estimation: The jump operators of a Davies generator depend on the coupling operators and on the Bohr frequencies (differences of the energies) of the Hamiltonian and the corresponding pairs of eigensubspaces.To realize these, we rely on energy estimation of the Hamiltonian.
By restricting to the promised subspace, we can guarantee using techniques from [Ral21] that for every energy λ there exists a single unique estimate λ.There remain two other kinds of error, both of which can be dealt with rigorously.
First, the energy estimation map is only implemented with a failure probability δ est that can be exponentially suppressed.By phrasing this error as a deviation in spectral norm of the Lindbladian jump operators, we can show that this error is not blown up when the Davies generator is simulated for a long period of time.
Second, the cost of energy estimation scales linearly with the precision of the estimate.
To analyze this, we leverage a trick from [PW09]: instead of preparing the thermal state for the original Hamiltonian, one can interpret the resulting process as preparing the thermal state for a rounded Hamiltonian.Since the norm of the difference of these Hamiltonians is small, the corresponding thermal states are close to each other.
• Preparation of an initial state on the promised subspace: We implement a Davies generator whose dynamics are trivial outside the promised subspace.In order to prepare a promised thermal state, we must feed the Davies generator with a state that is exclusively supported on the promised subspace.We achieve this by taking an arbitrary input state and measuring a two-outcome 'left-right' POVM.Depending on the POVM outcome, we know that either a 'left' or a 'right' rounding promise is satisfied by the post-measurement state.
• Approximation of the ideal Gibbs state: Our goal is to prepare a Gibbs state supported on the entire Hilbert space.But our protocol only prepares promised Gibbs states, which are only supported on the promised subspace.Thus, individual promised thermal states are generically far in trace distance from the ideal state.To deal with this, we show that there exists an ensemble of rounding promises such that the ensemble average of all the promised Gibbs states is an accurate approximation of the ideal thermal state.
The basic idea of this ensemble is that the probability of any particular energy being excluded can be made arbitrarily small.
• Leakage and attenuation errors in removal of rounding promise: The constructed coupling operators of the Davies generators for the well-rounded Hamiltonians on the promised subspaces have two important types of errors, namely, 'leakage' and 'attenuation.' The leakage error δ leak measures how much coupling operators and initial states 'leak' outside a promised space.Fortunately, δ leak can be made exponentially small by using quantum singular value transformation.
Blocks of coupling operators corresponding to some pairs of energies can be 'attenuated', i.e., multiplied by small positive numbers.However, as long as attenuation remains nonzero, the fixed-point remains the thermal state for the well-rounded Hamiltonian on the promised subspace.Unfortunately, the mixing speed can be negatively affected.
• Mixing assumption for projected/attenuated coupling operators: In principle, it could happen that the ideally projected coupling operators do not guarantee convergence to the thermal state of the well-rounded Hamiltonian on the promised subspace anymore.Moreover, even if they do, attenuation could increase the convergence time.However, we provide numerical evidence that these unfavorable situations do not typically occur.For the theoretical analysis of our quantum algorithm, we must assume that the attenuated coupling operators for the well-rounded Hamiltonian on the promised subspaces have similar mixing behavior as the original coupling operators.

Preliminaries
In this section, we first review some preliminaries about Gibbs states and Davies generators in order to establish notation and to rigorously define our goal: to prepare a Gibbs state by simulating the time evolution of a Davies generator.To do so, we leverage an algorithm for simulating general Lindblad evolution [LW22].
In order to implement the Davies generator, we require conditions under which the energy of a Hamiltonian can be estimated without producing superpositions of different energy estimates.So, in Section 2.2 we establish the notion of a rounding promise, and review a result from [Ral21] how a rounding promise can ensure that each energy is rounded to a unique energy when performing energy estimation.
Here are some notations and conventions used in this paper.We use Z + to denote the set of all positive integers.Let H be a Hilbert space.We use L(H) to denote the collection of all linear operators (matrices) of the form: A : H → H.For a matrix A, the spectral norm ∥A∥ is the largest singular value of A, and the trace norm ∥A∥ 1 is the sum of the singular values.For a superoperator Λ, the induced trace norm of Λ, denoted by ∥Λ∥ 1 , is defined as ∥Λ∥ 1 = max A̸ =0 ∥Λ(A)∥ 1 /∥A∥ 1 .If Λ is acting on L(H) for some Hilbert space H, then, the diamond-norm of Λ, denoted by ∥Λ∥ ⋄ , is defined as ∥Λ∥ ⋄ = ∥Λ ⊗ I∥ 1 , where I is the identity map acting on L(H).If A and B are matrices, then we say A ≥ B if A − B is positive semidefinite.Finally, a block encoding of A is a unitary matrix U A that, for some ancilla-count k, satisfies: (4)

Gibbs states and Davies generators
We now define some fundamental concepts required to state our quantum algorithm and to analyze its performance.
Definition 2. Let H be a Hamiltonian on the Hilbert space H with eigendecomposition Here, the Π i are projectors onto the subspace with energy λ i .We assume that the spectrum of H is contained in the interval [0, 1].For inverse temperature β > 0, the Gibbs state ρ β is the state such that The partition function Z β is the normalization factor given by Our quantum algorithm makes it possible to approximately prepare thermal states ρ β .It is based on the Davies generators defined below.The Davies generators describe dissipative dynamics that converge to thermal states.Definition 3. Let {S α } α be a collection of Hermitian operators acting on H.The Davies generator with respect to the system Hamiltonian H and the coupling operators S α is the Lindbladian L in the Schrödinger picture given by The jump operators S α (ω) are enumerated by the Bohr frequencies ω of H and are obtained from the coupling operators S α by The filter function G(ω) is a real-valued function satisfying G(ω) = e βω G(−ω).
The time evolution of quantum states is given by the quantum channels for t ≥ 0. We say that a quantum state ρ is a fixed point of the Davies generator L if It is known that if the fixed point of the Davies generator L is unique, then the Lindbladian time evolution is relaxing in the sense that for any initial state σ (see [Nig19] and the references therein).The converse clearly holds.
Definition 4. We call the coupling operators S α of a Davis generator L mixing if the thermal state ρ β is the unique fixed point of L. In this case, t mix denotes the mixing time of L, that is, the smallest time t such that T t (σ) is guaranteed to be sufficiently close to the desired thermal state ρ β for any initial state σ.
There are several sufficient conditions guaranteeing the uniqueness of the fixed point.For instance, it can be shown that the thermal state ρ β is the unique fixed point of the Davis generator L, if the matrix algebra generated by the coupling operators S α is the full matrix algebra (this statement follows from [Spo77]; see also the discussion in [Nig19] for an overview of other sufficient conditions).
The Lindbladian L of the Davies generator can also be written in terms of a collection of jump operators {L ω,α } ω,α : where L ω,α = G(ω)S α (ω).Having brought the Davies generator into this form, we can leverage existing results for Lindblad simulation [CW17].
In particular, say we can implement a block encoding of the operator O L that implements the jump operators as follows: Then, given access to O L , we have technical tools to simulate the time evolution e Lt .We use the simulation algorithm from [LW22], which simplifies the simulation algorithm of [CW17] and also generalizes their input models to block-encodings.In our adaptation of the theorem above from [LW22], we have replaced a quantity ||L be || with a bound O(m) which follows from the fact that ||L α,ω || ≤ 1.
So, the central task is to implement a block encoding of O L for the Davies generator, which essentially amounts to implementing block encodings of the S α (ω) operators.Once this is accomplished, we can simulate the evolution e t mix L on any input state, and obtain an approximation of the thermal state ρ β .

Rounding promises
Unfortunately, it is not possible to block encode the jump operators S α (ω) of Davies generators without prior knowledge of the spectrum of H. Ideally, we would like to implement a unitary that performs the isometry: that is, computes a binary representation of the energy λ i into a new register.However, there are two unavoidable limitations.First, λ i can only be estimated to precision ε using only O(1/ε) many resources.So, unless we select ε to be less than the smallest gap between any pair of energies, which requires at least O(dim(H)) many resources, we will be unable to distinguish certain energies.We find that, leveraging a trick from Poulin and Wocjan [PW09], this error can be dealt with formally.However, the second limitation is more challenging to deal with.Rall [Ral21] observed that for any energy estimation algorithm there will exist particular energies λ i such that a corresponding eigenstate |ψ i ⟩ will produce a superposition of estimates: where α, β are complex coefficients and λi , λ′ i are two estimates of the energy λ i .Essentially, certain energies λ i are indecisive about which direction they want to round, and end up being rounded up or down in superposition.
The superposition of rounding directions creates cross terms between the rounding results in the construction of approximate block encodings of the jump operators S α (ω), which cause errors in spectral norm of constant size no matter how high the precision of energy estimation, at least using SVT-techniques.A new approach is needed.
To remove the superposition of rounding errors, Rall [Ral21] introduced the notion of a 'rounding promise'.Observing that certain λ i are indecisive about their rounding, the rounding promise simply asserts that these λ i do not appear in H.This assumption on H is very unphysical.In our work, we introduce a new notion of a rounding promise that can be guaranteed without restricting the Hamiltonian, and makes a similar guarantee.The central idea is similar: certain energies λ i are disallowed.But instead of being an assumption on the Hamiltonian itself, our rounding promises define subspaces of the Hilbert space.
x+1 for all x.We always use the convention that the first interval starts at 0 and the last one ends at x , a is the projector onto the eigenspaces of x , b x ], that is: The promised subspace projector P (M ) is the projector onto the eigensubspaces of H whose eigenvalues lie in M , that is: The promised subspace P (M ) is Remark 7. When the rounding promise M is fixed, we often omit the superscript (M ) to abbreviate the notation.For instance, we write P instead of P (M ) for the promised subspace projector.
Rather than restricting the Hamiltonian itself, we have defined a subspace P (M ) on which energy estimation can be performed without superpositions of rounding errors.Furthermore, the rounding promise also conveniently specifies the estimates themselves: estimating λ amounts to identification of the index x such that λ ∈ [a x ].Indeed, we have the following result: Proposition 8 (Energy estimation given a rounding promise [Ral21]).Say M is a rounding promise, say κ is the length of the smallest gap.Suppose also that the number of intervals s (M ) satisfies s (M ) ≤ 2 n+1 , so each label can be thought of as a bit string x ∈ {0, 1} n+1 .
Then, for any δ est > 0, there exists a family of operators Px that approximate P x in the sense that: In fact, P x , Px and P (M ) commute.Moreover, say we have a block encoding of a Hamiltonian H.Then, for any δ est , there exists a quantum circuit that implements an isometry Ẽ(M) satisfying: where the G Px are isometries satisfying G † Px G Px = Px .This circuit can be implemented using O(n 2 κ −1 log δ −1 est ) applications of the block encoding of H and 1-and 2-qubit gates.
We prove this in Appendix B. The factor of O(n 2 ) can be removed with some additional care using techniques from [Ral21], but we keep it here to simplify the algorithm.This factor corresponds to the O(log 2 (β/ε)) in the performance in Theorem 1.

Algorithm overview
While we cannot use Davies generators to exactly prepare the Gibbs state ρ β on the entire space, we can prepare the promised Gibbs states ρ (M ) β which are restricted to a particular promised subspace P (M ) .Definition 9. Let M be a rounding promise.For each x, let m is the density matrix supported only on the promised subspace P (M ) such that The promised partition function is the normalization factor With this idea in place, we can give an informal high-level overview of our algorithm.Each of the major challenges in the algorithm's construction is treated in a section of this paper.
First, promised thermal states are in general not close to the ideal thermal state.However, we find that a suitable ensemble of promised thermal states is satisfactory.
Claim 10 (Section 4).The desired thermal state ρ β for the Hamiltonian H can be approximated by an average of promised thermal states ρ (M j ) β for suitably chosen rounding promises M j , where j ∈ {0, . . ., 2 r − 1}.In particular, let ρ * be the ensemble average over the promised thermal states of the M j .Then, if we perform energy estimation to n bits of precision: The main idea is that the rounding promises M j have to be chosen such that each eigenvalue λ i of the Hamiltonian H is contained in at least 2 r − 1 rounding promises M j .That way any individual eigenspace can be missing with probability at most 2 −r .
Second, we must restrict the dynamics of a Davies generator to a promised subspace.This 'promised Davies generator' will then be used to prepare the promised thermal states.

Definition 11 (Promised Davies generator). Say L is a Davies generator on the full
Hilbert space with coupling operators S α .Say M is a rounding promise, and we are given an 'attenuation operator' A (M ) that commutes with the Hamiltonian and satisfies: Then, the promised Davies generator L (M ) is a Lindblad operator defined by the jump operators L ν,α given by: Its Bohr frequencies ν are differences of the form m y , where m Finally, the coupling operators S (M ) α are given by In Section 4 we will construct the attenuation operator A (M ) as well as show that the promised Davies generator has the following properties: Claim 12 (Section 5).Say L is a Davies generator on the full Hilbert space with coupling operators S α .Then, for any rounding promise M , there exists a promised Davies generator L (M ) with the following properties: • If an input state σ is supported only on P (M ) , the output state e tL (M ) (σ) will be as well.
• The unique fixed point of L (M ) is the promised thermal state ρ (M ) β .
• Numerical simulations indicate that the mixing time tmix of L (M ) is not too much slower than that of L.
The purpose of the attenuation operator A (M ) is to ensure that the output of the Lindblad evolution according to L (M ) is confined to the promised subspace.Then, so long as time evolution starts in an initial state σ(M) that is approximately supported on the promised subspace P (M ) , then the output state will be as well.
The attenuation operator can make the mixing time of the Lindblad evolution slower.First, even in the ideal case when A (M ) = P (M ) , the elimination of certain energy eigenspaces may result in slower mixing.Second, it is not possible to project onto the promised subspace perfectly, for similar reasons that it is not possible to estimate energies without superpositions of rounding errors.We solve this problem by attenuating some of the eigenspaces in the promised subspace as well.This may also result in slower mixing.We show via numerical study that neither of these effects make the mixing time too much worse in practice.
Remark 13 (Well-rounded Hamiltonian).Note that the promised Davies generator L (M ) is block diagonal with respect to the orthogonal decomposition H = P (M ) ⊕ P (M ) ⊥ .It acts as 0 on P (M ) ⊥ , and acts on P (M ) as the Davies generator with respect to the promised Hamiltonian and coupling operators S (M ) α (where we view all operators as restricted to the promised subspace P (M ) ).
Third and finally, we show how to construct an approximate block encoding O L (M ) for each promised Davies generator.This lets us leverage Proposition 5 to simulate the Lindblad dynamics and prepare the promised thermal states.
Claim 14 (Section 6).Say we are given block encodings of some coupling operators S α , and one of the rounding promises M ∈ {M 0 , . . ., M 2 r −1 }.Suppose we perform energy estimation to n bits of precision.Then, there exists a quantum algorithm that implements time evolution for time t according to L (M ) to any precision δ L > 0 in the diamond norm using: invocations to the block encoding of H, where γ is an attenuation coefficient introduced in Section 5. We accomplish this by constructing and simulating time evolution according to an approximate promised Davies generator L(M) .
4 Averaging together promised thermal states By selecting a rounding promise M restricting to a promised subspace P (M ) , we have gained the ability to approximately prepare promised Gibbs states ρ β (M ) .These Gibbs states have support only on P (M ) , and thus may be far from the true thermal state ρ β .In this section we show how to construct several rounding promises M j ∈ {M 0 , . . ., M 2 r −1 }, such that the ensemble average over the ρ β (M j ) 's is provably close to ρ β .Furthermore, to prepare an approximation of ρ β (M j ) , we require an initial state σ that is also approximately supported only on P (M j ) .This state is then fed as input to the promised Davies generator, whose dynamics are trivial outside of P (M j ) .We also show how to prepare these initial states in this section: it is achieved by measuring a POVM called the 'left-right POVM' on an arbitrary initial state.
The left-right POVM splits the computation into two branches, each corresponding to a family of left rounding promises L, L j and right rounding promises R, R j .We use M as a symbol to denote either L or R depending on which branch we are in.An overview of the algorithm as a whole is as follows: Step 1.Take an arbitrary initial state σ, and measure the left-right POVM defined in Section 4.1.Depending on the measurement outcome, the resulting state σ(M) will be approximately supported on P ( M ) , where M ∈ { L, R} is one of two 'fine-grained' rounding promises.
Step 2. Pick an index j ∈ {0, ..., 2 r − 1} uniformly at random.This index determines a 'coarse-grained' rounding promise M j , which is either L j or R j depending on the measurement outcome of the left-right POVM.These are defined in Section 4.2.Since the M j are coarse-grainings of M , the input state σ(M) is supported only on P (M j ) for any j.
Step 3. Use the promised Davies generator to approximately prepare ρ β (M j ) .This is discussed in Sections 5 and 6.

Analysis.
Let ρ * M be the ensemble average over the index j that we selected in step 2 and used to prepare ρ β (M j ) .We show that for both M ∈ {L, R}, the density matrix ρ * M is close in trace distance to the ideal thermal state ρ β .We perform this analysis in Section 4.3.This protocol is represented diagrammatically in Fig. 1.

Initial State
Step 1: Step 3: Apply Promised Davies Generator Figure 1: Sketch of a protocol that leverages promised Davies generators to prepare an approximation ρ * M of an ideal thermal state ρ β (for M ∈ {L, R}).In the figure above, ρ ∈ M is a shorthand for 'ρ is approximately supported entirely on P ( M ) ', which we define more rigorously in Definition 15.

The Left-Right POVM
The Left-Right POVM has the purpose of producing a quantum state σ(M) which is guaranteed to be approximately supported on one of the two fine-grained rounding promises M ∈ { L, R} depending on the measurement outcome.In this section we show how to implement this POVM, and prove that the post-measurement state has the desired property.We define this notion rigorously now: Definition 15.Say M is a rounding promise and σ a density matrix.We say σ is δ supapproximately supported on M if In other words, if we measure the projector P (M ) in order to project into the promised subspace P (M ) , we succeed with probability ≥ 1 − δ sup .
The construction depends on two parameters: n and r.Here n is the number of bits of precision for energy estimation, and r determines the number of coarse-grained rounding promises, which is 2 r .The implementation of the left-right POVM as well as energy measurement using Proposition 8 will require O(2 n+r ) applications of the block encoding of H.This is because all the rounding promises in this construction have a minimum gap size κ = 2 −n−r−2 .
The quantities 2 −n and 2 −r correspond to two different errors on the final state.The quantity 2 −n corresponds to the accuracy of energy estimation, and 2 −r is the probability with which any particular energy is excluded from the ensemble.We will show in Theorem 24 that if output state is ρ * , then we have the guarantee is In our final construction in Section 6 we will select 2 n ∼ β/ε 2 and 2 r ∼ 1/ε, achieving an accuracy of ∼ ε.
As we present the construction, we recommend following along with Fig. 2. The main idea is that we would like to eliminate small regions of the spectrum via the Left-Right POVM, which is defined by an operator P LR whose spectrum is sketched in Fig. 2. When P LR has no support on an eigenstate, and we observe the POVM outcome corresponding to P LR , then the output state has no support on that eigenstate.Consequently the output state satisfies the rounding promise L. Similarly, the rounding promise R is defined by the eigenstates on which P LR has eigenvalue 1.The key property of L and R that we require for the remainder of the construction is that there are 2 n+r many evenly spaced gaps in the spectrum.Later, we will define coarse-grained rounding promises that close all but 2 n of these gaps at random, so that the probability of any individual energy being excluded by the coarse-grained promise is at most 2 −r .
Definition 16 (The fine-grained rounding promises).For n, r ∈ Z + , let the rounding promises L and R be obtained by taking the full interval [0, 1] and removing some subintervals (a ′ , b ′ ).Specifically: Let the [a x , b x ] denote the resulting connected components.These closed intervals define the rounding promise in L and R respectively.
To prepare a state supported only on L or R, we need to take our input state σ and remove all the support on the eigenspaces in the deleted regions above.To this end, we construct a block encoding of an operator P LR that makes P LR σP LR be supported only on L and (I − P LR )σ(I − P LR ) be supported only on R.
Lemma 17 (Left-right projection operator).For any n, r ∈ Z + and any δ sup > 0, there exists a Hermitian operator P LR that satisfies P 2 LR Π i ≤ δ sup /3 for λ i ̸ ∈ L, and satisfies This operator has a block encoding that can be implemented using O((n + r)2 n+r log δ −1 sup ) invocations of a block encoding of H.
Proof.See Appendix B. We will ensure P LR commutes with the Hamiltonian, and takes the form P LR = i p LR (λ i )Π i for a function p LR that is plotted in Fig. 2.
It remains to show how to turn P LR into a POVM, and how to implement that POVM using a quantum circuit.The starting point for the implementation is the block-measurement theorem from [Ral21].This theorem takes a block encoding of an operator A satisfying |A 2 − Π| ≤ δ for some projector Π, and uses it to implement a quantum channel 4 √ 2δ-close in ⋄-norm to the isometry: This does not quite suffice for us, since P LR is not a projector.But, if we allow ourselves to introduce a constant amount of postselection, then we can implement a very similar operation that suffices for our purposes.
Lemma 18 (Postselective block-measurement).Say A is a Hermitian matrix with a block encoding U A .Then there exists a quantum circuit with postselection using one U A and one , where: The postselection succeeds with probability at least 1/2.

Proof. Say the block encoding U
We consider the following circuit: where the CNOT above denotes the operator When the postselection succeeds, this implements V A as desired.It remains to bound the postselection probability, which is: = Tr(σ Since A is Hermitian and If we measure the top qubit after applying V A , we will implement a POVM defined by the operator A. All the tools are in place to define and implement the left-right POVM and demonstrate that it guarantees that the post-measurement state σ(M) has the desired property of being approximately supported on M .One caveat is that the procedure only works if the probability of the observed outcome is bounded away from 0 by a constant.This is because the normalization of the state that occurs post-measurement may blow up the support outside the promised subspace.Fortunately, such a lower bound on the outcome probability is easily attained as we will see later.
Proposition 19 (Left-right POVM).For any n, r ∈ Z + and any δ sup > 0, there exists a two-outcome POVM with the following property.We label the two outcomes 'L' and 'R', and for any input state σ let the output state be called σ(M) for M ∈ {L, R}.Say p M is the probability of the M outcome, and suppose we obtain an outcome with p M ≥ 1/3.Then the output state σ(M) is δ sup -approximately supported on M .
Furthermore, there exists a quantum circuit that implements the POVM with success probability ≥ 1/2 using O(2 n+r ) invocations of a block encoding of H.
Proof.The circuit is just an invocation of postselective block-measurement from Lemma 18 with the left-right projection operator P LR from Lemma 17.After successfully applying the isometry V P LR we measure the label qubit and output 'L' in the 0 case and 'R' in the 1 case.This immediately gives the bound on the circuit complexity and success probability.It remains to show that the output state has the desired property.
The probabilities of the measurement outcomes are: The corresponding output states are: Another way of writing the results of Lemma 17 is: Having set the scene, we are ready the compute the probability of the output register being approximately supported on P ( M ) .We use the inequality Tr(XY ) ≤ Tr(X)Tr(Y ) that holds for all positive semidefinite matrices X and Y .We obtain: ≤ Tr(σ) • (δ sup /3) 2 ≤ δ sup /3 and (45) = Tr(σ So for M ∈ {L, R}, we have: Since we assumed p M ≥ 1/3, we have Tr σ(M) P ( M ) ≤ δ sup as desired.
Achieving p M ≥ 1/3 is very easy: we just repeat the protocol a couple of times and pick the outcome we saw the most frequently.
Corollary 20 (Selecting a high-probability outcome).For any δ fail , δ sup > 0 there exists a procedure that produces a quantum state σ(M) that, for some random M ∈ { L, R}, is δ supapproximately supported on M with failure probability at most δ fail .Here, 'failure' refers to selecting an outcome M with p M < 1/3.This procedure requires at most 20 • log δ −1 fail many re-preparations of the input state σ and implementations of the left-right POVM from Proposition 19.
Proof.The procedure is as follows.Let N be the smallest odd number greater than 20 log(δ fail ).We repeat the left-right POVM from Proposition 19 a total of N times, preparing a new initial state σ each time.We let M be the most frequently observed outcome, and return the output state of any of our attempts that measured the M result.
We fail if p M < 1/3.Let K be the number of times we observed the M outcome.Then, using the Chernoff-Hoeffding theorem, the probability that we observe K/N > 1/2 despite the fact that p M < 1/3 is: (50)

Coarse-grained rounding promises
After having obtained a state σ(M) that (on average) is approximately supported on a finegrained rounding promise M , the next step of the algorithm is to select a coarse-grained rounding promise M j ∈ {M 0 , ..., M 2 r −1 }.Because M ⊂ M j for all M j , we automatically have that σ(M) is approximately supported on each of the M j .
The purpose of the coarse-grained rounding promises M j is to ensure that the ensemble average ρ * M of the promised thermal states ρ β (M j ) is close to the ideal thermal state ρ β .
This would not be true if we just prepared the two states ρ β ( M ) instead and considered their ensemble average, as we will explain now.
Recall that the main mechanism of the rounding promise is to eliminate certain energy eigenspaces.On the one hand, this enables energy estimation without superpositions of rounding errors.On the other hand, since certain energies are eliminated, an individual promised thermal state could never guaranteed be close in trace norm to ρ β .
In an ensemble average ρ * M over promised thermal states ρ β (M j ) however, we can guarantee that the probability of each individual energy being eliminated is small (≤ 2 −r ).We achieve this condition by selecting each of the M j with probability 1/2 r , and showing each energy is eliminated by at most one of the M j 's.We will show in the next section that this suffices to guarantee that ∥ρ

An ensemble over ρ ( M ) β
does not have this property: if an energy λ i is eliminated by, say, L, then the probability of the eigenspace being eliminated in the ensemble is p L .If we leave the initial state arbitrary, then p L could be as large as 1.But even if we set σ to the maximally mixed state, thereby ensuring p L = p R = 1/2, the probability of eliminating the energy is still constant, resulting in a constant-size error.To guarantee an error of at most ε, we require O(ε −1 ) many rounding promises.We will show that they have width at most 2 −n , so the error from energy rounding will be at most ∼ 2 −n .
We require one more property of the M j 's: recall from Definition 9 that the promised thermal state's eigenvalues are ∝ exp −βm , where m x is the midpoint of the interval [a x , b x ].This is because when we perform energy estimation according to Proposition 8, all of the energies in the interval [a x , b x ] are rounded to m x .To avoid introducing too much error here, we must ensure that the intervals are not too wide.
The construction of the M j is also visualized in Fig. 2. The intuition for the construction is as follows.We want to achieve M ⊂ M j , so we construct the M j my merging gaps between the intervals of M .As discussed above, the widths of the intervals in M are ∼ 2 −n−r , and we want the widths of the intervals of M j to be ∼ 2 −n .This can be achieved by merging all but every 2 r 'th gap.The promise M j is defined by keeping every gap with index = j modulo 2 r .
Definition 21 (The coarse-grained rounding promises).For n, r ∈ Z + and M ∈ {L, R}, let the rounding promises M j ∈ {M 0 , ..., M 2 r −1 } be constructed by taking the fine-grained rounding promise M from Definition 16 and merging some of the gaps.
By merging the x'th gap, we mean the following transformation on a rounding promise: if the promise contains the intervals Then, M j is obtained by starting with M and merging all the x'th gaps where x ̸ = j mod 2 n+r .
Lemma 22 (Properties of the coarse-grained rounding promises).First, if σ(M) is δ sup -approximately supported on M , then it is also approximately supported on all the M j 's.
Second, say λ ∈ [0, 1] is any energy.Then at most one of the M j 's does not contain it.Here is another way of saying this: let 1 λ∈M j be 1 if λ ∈ M j and 0 otherwise.Then: (51) Proof.For the first claim, since the gap merging procedure ensures that M ⊂ M j , we have that P ( M ) ≤ P (M j ) .Consequently, following Definition 15: For the second claim, we first consider what would happen if we merged all the gaps in M : since the gap-merging procedure includes cases for the first and last interval, this would just result in the promise {[0, 1]}.So that means that for each M j the only energies that are missing are those present in the gaps (b x , a x+1 ) where x = j mod 2 n+r .For every gap index x, there is a unique j for which this equation holds.
For every λ there are two cases.If we have λ ∈ M , then since M ⊂ M j we also have λ ∈ M j for all M j .Otherwise, since λ ̸ ∈ M , λ must be contained in some gap of M .By the above, that gap is contained in all but one of the M j , so λ must be contained in all but one of the M j 's.
Finally, we must show that the intervals in the M j are not too wide.The intervals in M all have width 3/4 • 2 −n−r , and are 1/4 • 2 −n−r apart.If we close all but every 2 r 'th gap, then the intervals will have width at most: (53)

Ensemble analysis
In this subsection we prove that the ensembles ρ * M are close to the ideal thermal state ρ β .The error comes from two sources: from the accuracy of energy estimation (via n), and from number of rounding promises (via r).We first deal with the error from the accuracy of energy estimation by defining a new notion of an 'exact' promised thermal state ρ(M j ) β .
Definition 23.Let M be a rounding promise.Let the exact promised partition function Ẑ(M) and exact promised thermal state ρ(M) β be defined similarly as their non-exact versions, just leaving λ i intact rather than swapping them out with the midpoints m x of the intervals containing them.

Ẑ(M)
ρ(M) These exact promised thermal states are still restricted to the promised subspace, but unlike the regular promised thermal states ρ β (M j ) their energies are unrounded and left intact.
Thus, the difference between ρ(M j ) β and ρ β (M j ) captures the error from energy estimation only.
We will show that ρ(M j ) Second, we bound the error stemming from the elimination of certain energies in the rounding promises.In the previous subsection, we established that the probability that any particular energy λ is missing from M j is at most 1/2 r .Through careful consideration of the promised partition functions Z β (M j ) , exact promised partition functions Ẑ(M j ) β and the ideal partition function Z β , we use this property to show that ρ * M − ρ β ≤ 2 • 2 −r .Combining these two facts yields the main theorem of this subsection: Theorem 24 (Accuracy of the final ensemble).Consider the following ensembles of density matrices: Then, we have Proof.Combining Lemmas 27 and 28 with the triangle inequality, we obtain We start with the analysis of the rounding errors via the exact promised thermal states.To prove that these are close together, we leverage a fact from is the fidelity: This tool is convenient because a naive analysis of the trace distance yields an contribution from every eigenvalue, which accumulates into a total error proportional to the dimension.By considering the spectral distance in the Hamiltonian instead, we avoid the exponential blowup in the error.However, this trick is a bit tricky in our case since, when viewed as operators on the full Hilbert space, ρ β (M j ) is not full rank.But states of the form e −βH /Tr(e −βH ) are always full rank whenever H is.
To remedy this issue, we recall an observation made in Remark 13, where we noticed that if we view ρ β (M j ) as a density matrix over P (M j ) only, then it can be written as ) where H (M j ) is the promised Hamiltonian.We will leverage this viewpoint where we restrict to P (M j ) for Definition 26 and Lemma 25.
Definition 26.Let M be a rounding promise.Let the exact promised Hamiltonian Ĥ(M) be the part of H with support on the promised subspace P (M ) .But unlike the non-exact promised Hamiltonian H (M ) , the eigenvalues are unchanged.It is given by: Here, we view this operator as Ĥ(M) ∈ L(P (M ) ).
Lemma 27 (Rounding errors from energy estimation).For M ∈ {L, R} consider the following ensemble over exact promised thermal states: Proof.Recall from Lemma 22 that each of the intervals [a Following Lemma 25, we have F (ρ ) ≥ e −β2 −n , so, exploiting the fact that Finally, we average this error over the two ensembles: We move on to the error stemming from the elimination of energies in the promises M j .
Lemma 28 (The average of exact promised thermal states approximates the thermal state).For M ∈ {L, R}, we have Proof.Letting d i := Tr(Π i ) be the occupation numbers, we have: Here, for each λ i , we split the rounding promises M j into two categories: those not containing λ i , for which ρ(M) β has eigenvalue 0, and those containing λ i , for which ρ(M) After applying the triangle inequality to these two groups, these correspond to the first and second terms in the above sum, respectively.We will show that both of these terms are bounded by 2 −r , so the overall bound is 2 • 2 −r as desired.
The first term is bounded by: So all that is left is the second term.Observe that Ẑ(M j ) This allows us to remove the absolute value: Now we can bound the two terms corresponding to the difference 1/ Ẑ(M j ) The term with −1/Z β can be bounded by leveraging that for any λ i , at most one M (j) doesn't contain it, as shown in Lemma 22.That is, 2 −r j 1 So a bound on the second term is 1 − (1 − 2 −r ) = 2 −r .So, each of the two terms is at most 2 −r , so the total error is at most 2 • 2 −r .

Lindblad dynamics on the promised subspace
In the previous section we required the following capability: to take a state σ (M ) that is supported on a rounding promise M , and to transform it to the promised thermal state ρ β (M ) .Once we have this capability, we can prepare good approximations to the true thermal state ρ β .
We attain this capability by running a promised Davies generator L (M ) on the system, which we defined in Definition 11.We constructed L (M ) by starting with a regular Davies generator L and truncating it so that its dynamics are restricted to P. The basic idea is to consider the coupling operators S α of the original Davies generator, and to conjugate them with an attenuation operator A (M ) to obtain: Then, L (M ) is defined the exact same way as L, just with S α (M ) instead of S α .The other difference is that L is defined with respect to the ideal Hamiltonian H = i λ i Π i , and L (M ) with respect to the promised Hamiltonian In this section we detail the construction of the attenuation operator A (M ) and how that the resulting promised Davies generator has some desirable properties.Ideally, we would like to have A (M ) = P (M ) , but due to limitations of singular value transformation, we cannot achieve this.This results in two challenges.
The first challenge is leakage: we cannot perfectly eliminate the transfer of probability mass from P to P ⊥ ; we can only suppress it by a factor δ leak .Fortunately, this error is both exponentially suppressible, and its ramifications on the final output error can be quantified formally.We give a construction of an approximate attenuation operator Ã(M) that has support at most δ leak on P ⊥ .The operator A (M ) itself has no support on P ⊥ at all, ensuring that L (M ) does not leak.The leakage error can be dealt with by studying its approximate implementation L(M) which we construct in Section 6.
The second challenge, which is more significant, is attenuation.Since A (M ) is implemented via singular value transformation, it takes the form A (M ) = i a (M ) (λ i )Π i for some continuous function a (M ) : [0, 1] → [0, 1].We have already achieved that a (M ) (λ) = 0 for λ ̸ ∈ M .Ideally, we would like to also have a (M ) (λ) = 1 for λ ∈ M , but this means that a (M ) has a discontinuity at the boundaries of M .The function a (M ) must be continuous, so we need to smoothly interpolate from 0 to 1 near the edges of the intervals of M .Formally, we can define a truncated rounding promise M T whose intervals are a little bit narrower than those of M , and then ensure a (M ) (λ) = 1 for λ ∈ M T .However, this choice means that the operator also attenuates certain eigenvectors of H whose eigenvalues are close to the edges of the intervals of M .This attenuation might slow the rate of convergence of L (M ) .
Having highlighted the trade-offs of the construction, we formally define the attenuation operator A (M ) and its approximation Ã(M) .
Lemma 29 (Attenuation operators).Consider any rounding promise M with minimum gap width κ whose intervals are also at least κ wide.Then, for any leakage error δ leak > 0, and attenuation factor γ there exists an attenuation operator A (M ) and an approximate attenuation operator Ã(M) .They both commute with the Hamiltonian, and take the form: for some functions a (M ) , ã(M) : Consequently, we have ∥A (M ) − Ã(M) ∥ ≤ δ leak , and ∥A (M ) (I − P (M ) )∥ = 0 as required by Definition 11.Furthermore, there exists a block encoding of The proof is in Appendix B. We will show how to deal with leakage errors in Section 6.This section is dedicated to assessing the impact of truncation and attenuation on the dynamics of L (M ) .Our goal is to show that L (M ) is mixing (recall Definition 4) and that the mixing time t mix is not too much slower than that of L.
On a positive note, we emphasize that truncation to P as well as attenuation cannot cause L (M ) to fail to converge onto the promised thermal state unless the coupling operators S α are extremely contrived.For example, one way L could be mixing while L (M ) could fail to be so is if the S α encode a particular random walk over the energy eigenspaces.Specifically, if the random walk is over a bipartite graph where one half of the eigenspaces are in P and the other half are in P ⊥ , then a Davies generator based on S α may be mixing, while S (M ) α will not be.But construction of such a coupling operator requires careful knowledge over the energy eigenspaces.The chances of accidentally selecting such an S α in practice are slim.
The main problem is that the projection of S α to S α (M ) will slow the rate of mixing, and that attenuation only makes this problem worse.We assess the severity of this problem by performing numerical simulations on a particular physical system: the transverse field Ising model.Let (n 1 , n 2 ) be the width and height of 2d-grid of qubits, and let v be the strength of the transverse field.Let ⃗ i be indices in an n 1 × n 2 grid, let X ⃗ i and Z ⃗ i denote Pauli X and Z on the qubit at grid position ⃗ i, and let ⃗ i ∼ ⃗ j denote that the two grid positions are adjacent.Then, the Hamiltonian we consider in this section is: We study the dynamics of Davies generators defined with respect to this Hamiltonian, as well as the coupling operators X ⃗ i and Z ⃗ i for all grid positions ⃗ i.To assess the mixing time t mix we compute the spectral gap ∆ of the Lindbladian.
In Fig. 3 we plot spectral gap as a function of the attenuation factor γ for various parameters.In all our experiments, we verified that the ideal Davies generator L had the ideal thermal state ρ β as its unique fixed point, and that the promised Davies generators L (M j ) had the promised thermal state ρ (M j ) β as their unique fixed point.We analyze their exact versions with A (M ) rather than Ã(M) and P x rather than Px .
We find that in the ideal case of γ = 0, the ideal Davies generator and the promised Davies generators have about the same convergence rate -the promised Davies generators are a little slower.So while there is some slowdown in mixing time from projection, the slowdown is largely due to attenuation.We find that setting reasonably small values of γ achieves convergence rates similar to the γ = 0 case.Furthermore, this slowdown only appears to be significant for certain rounding promises M j .One interpretation of our technique in this paper is that we trade mixing time for a provable accuracy guarantee on the final thermal state.So, some amount of slowdown is acceptable in exchange for the rigor.80).We set β = 10.0 throughout.We see that the convergence rates of the promised Davies generators are a little slower than the ideal Davies generator as γ → 0. While the convergence rate slows down significantly with large γ, selecting reasonably small values of γ achieves a convergence time just as fast as the unphysical case of γ = 0.There is also significant variation among the promises M j .The software that produced these data is available at: https://github.com/qiskit-community/promised-davies-generator.

Implementing Lindblad dynamics
The goal of this section is to show how to simulate the evolution of the promised Davies generator for time t given a block encoding of its jump operators.This block encoding can be implemented by energy estimation.This construction also involves the coupling operator projected onto the promised subspace.Both energy estimation and the projection are not perfect, and they lead to errors in the block-encoded jump operators.We begin with a lemma that shows that a small perturbation of the jump operator will not change the Lindbladian by much.Then we give the detailed construction and prove the main theorem.
For simplicity, in this section we assume that there is only one coupling operator so that we can drop the index α in for the jump operators L ν,α and the coupling operators S α , and simply write L ν for the jump operators of a promised Davies generator (Definition 11), and write S for the coupling operator.Our analysis can be easily generalized to the case of more coupling operators: the number of jump operators just grows by the factor of the number of the coupling operators.
Furthermore, it is reasonable to suppose that ∥S∥ ≤ 1.First, this is generally the case in practice: the coupling operator S is usually selected to be a unitary transformation that scrambles the energy eigenbasis.Even if S were a block encoded operator, it would require a unitary implementation, so the only way to achieve ∥S∥ ≥ 1 would be via some virtual scale factor.Second, in the analysis of Hamiltonian simulation e iHt , it is usually assumed that ∥H∥ ≤ 1, since we can always absorb a rescaling of the Hamiltonian into the time t.The exact same argument applies for the simulation of Davies generators.
To begin with, we show that the distance between Lindbladians in terms of the diamond norm can be bounded by the distances of their jump operators in terms of the spectral norm.
Lemma 30 (Approximating Lindblad evolution via jump operators).Say L ν ∈ L(H (M ) ) are jump operators defining a Lindbladian L, and say there are m many of them.Say some Lν ∈ L(H) similarly define a Lindbladian L, and say these satisfy: Say we also have ∥L ν ∥ ≤ 1.Then, for any t > 0: Proof.First of all, by Eq. ( 81) and the triangle inequality, we have Now, recall that Let ρ be the matrix that achieves the induced trace norm of ∥L − L∥ 1 , i.e., Now, we have all the tools to prove the main theorem and present our quantum algorithm for simulating the Davies generator.

Theorem 32 (Implementation of the Davies generator given a rounding promise).
For some rounding promise M , say L (M ) is an promised Davies generator as defined in Definition 11.Then, for any δ L , t > 0, there exists a quantum algorithm implementing a channel δ L -close in diamond norm to e tL (M ) .If κ is the minimum gap of M , M has s (M ) many intervals, and γ is the desired attenuation factor for the S(M) , then this algorithm uses queries to the block encoding of S, queries to the block encoding of H, and additional 1-and 2-qubit gates.
Proof.We aim to invoke Proposition 5, so we need an oracle O L. This lets us approximately apply the map e t L. We then show that the Lν defining L are close to the L ν defining L (M ) , allowing us to invoke Lemma 30.
Recall the approximate energy estimation isometry Ẽ(M) from Proposition 8, which extracts estimates according to Px satisfying We define the block-encoded operator: where the controlled ± refers to adding/subtracting m x (M ) to/from the value of the target register (recall Remark 13, the control register has the value of x), and the kets on the right side of the circuit denote postselection onto that state.Next, obtain a block encoding of: Finally, let O L be given by: Now, using Lemma 31, we obtain x,y m Py − x,y m Combining Eqs. ( 113) and ( 116), we have So we have the desired property with (s (M ) ) 2 (δ leak + 2δ est ).Letting δ L := (s (M ) ) 2 (δ leak + 2δ est ) and invoking Lemma 30 with m = (s (M ) ) 2 , we have where the first inequality follows from the observation that for all integers k ≥ 0, Here the first inequality follows from the subadditivity of the diamond norm [Wat18, Propositon 3.48], and the last inequality is due to Taylor expansion.We assume S is given as a block encoding with normalizing constant 1 since ∥S∥ = 1.To use Proposition 5, it suffices to set , and to δ L -approximately simulate e L (M ) t by simulating e L(M) t .Observe that ∥L ν ∥ ≤ 1.Let k be the number of system qubits.The simulation algorithm costs additional 1-and 2-qubit gates.Note that the number of queries to O L is also the number of innovations to S.
To implement O L, we invoke Proposition 8 with precision parameter δ est and Lemma 29 with precision parameter δ leak .Each application of O L has the following number of invocations of U H : To achieve the final number of queries to the block encoding of H, we observe that the rounding promises M j have κ −1 ∈ O(2 n+r ) and s (M ) ∈ O(2 n ).Plugging these in, and accounting for the cost of the left-right POVM from Lemma 17 we obtain: This establishes Claim 14.Now we recall Theorem 24 which states that the final output state ρ * M satisfies: To achieve ∥ρ * M − ρ β ∥ 1 ≤ ε, we select n = log 2 (β(ε/2) −2 ) and r = log 2 (4/ε).Plugging these into the query complexity, we get: This establishes Theorem 1.

Some open questions
Our result achieved a time complexity that scales linearly in the mixing time t mix .However, the performance with respect to the inverse temperature β and accuracy ε is Õ(β 3 ε −7 ) which has plenty of room for improvement.One potential path that may yield an accuracy of Õ(βε −2 ) is to remove the linear dependence on the number of jump operators in the Lindblad simulation algorithm from Proposition 5. Proposition 5 demands block encodings of the individual jump operators L j , but the oracle O L we prepare may actually be much more powerful than this.To see this, first note that if we are given access to the isometry j |j⟩⊗K j , then implementing the channel is trivial as we already have the Stinespring dilation.In our algorithm, we implemented an oracle in the form of j |j⟩ ⊗ L j , which is close to the Stinespring dilation we want because for the infinitesimal approximation channel in [CW17], all but one Kraus operators are proportional to L j .That one special Kraus operator involves all the L j 's.Does there exist any special treatment of this special Kraus operator so that we can leverage the special structure of the oracle j |j⟩⊗L j to get rid of the O(m) dependence? 1 A central goal in our work is to attain a rigorous bound on the accuracy of the final output state.To this end, we assume that we are given a lower bound on the mixing time t mix , so that we know for how long the dynamics of the Davies generator must be simulated to achieve a high-accuracy output state.But this assumption is rarely the case in practice.Furthermore, the attenuation discussed in Section 5 may slow the mixing time and exasperate this problem.Is there a technique that can detect if the Davies generator has been run for long enough?One approach might be to purify the dynamics of the Lindblad simulation, and then use amplitude estimation to compare the resulting pure output state to an ideal thermal state purification.A reflection operator around a purification of the thermal state could be obtained via the techniques from [WT23].
The only piece of our method that eludes rigorous mathematical treatment is the impact of attenuation on mixing time tmix .Certainly the dependence of γ −1 in the circuit complexity the synthesis attenuation operators A (M ) is optimal, due to lower bounds on approximation of threshold functions with polynomials.But perhaps given additional knowledge about the Hamiltonian, there may be other approaches for selecting coupling operators S α that do not leak out of a given promised subspace.
∀λ ∈ [−2, 2] we have 0 ≤ Θ κ,δ (λ) ≤ 1 and: Given a polynomial approximation of a single transition from f (λ) = 0 to f (λ) = 1, we can approximate an arbitrary sequence of transitions by shifting and adding multiple such polynomials together.Above, we guaranteed that the approximation of the step function holds for λ ∈ [−2, 2], so that we can shift the polynomial by up to 1 and still have good approximation on the interval [−1, 1].
Lemma 34 (Construction of approximate projection polynomials).Say we have a collection of intervals {[a x , b x ]} in [0, 1] with a x < b x and b x < a x+1 , and each interval is labeled with a bit c x ∈ {0, 1}.Say furthermore that each of the intervals is at least κ far apart, that is, b x + κ < a x+1 .
Proof.We will construct P (λ) by adding together several Θ κx,δ/K for various values of κ x , all of which satisfy κ x ≥ κ.Since adding several polynomials doesn't change the degree, the resulting polynomial P (λ) has degree O(κ −1 log Kδ −1 ).
There will be one Θ κx,δ/K for each flip, that is, each pair of intervals with c x ̸ = c x+1 .Let t x := (b x + a x+1 )/2 be the midpoint between two intervals, and let κ x = a x+1 − b x ≥ κ be the distance between them.We construct: It remains to prove that P (λ) satisfies the desired property.By Lemma 33 and the choice of t x , κ x it is guaranteed that Θ κx,δ/K (λ − t x ), is always either ≤ δ/K or ≥ 1 − δ/K outside of the region [b x , a i+x ].So, for any interval [a x , b x ] we have that any of the Θ κx,δ/K (λ − t x ) is within δ/K of 0 or 1.Thus, for the purposes of analyzing P (λ), let us pretend for the rest of the proof that they are exactly 0 or 1.In doing so, we will be wrong by at most K • δ/K = δ.
Let us consider any particular interval [a x , b x ].We want to argue that for any λ in this interval, P (λ) = c i .We achieve this through induction in K.The base case is very easy: if K = 0, then all the c x are equal, so we can just set P (λ) = c 1 .Now, consider only the first K − 1 switches, and let P ′ (λ) be the polynomial only from these.If It remains to invoke singular value transformation in order to construct operators with this polynomial as their spectrum.The polynomial we have constructed has mixed parity, which is acceptable because we perform singular value transformation on a Hermitian operator.The method from [GSLW19] splits the polynomial into its even and odd parts and combines them together via a linear combination of block encodings.This introduces an extra factor of 1 2 .Since we only care about making eigenvalues either close to 0 or 1, we can use a simple version of oblivious amplitude amplification via the Chebyshev polynomial T 3 to remove this extra factor.
Lemma 35 (Construction of projectors via singular value transformation).Say H is a Hermitian matrix with −1 ≤ H ≤ 1, and we are given a block encoding U H of H. Say H has the eigendecomposition H = i λ i Π i .Say P (λ) is a degree-d polynomial that, for some δ and for certain regions Then, there exists a quantum circuit U P (H) which is a block encoding of a matrix P (H) defined by: where P is a function that satisfies λ ∈ [a Proof.Our starting point is Theorem 56 of [GSLW19], which allows us to construct a block encoding of To get rid of the factor of 1 2 , let T 3 (λ) be the third Chebyshev polynomial of the first kind, and use Corollary 18 of [GSLW19] to construct a block encoding of: Observe that −T 3 (0) = 0 and −T 3 ( 1 2 ) = 1, and also that −T 3 (δ/2) ≤ 2δ and −T 3 ((1 − δ)/2) ≥ 1 − 2δ.So if we let P (λ) := −T 3 1 2 P (λ) then −T 3 1 2 P (H) is the P (H) we wanted to construct.Since P (λ) has degree d, and −T 3 (λ) has degree O(1), the total circuit complexity and number of invocations to controlled-U H is O(d) overall.Now we are ready to synthesize the operators we need in order to prove Proposition 8, Lemma 17, and Lemma 29.
We start with the energy estimation result from Proposition 8, which is adapted from [Ral21].This work gives a construction for energy estimation that attempts to minimize the constant factors in polynomial degree as well as ancilla count.We do not concern ourselves with these, and just give an asymptotic bound that also has an extra factor of n in the complexity.A similar method also appears in [MRTC21].
Proof of Proposition 8. Say U A is a block encoding with k ancilla qubits of a Hermitian operator A. For y ∈ {0, 1} k we define A y := (⟨y| ⊗ I)U A ( |0 k ⟩ ⊗ I), so that: Then, let G A := |0 k ⟩⊗A and ḠA := y∈{0,1} k \{0 k } |y⟩⊗A y and apply a generalized Toffoli gate to synthesize the isometry: Our goal is to compute |x⟩, where x is the index of the interval [a x , b x ] containing λ.We will repeatedly use the above isometry to compute |x⟩ one bit at a time.To do so, we make use of Lemmas 34 and 35 to synthesize δ j -accurate projectors for each bit of x -call them P (j) for the approximate projector for the j'th bit.
If we apply the above construction involving U A above for each of P (j) , we construct an isometry Ẽ(M) that achieves Ẽ(M) = x |x⟩ ⊗ G Px with the desired properties.We observe that: Since the P (j) commute, we can evaluate: It remains to show that Px P (M ) − P x ≤ δ est .The j-th bit has 2 j many flips.We take δ j := δ est /2 • 2 −j such that P (j) P (M ) − P (j) ≤ δ j , with P (j) being the ideal projector onto the j'th bit.That way the total error in spectral norm is at most δ est • n+1 j=1 2 −j ≤ δ est .Each projector has complexity O κ −1 log 2 j δ −1 j .The total complexity is: Next, we construct the operator that underlies the left-right POVM.To follow this argument, we refer again to Fig. 2 which gives a sketch of L, R as well as the function approximated by Lemma 34.
Proof of Lemma 17.The operator P LR is constructed using Lemmas 34 and 35: all we need to do is select the intervals [a x , b x ] and the labels c x .Recall that the goal was to ensure that ∥P 2 LR Π i ∥ ≤ δ sup for λ i ̸ ∈ L and ∥(I − P 2 LR )Π i ∥ ≤ δ sup for λ i ̸ ∈ R. Recall Fig. 2. The intervals of interest are the gaps of L and R, for which we set c x = 0 and c x = 1 respectively.If we invoke Lemma 34 with precision δ = δ sup /6, then for λ i ̸ ∈ L we have ∥P 2 LR Π i ∥ ≤ δ 2 ≤ δ sup /3 and for λ i ̸ ∈ R we have ∥(I − P 2 LR )Π i ∥ ≤ 1 − (1 − δ) 2 ≤ δ sup /3 as desired.From the construction of L and R in Definition 16, we see that each of the intervals are exactly 2 −n−r−2 apart, and there are 2 n+r+1 of them.Since the intervals alternate with c x = 0 and c x = 1, there are O(2 n+r ) many switches.So, the implementation requires (n + r)2 n+r log δ −1 sup many invocations of the block encoding of H.
Finally, we construct the attenuation operator A (M ) and its approximation Ã(M) .We demand that A (M ) vanishes outside of the rounding promise M , and is close to 1 in inside of a truncated rounding promise M T .For all the other eigenspaces, the approximation Ã(M) and A (M ) agree.This means that we can use the construction of Ã(M) to define what A (M ) should be outside these regions.
Proof of Lemma 29.As usual, we use Lemmas 34 and 35 to construct the block encoding of Ã(M) .This construction results in a function ã(M) which lets us define a (M ) and hence A (M ) via the requirements in Lemma 29.
We need to ensure that ã(M) (λ) ≤ δ leak for λ ̸ ∈ M , and that ã(M) (λ) ≥ 1 − δ leak for λ ∈ M T .This immediately shows what intervals to use for Lemma 29: the intervals with c x = 0 are the gaps of M , and the intervals with c x = 1 are the intervals of M T .
There are 2 'flips' per interval of M , so there are O(s (M ) ) flips total.Since M T was obtained by taking M and shrinking the intervals by κγ on each side, the gap width for the purposes of Lemma 29 is κγ.Accordingly, the polynomial degree and hence the query complexity are O κ −1 γ −1 log s (M ) δ −1 .

C The Approximate Lindbladian
We describe here a method that does not try to force any rounding promises to make energy estimation unambiguous.Rather this method works directly with approximate jump operators that arise from imperfect energy estimation of the Hamiltonian H on the entire Hilbert space H.
We consider a general energy estimation unitary V (±) that acts as where λ i 's are the energies of H and Π i 's the projectors onto the corresponding eigensubspaces, and λx 's denote the possible outcomes produced by the energy estimation method.
The summation over x means that any energy estimation method necessarily produces superpositions of energy estimates.Typically, the magnitude of f (λ i , λx ) decreases with increasing distance between the true energy λ i and the approximate energies λx .For instance, when energy estimation is based on phase estimation the unitary U = exp(iH), the estimates are simply n-bit binary fractions λx = x/2 n and the amplitudes f (λ i , λx ) are given by The spread of f (λ i , λx ) can be made narrower with the help of median amplification.We now construct an oracle O L encoding jump operators of an approximate Lindbladian L using a general estimation unitary V (±) as in Eq. ( 145).We consider the case that there is only one coupling operator S. O L is given by the circuit: where G denotes a block encoding of: To describe the jump operators that are encoded by the above circuit, we define the approximate Bohr frequencies ν to be differences of the form λx − λy .The approximate jump operators L(ν) are given by where S(ν) = x,y λx− λy=ν i,j f ( λx , λ i )f ( λx , λ j )Π i SΠ j (150) = x,y λx− λy=ν A(x)SA(y) † . (151) Unfortunately, the operators that "sandwich" the coupling operator S are not projectors as in the case for Davies generators.Thus, it is not possible to interpret the operators L(ν) as jump operators of a Davies generator with respect to some Hamiltonian that is close to the original Hamiltonian.Therefore, it is much more difficult to determine the fixed point of the corresponding Lindbladian L, which is given by L(ρ) = ν G(ν) S(ν)ρ S(ν) † + 1 2 S(ν) † S(ν)ρ + ρ S(ν) † S(ν) . (153) This 'approximate Davies generator' L is challenging to analyze.It is intuitive that this Davies generator's steady state must be somewhat close to the ideal thermal state: the error due to the finite precision of energy estimation can be dealt with in the same way as for the promised Davies generators, and the 'rounding errors' stemming from energies located at x/2 n + 1/2 n+1 also only shift the accuracy of the estimate slightly.
To our knowledge, no method exists for rigorously proving a bound on the distance between this Lindblad operator's steady state and the true thermal state.A potential candidate for a proof technique via 'approximate detailed balance' appeared in [TOV + 11, CB21], where it was leveraged to prove the accuracy of the quantum Metropolis algorithm.But it is not clear how to translate this technique to Davies generators.
How severe are these rounding errors?In practice, it may be the case that the steady state of approximate Davies generator is close to the ideal thermal state.Here, we present some evidence that it may be difficult to prove such a claim without an assumption on the Hamiltonian.We construct an 'adversarial' Hamiltonian that places its eigenvalues exactly at the locations where the rounding errors are maximized.Specifically, its energies are located at: where 0 ≤ i < 2 n is an integer and α ∈ [0, 1] is an 'adversariality parameter'.When α = 0, then energy estimation is perfect and the operators A(x) are projectors exactly.But as α increases, the eigenvalues are shifted so that they are rounded up or down with increasing entropy.This Hamiltonian is specifically designed to capture these rounding errors alone: the energies are spaced out evenly with exactly the precision of the energy estimation protocol.Fig. 4 shows the error in preparing the ideal thermal state when using the approximate Davies generator with this adversarial Hamiltonian.We see that at α = 0 the method computes the thermal state exactly.But as α increases we see significant errors.If we amplify the accuracy of the energy measurement using median amplification, then the errors appear only for larger α.But no amount of amplification can remove the error for α = 1.
Since the approximate Davies generator is significantly simpler and may be less costly to implement than our scheme of random promised Davies generators, a proof technique establishing rigorous accuracy bounds on this Davies generator may result in a substantially improved algorithm for thermal state preparation.Figure 4: Accuracy of thermal state preparation using the approximate Davies generator in Eq. ( 153) with an adversarial Hamiltonian with eigenvalues as in Eq. ( 154).We perform energy estimation to 3 bits of precision on a 4 qubit system.Energy estimation is performed with median amplification, where M denotes the number of estimates over which the median is performed.We observe significant errors for large α.
When M is increased, the large errors appear only for the largest α.The software that produced these data is available at: https://github.com/qiskit-community/promised-davies-generator.
(M ) x+1 ) that spans the gap between two adjacentAccepted in Quantum 2023-10-03, click title to verify.Published under CC-BY 4.0.If we write λ ∈ M , we mean that λ is contained in the union of all the intervals.The promised eigenspace projector P (M ) x of the intervals [a x , b y ] and [a y , b y ] ∈ M , respectively.Its jump operators S (M ) α (ν) are given by

Figure 2 :
Figure 2: Sketch of the rounding promises from this section with n = 2, and r = 1.The diagram features a sketch of the function p LR (λ) that defines the projector P LR from Lemma 17, a sketch of the fine-grained rounding promises from Definition 16, and a sketch of the coarse-grained rounding promises from Definition 21.The numbers above the intervals in the coarse-grained promises indicate the index x.
[a x , b x ] and [a x+1 , b x+1 ], then we take any intervals [a ′ , b x ], [a x+1 , b ′ ] for some a ′ , b ′ , and replace them with [a ′ , b ′ ].If [a ′ , b x ] is the final interval, then we replace it with [a ′ , 1].Finally, if x = 0, then we take the first interval [a ′ , b ′ ] and replace it with [0, b ′ ].

Figure 3 :
Figure 3: Analysis of the spectral gap of the Davies generators of the transverse field Ising model from Eq. (80).We set β = 10.0 throughout.We see that the convergence rates of the promised Davies generators are a little slower than the ideal Davies generator as γ → 0. While the convergence rate slows down significantly with large γ, selecting reasonably small values of γ achieves a convergence time just as fast as the unphysical case of γ = 0.There is also significant variation among the promises M j .The software that produced these data is available at: https://github.com/qiskit-community/promised-davies-generator.