Single-ancilla ground state preparation via Lindbladians

We design a quantum algorithm for ground state preparation in the early fault tolerant regime. As a Monte Carlo-style quantum algorithm, our method features a Lindbladian where the target state is stationary, and its evolution can be efficiently implemented using just one ancilla qubit. Our algorithm can prepare the ground state even when the initial state has zero overlap with the ground state, bypassing the most significant limitation of methods like quantum phase estimation. As a variant, we also propose a discrete-time algorithm, demonstrating even better efficiency and providing a near-optimal simulation cost depending on the desired evolution time and precision. Numerical simulation using Ising models and Hubbard models demonstrates the efficacy and applicability of our method.

Zhiyan Ding

Chi-Fang Chen

Lin Lin

A promising application of quantum computers is to simulate ground state properties of quantum many-body systems [3, 26, 47, 16, 37, 29]. To concretely evaluate the end-to-end algorithmic cost, however, the state preparation problem rises as a major conceptual bottleneck [36, 27].

From a complexity theory standpoint, few-body Hamiltonian ground states can be QMA-hard to prepare [24, 2, 23, 35, 1].

Thus, we do not expect quantum computers to efficiently prepare ground states for every few-body Hamiltonian. While this worst-case hardness may not apply to practically relevant systems, it explains the theoretical obstacles towards proving algorithmic guarantees. Indeed, justifying the efficacy of most existing ground state algorithms requires additional assumptions. The most common and transparent assumption is the existence of an easy-to-prepare (pure or mixed) quantum state $\rho_0$ with a good overlap with the ground state $\ket{\psi_0}$, i.e., $p_0=\braket{\psi_{0}|\rho_0|\psi_{0}}=\Omega(1/\mathrm{poly}(n))$

where $n$ is the system size. This assumption allows us to provably solve ground state preparation problems with a cost scaling with $\poly(1/p_0)$ (among other dependences). For instance, the cost associated with quantum phase estimation (QPE) scales as $\Or(1/p^2_0)$ or $\Or(1/p_0)$ depending on the implementation [34, 24, 30], while nearly optimal post-QPE algorithms exhibit a scaling of $\Or(1/\sqrt{p_0})$ [29, 15].

Unfortunately, the above strategy demands an instance-dependent choice of good trial states and can be nontrivial to justify in practically relevant systems [27].

From a thermodynamics standpoint, however, ground state (and more generally, low-energy states) should be easy to prepare: just put the sample into
a low-temperature fridge.

Drawing from this intuition, several studies  [21, 38, 40, 6, 32, 33] explore the idea of creating an appropriate cold bath and establishing an effective system-bath coupling to efficiently cool the system to a low-energy state. Nonetheless, these methods often require a large bath and lack the theoretical framework to rigorously study their efficacy.

Guided by thermodynamic intuition, a new wave of Monte Carlo style quantum algorithm

has been proposed to extract the precise working principle behind cooling [44, 48, 42, 7, 9, 13]. The idea is to algorithmically emulate a bath by designing a Lindbladian $\mc{L}$ whose fixed point is provably the Gibbs state $\rho_{\beta}  = e^{-\beta H}/\Tr[e^{-\beta H}] $. That is, we evolve the density operator according to the master equation

\begin{align}
\frac{\ud \rho}{\ud t} = \mc{L}[\rho] \quad \text{where}\quad \mc{L}[\rho_{\beta}] \approx 0.
\end{align}

This quantum algorithm can be regarded as a cousin of the classical Markov chain
Monte Carlo (MCMC) method. As long as we prescribe the fixed point to be our
target  state, the algorithmic cost per sample is

\begin{equation}\label{eq:MC_cost}
(\text{simulation cost per unit time}) \times (\text{mixing time}).
% \begin{aligned}
%    &(\text{cost per sample})\\
%    = &(\text{simulation cost per time}) \times (\text{mixing time})
% \end{aligned}
\end{equation}

Thus, the hardness of low-energy problems ($\beta\gg 1$), in this framework, reduces to the mixing time of certain “quantum Markov chains”. The mixing time can be
very much case-dependent for classical Markov chains; we expect similar rich
behavior to transfer to the quantum case, and in particular, we do not expect
this algorithm to efficiently solve QMA-hard problems (just as we do not expect
classical Monte Carlo methods to solve NP-hard problems efficiently). While the
efficacy of the QPE and the Monte Carlo approaches both rely on additional
assumptions, we hope the quantum Markov chain can be constructed more systematically than the ansatz state. Furthermore, based on the experimental success of cooling, we might expect the mixing time for physically relevant systems to be reasonably short, i.e.,  scaling polynomially with the system size.

For a Hamiltonian $H = \sum_{i} \lambda_i
\ket{\psi_i}\bra{\psi_i}$ with spectral gap $\Delta = \lambda_1 - \lambda_0>0$,
the ground state is very close to the Gibbs state $\rho_{\beta}$ with an inverse temperature $\beta=\Omega(\Delta^{-1})$[The prefactor may depend logarithmically on the Hilbert space dimension.].
However, the preparation of the Gibbs state requires imposing certain quantum
detailed balance condition [25, 5, 9],
a delicate balance of the rates between the heating and cooling transitions.
To accommodate specific energy measurements, energy label storage, and certain coherent arithmetic operations, the algorithm uses techniques such as Quantum Fourier Transform and controlled Hamiltonian simulation, along with other complex control logic [44, 48, 9, 39].

The approach involving Gibbs state preparation presents another algorithmic challenge. While simulating the Lindblad dynamics up to time
$T$ may achieve a complexity of $\Or( T \log(T/\epsilon)/\log\log(T/\epsilon) )$ [12, 28] – achieving near-optimal complexity in $T$ and $\epsilon$ – these methods rely on yet another layer of intricate block encodings and precise control circuits, assuming fully fault-tolerant quantum computers. In the early fault-tolerant regime with constrained quantum resources, such as a limited number of logical qubits, a significant simplification would be essential.

Returning to our discussion on ground state preparation, it is important to recognize that the ground state represents the lowest-energy state where the detailed balance condition becomes singular (i.e., all heating processes are prohibited). As a result, it might not be necessary to utilize all the previously mentioned techniques to address general detailed balance conditions. This line of reasoning naturally leads us to our guiding question:

Can we devise a Monte Carlo-style quantum algorithm for ground states in the early fault-tolerant regime while maintaining near-optimal complexity in the Lindblad evolution time $T$ and precision $\epsilon$?

Main result

In this work, we introduce a Lindbladian and develop algorithms that satisfy the following salient features:

• Correctness: The ground state is a fixed point of a Lindblad evolution defined by a single jump operator.

• Convergence: The algorithm converges to the ground state even when starting with zero initial overlap, i.e., $p_0=0$, after the mixing time. While the mixing time can vary across systems, we give a random matrix example where we can grasp this quantity.

• Efficient simulation: The continuous-time Lindblad evolution can be simulated using one ancilla qubit and minimal control logic. With a discrete-time reformulation of the Lindblad evolution, the cost attains near-optimal complexity in both the evolution time and precision.

The main Lindbladian is defined as follows

\begin{equation}\label{eqn:Lindblad_dynamics}
\frac{\mathrm{d}}{\mathrm{d} t}\rho = \mc{L}[\rho] = \underset{=:\mc{L}_H[\rho]}{\underbrace{-i [H, \rho]}}+ \underset{=:\mc{L}_K[\rho]}{\underbrace{K \rho K^{\dag}-\frac12 \{K^{\dag}K,\rho\}}}\,,
\end{equation}

with one jump operator

\begin{align}
K&:=\sum_{i,j\in[N]} \hat{f}(\lambda_i-\lambda_j) \ket{\psi_i}\braket{\psi_i|A|\psi_j}\bra{\psi_j}\label{eqn:jump_freq}\\
&= \int_{-\infty}^{\infty} f(s) A(s) \ud s.%\quad\quad{\text{(time domain)}}
\label{eqn:jump_time}
\end{align}

Here $A$ is a Hermitian coupling operator (or coupling matrix) that acts on the system with its Heisenberg evolution $A(s)=e^{i Hs} A e^{-iHs}$, and represents the interaction between the system and the environment. The time domain function $f(s):=\frac{1}{2\pi}\int_{\RR} \hat{f}(\omega)e^{-i\omega s}\ud \omega$ is the inverse Fourier transform of the filter function $\hat{f}$ in the frequency domain. In the following, we elaborate on the desirable features of the proposed algorithm.

Left: Sketch of the circuit structure for one iteration of the main algorithm. The circuit utilizes only one ancilla qubit. In this graph, $C_1,\cdots,C_M$ are local interaction operators that can be designed to act on only one system qubit (or a small number of system qubits).
Right: a visual representation of how the dissipative part $\exp(\mathcal{L}_K\tau)$ transforms the density matrix in each iteration towards the ground state..

The right-hand side of (<ref>) consists of two terms: The coherent part features the Hamiltonian system $[H,\rho]$ that generates unitary dynamics; the dissipative part $\mc{L}_K$ is parameterized by the jump operator $K$. Since the Hamiltonian trivially fixes its own ground state, it suffices to focus on the dissipative part; the fixed-point property remains valid with or without the coherent part[Nevertheless, our numerical findings suggest that incorporating the coherent part often leads to a significant reduction in the mixing time.].  In the frequency domain, the jump operator $K$ is expressed as the coupling matrix $A$ in the energy basis weighted by the filter function $\hat{f}(\omega)$, which depends on the difference in energy $\omega = \lambda_i-\lambda_j$. This function should be real, nonnegative, and satisfy the following condition (detailed assumptions are given in <ref> in Supplemental Materials):

\begin{equation}\label{eqn:f_equal_to_zero}
\hat{f}(\omega)=0\quad \text{for each}\quad \omega\geq 0\,.
\end{equation}

Under <ref>, we obtain $\braket{\psi_i|K|\psi_j}=0$ when $\lambda_i\geq \lambda_j$. Intuitively, the jump operator $K$ forbids energy increments and favors a decrease in energy.
Remarkably, such a simple condition guarantees that the ground state is an exact fixed point of the Lindbladian. To see this, we first notice that

\begin{equation}\label{eqn:K_fix_ground}
K\ket{\psi_0}=\sum_{i} \hat{f}(\lambda_i-\lambda_0) \ket{\psi_i}\braket{\psi_i|A|\psi_0}=0\,,
\end{equation}

which implies $\mc{L}_K[\ket{\psi_0}\bra{\psi_0}]=0$.

Together with $\mc{L}_H[\ket{\psi_0}\bra{\psi_0}]=0$, we  conclude that $\ket{\psi_0}\bra{\psi_0}$ is a fixed point of the Lindblad dynamics (<ref>). That is,

\begin{align}
\quad e^{(\mc{L}_H+\mc{L}_K) t}[\ket{\psi_0}\bra{\psi_0}] =\ket{\psi_0}\bra{\psi_0}, \quad t \ge 0.
\end{align}

Of course, there might exist other fixed points (i.e., the map may not be ergodic), and the proper choice of $A$ that ensures ergodicity should be case-dependent.

In practice, the coupling matrix $A$ can often be chosen as a simple operator (e.g., a local Pauli operator); it is difficult to prove convergence in a general setting. We consider a simplified scenario where $A$ is represented by a random matrix with independent entries in the energy eigenbasis of $H$. In this case, we give a partial argument for ergodicity.

Let $\rho(t)$ be the solution to the Lindblad dynamics (<ref>) and $A_{i,j}=\left\langle \psi_i\right|A\ket{\psi_j}$. Assume for any $t\geq 0$, $A$ is independently drawn from a probability distribution[To accurately represent its dependence on time, $A$ should be denoted as $A_t$. However, for simplicity and consistency with other notations, we will omit the subscript $t$ in this theorem.] on the set of Hermitian matrices such that $\mathbb{E}(A_{i,j})=0$.
Then starting from a diagonal matrix  $\rho(0)$ in the eigenbasis of $H$,
$\rho^\star=\ket{\psi_0}\bra{\psi_0}$ is the unique fixed point of the Lindblad dynamics in the expectation sense. In particular, given any observable $O$, $\lim_{t\rightarrow\infty} \mathbb{E}(\mathrm{Tr}(O\rho(t)))=\bra{\psi_0}O\ket{\psi_0}$, where the expectation is taken on the randomness of $A$.

See the complete result in the Supplemental Material (<ref>). Strictly speaking, the expected operator $\mathbb{E}(\rho(t))$ is not the density operator $\rho(t)$ that we store in the quantum memory but still gives some optimistic intuition about ergodicity (see Supplemental Section <ref>).

Technically, taking expectations over independent entries of $A$ substantially simplifies the transition matrix. This independent assumption can also be seen as a version of the Eigenstate Thermalization Hypothesis (ETH) [43, 14], which incorporates additional assumptions on the variance of $A_{i,j}$. In fact, ETH has been employed to explain finite-time thermalization in chaotic open quantum systems [42, 7]. Under stronger ETH-type assumptions, one may be able to prove the convergence for $\rho(t)$ instead of $\mathbb{E}(\rho(t))$ as in [7], but we merely focus on the much simpler object $\mathbb{E}(\rho(t))$ without distracting from the presentation of the algorithm.

The convergence rate of a Markov chain process is characterized by the mixing time, which describes the time scale at which any two input states become close to each other (see <ref>).

Analogous to classical Monte Carlo sampling, we do not know a priori the mixing time

associated with the jump operator $K$ and the Hamiltonian $H$; we expect this to be system-dependent. In fact, proving rapid mixing for Lindbladians

is a highly challenging problem [22, 4], and it is not the primary objective of our work.
However, in an exercise to understand the convergence of our method (in Appendix <ref>), we show that under additional assumptions on the coupling matrix $A$ and the eigenvalue distributions of $H$, the Lindblad dynamics discussed in our study can achieve polynomial mixing time.

Simulating continuous-time dynamics

Quantum computers in the early fault-tolerant regime often have few ancilla qubits and minimal control logic.
Unlike Hamiltonian simulation, where product formulas (or Trotterization) are well known to satisfy the above constraints, the case of Lindbladian simulation has received comparatively less attention. In this study, we present an efficient algorithm for simulating (<ref>) that uses a single ancilla qubit and simple controlled gates as follows:

For any Hamiltonian $H$ and coupling operator $A$, there exists a quantum algorithm simulating the continuous Lindblad dynamics  <ref> using one ancilla qubit. For simulation time $T$ and precision $\epsilon$, the total cost in terms of Hamiltonian simulation time is $ T_{H,\mathrm{total}}=\widetilde{\Theta}\left((1+\|H\|)\Delta^{-1}T^{2+o(1)}\epsilon^{-1-o(1)}\right)$.

The structure of our simulation circuit is sketched in <ref>, and its detailed derivation is given in <ref> in the Supplemental Materials. In the above results, the cost of simulating the coupling matrix $A$ is mild because $A$ is often much simpler to simulate than the (global) Hamiltonian $H$.
A more comprehensive version of the above theorem (<ref>), including the cost based on the number of controlled-$A$ gate queries, can be found in the Supplemental material (Appendix <ref>). The cost also depends on the specific form of the filter function $\hat{f},f$ that we choose. More specifically, when the function $f$ have bounded $N$-th derivatives, the term $o(1)=N^{-1}$, which diminishes as the smoothness of $f$ increases (as $N$ approaches infinity). Additionally, the constant in $T_{H,\mathrm{total}}$ also depends on the upper bound of $\left|\mathrm{d}^N \hat{f}/\mathrm{d}\omega^N\right|$.

We expect the most resource-intensive part of the algorithm to be Hamiltonian simulation. Thus, we quantify the cost through the total Hamiltonian simulation time, represented as $T_{H,\mathrm{total}}$. For an end-to-end cost analysis, one may further Trotterize the Hamiltonian simulation subroutine and analyze its discretization error, and there is no hidden block-encoding or extra ancilla involved (i.e., no controlled Hamiltonian simulation and only controlled-$A$, which is much easier to implement).

The complexity of our algorithm approximately scales quadratically with the simulation time $T$ (a reasonable choice is the mixing time $T\sim \tmix$) and inversely with the error $\epsilon$. This scaling reminds us of a first-order product formula, which has a second-order error $\exp((\mathcal{L}_H+\mathcal{L}_K)\tau))=\exp(\mathcal{L}_H\tau)\exp(\mathcal{L}_K\tau)+\mathcal{O}(\tau^2)$. We emphasize that the Lindbladian simulation qualitatively differs from the Hamiltonian simulation: even to arrive at the $\wt{\Or}(T^2/\epsilon)$ scaling requires a nontrivial argument to approximately implement $\exp(\mathcal{L}_K\tau)$.

Simulating discrete-time dynamics

Is it possible to reduce the “first-order” scaling $T^2/\epsilon$ of the cost?
For Hamiltonian simulation, we can use higher-order formulas, but here we cannot invert dissipation (i.e., simulating $\exp(-\mathcal{L}_K\tau)$), which leads to fundamental constraints on the permissible product formula sequences. This forces us to use a small time step size. What if we shift our perspective and interpret

\begin{equation}\label{eqn:discrete_informal}
\rho_{n+1}=\exp(\mathcal{L}_H\tau)\exp(\mathcal{L}_K\tau)\rho_n=\mathcal{N}_{\tau}(\rho_n)
\end{equation}

not as an approximation for the continuous dynamics expressed by <ref>, but by itself a discrete-time dynamics? Indeed, our primary goal is not the precise simulation of Lindbladian dynamics but rather to prepare the ground state.
A key observation is that the ground state is a stationary state of the discretized dynamics for any value of $\tau$. Unlike the continuous dynamics simulation, which requires  $\tau=\Or(\epsilon/\norm{H}T)$, we could now potentially choose a large step size $\tau$. In such a scenario, <ref> might not approximate the continuous Lindblad dynamics due to the discretization error, yet it still fixes the ground state.
Simulating this discrete-time dynamics for $M$ discrete steps corresponds to a simulation for an effective time $T=M\tau$. By choosing $\tau=\Or(1)$, we obtain a scheme with a significantly improved scaling with respect to both the simulation time $T$ and precision $\epsilon$. Remarkably, the total Hamiltonian simulation time of the discrete-time dynamics no longer depends on $\Delta/\norm{H}$ as in <ref>, but only on $\Delta$. This can be significantly faster when simulating large-scale systems or Hamiltonians involving unbounded operators (such as differential operators in first quantization).

For any Hamiltonian $H$ and coupling operator $A$, there exists a quantum algorithm simulating the discrete-time Lindblad dynamics  <ref> using one ancilla qubit. For simulation time $T$ and precision $\epsilon$, the total Hamiltonian simulation time is  $T^\star_{H,\mathrm{total}}=\widetilde{\Theta}(\Delta^{-1}T^{1+o(1)}\epsilon^{-o(1)})$.

See the Supplemental Materials (<ref>) for the formal statement. Compared with the continuous-time Lindblad simulation, this algorithm exhibits a total Hamiltonian simulation time that scales nearly linearly with $T$ and sub-polynomially with $1/\epsilon$.[While adopting a high-order Trotter formula does not reduce the cost of continuous-time simulation, it does benefit the discrete-time dynamics.] In <ref>, the appearance of $o(1)$ stems from adopting a high-order Trotter formulation of order $p$, which contributes an exponent of the form $1/p=o(1)$ as $p\to \infty$.

The algorithm in <ref> approximates the discrete-time Lindblad dynamics (<ref>). Theoretically, the “effective” mixing time $\tmix'=M_{\mathrm{mix}}\tau$
of the discrete dynamics may not be the same as the continuous-time Lindblad dynamics $\tmix$. On the other hand, numerical observations indicate that the mixing time of both dynamics can be comparable, but the discrete-time version requires significantly less total Hamiltonian simulation time, reducing the total cost.

Numerics results

We numerically illustrate the performance of our algorithms using the one-dimensional transverse field Ising model (TFIM) model defined on $6$ sites (see the detailed setup, including a concrete example of $\hat{f}$ and $f$, and additional numerical results on TFIM and Hubbard models, in Appendix <ref>).

In our numerical simulations, we set $\tau=1$ and $r=2$ for discrete-time Lindblad simulation ($r$ is the number of segments that we use to approximately simulate $\exp(\mathcal{L}_K\tau)$ in the discrete-time Lindblad simulation algorithm, see <ref>),
while for continuous-time Lindblad simulation, we use $\tau=0.1$. We start with an initial state with zero overlap ($\bra{\psi_0}\rho_0\ket{\psi_0}\approx10^{-17}$), repeat each Lindblad dynamics simulation 100 times, and compute the average energy and overlap with the ground state. The results are depicted in Figure <ref>, which demonstrates the effectiveness of the Lindblad dynamics in reducing energy to the ground state and generating a state that closely aligns with the ground state. The trajectories of the continuous and discrete-time dynamics noticeably differ from each other, yet they can both prepare the ground state, and the mixing times are comparable. Additionally, the discrete-time Lindblad dynamics ($\tau=1,r=2$) demands almost one-tenth of the total Hamiltonian simulation time required by the approximated continuous-time Lindblad dynamics ($\tau=0.1$).

[Hamiltonian simulation time vs energy]
[Hamiltonian simulation time vs overlap]

[Lindblad simulation time vs overlap]
Performance of continuous versus discrete-time Lindblad dynamics for preparing the ground state for the TFIM with $6$ sites. The continuous-time simulation uses a small time step $\tau=0.1$. The trajectory of the discrete-time dynamics (using a large time step $\tau=1$) deviates from that of the continuous-time Lindblad dynamics, but it successfully prepares the ground state and is more efficient than the continuous-time dynamics. Here, the Hamiltonian simulation time refers to the sum of the Hamiltonian simulation $t$ in all $e^{\pm iHt}$ subroutines used in the circuit.

Key ideas

The design of our specific Monte Carlo-style quantum algorithms draws inspiration from the recent work for the preparation of thermal states [9], which considers a set of jump operators labeled by $a,\omega$ as $\{\sqrt{\gamma(\omega)}\hat{A}^a(\omega)\}_{a,\omega}$.
Intuitively, $\{\hat{A}^a(\omega)\}_{a,\omega}$ corresponds to the transitions of $A^a$ with an energy difference of approximately $\omega$. The weights $\gamma(\omega)$ depend on the temperature and are designed to satisfy the approximate quantum detailed balance condition. A key observation of this work is that such a condition for ground state preparation can be considerably simplified using a single jump operator.

To reduce the end-to-end ancilla usage, we first  define a dilation $\tilde{K}:=\ket{1}\bra{0}\otimes K+\ket{0}\bra{1}\otimes K^{\dag}$. Because $K$ is the single jump operator, we have
$e^{\mc{L}_{K} t}[\rho] \approx \Tr_a e^{-i \wt{K}\sqrt{t}}\left[\ket{0}\bra{0}\otimes\rho \right]e^{i \wt{K}\sqrt{t}}+\mathcal{O}(t^2)$, which reduces the Lindbladian simulation to the dilated Hamiltonian simulation.  Here, $\Tr_a$ denotes tracing out the ancilla qubit.

To implement $e^{-i \wt{K}\sqrt{t}}$, we observe that $K$ is expressed in integral form, and its integrand also involves Hamiltonian simulation. Consequently, the direct discretization and Trotterization of the right-hand side of (<ref>) pose challenges since the exponent in each Trotter step contains another Hamiltonian simulation $e^{-iHs}$.

The technical challenges rise due to the currently forbidding cost of implementing general block encodings and a limited number of control ancillas. This is analogous to the story between the two very different algorithms for Hamiltonian simulation: linear combinations of unitary (LCU) [11] / quantum singular value transformation (QSVT) [18], and product formulas [46, 10]. The former is theoretically transparent and helps us understand the asymptotic complexity of the Hamiltonian simulation, but practitioners often turn to the latter because of its simple implementation.

Our method to overcome this difficulty is motivated by the following observation. When $\tau$ is sufficiently small,

\begin{equation}
\begin{split}
&e^{i \tau\int_{0}^S f(s) A(s) \ud s}
= \mc{T} e^{i \int_0^S \tau f(s) A(s) \ud s}+\mathcal{O}(\tau^2)\\
=& e^{i H S} \mc{T} e^{i  \int_0^S  (\tau f(s) A-H)\ud s}+\mathcal{O}(\tau^2).\label{eq:time-order}
\end{split}
\end{equation}

where $\mc{T}$ denotes time-ordering for the exponential.

We recognize the time-ordered exponential as the interaction picture over Hamiltonian $H$, which leads to the second expression. Even in cases where the filter $f(s)$ is not real, applying the dilatation and trace-out technique allows us to realize $\tau f(s)A-H$ by a time-dependent Hamiltonian simulation. Then, we can further discretize and Trotterize it into simple gates and only one ancilla qubit is needed. More carefully, our actual implementation utilizes a second-order Trotter formula (not the first-order formula in (<ref>)), which ensures that the local truncation error of simulating $\exp(-i\wt{K}\sqrt{t})$ is of the order $\mathcal{O}(t^2)$ rather than $\mathcal{O}(t)$. For a more comprehensive description of the implementation, see Appendix <ref>.

This paper presents a Monte Carlo approach for ground state preparation of a quantum Hamiltonian using a single ancilla qubit. Merely a single jump operator ensures that the system's ground state is a fixed point of the Lindbaldian dynamics. Unlike many ground state preparation algorithms, such as QPE, the performance of this algorithm is mainly determined by the mixing time and may work even with zero initial overlap. Similar to classical Monte Carlo methods, the relationship between mixing time and system size, the practical selection of the coupling matrix $A$, and the optimal choice of the filtering function $f$ may be system-dependent.
For instance, optimizing these choices to prepare the ground state of strongly correlated chemical systems like FeMoCo [27], which are currently inaccessible via adiabatic state preparation, might offer a compelling avenue of research.
This material is based upon work supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator (Z.D.).

Additional support is acknowledged from the Challenge Institute for Quantum Computation (CIQC) funded by National Science Foundation (NSF) through grant number OMA-2016245, the  Applied Mathematics Program of the US Department of Energy (DOE) Office of Advanced Scientific Computing Research under contract number DE-AC02-05CH1123, and a Google
Quantum Research Award (L.L.). L.L. is a Simons investigator.   We thank Joao Basso, Garnet Chan, Toby Cubitt, Yulong Dong, Xiantao Li, Subhayan Roy Moulik, Yu Tong for helpful discussions.

[1]
D. Aharonov, D. Gottesman, S. Irani, and J. Kempe.
The power of quantum systems on a line.
Comm. Math. Phys., 287(1):41–65, 2009.

[2]
D. Aharonov and T. Naveh.
Quantum NP-a survey.
arXiv preprint quant-ph/0210077, 2002.

[3]
A. Aspuru-Guzik, A. D. Dutoi, P. J. Love, and M. Head-Gordon.
Simulated quantum computation of molecular energies.
Science, 309(5741):1704–1707, 2005.

[4]
Á. Capel, C. Rouzé, and D. S. França.
The modified logarithmic sobolev inequality for quantum spin systems:
classical and commuting nearest neighbour interactions.
arXiv preprint arXiv:2009.11817, 2020.

[5]
E. A. Carlen and J. Maas.
Gradient flow and entropy inequalities for quantum markov semigroups
with detailed balance.
Journal of Functional Analysis, 273(5):1810–1869, 2017.

[6]
M. Cattaneo, M. A. Rossi, G. García-Pérez, R. Zambrini, and
S. Maniscalco.
Quantum simulation of dissipative collective effects on noisy quantum
PRX Quantum, 4:010324, 2023.

[7]
C.-F. Chen and F. G. S. L. Brandão.
Fast thermalization from the eigenstate thermalization hypothesis.
arxiv/2112.07646, 2023.

[8]
C.-F. Chen, H.-Y. Huang, R. Kueng, and J. A. Tropp.
Concentration for random product formulas.
PRX Quantum, 2:040305, 2021.

[9]
C.-F. Chen, M. J. Kastoryano, F. G. S. L. Brandão, and A. Gilyén.
Quantum thermal state preparation.
arXiv/2303.18224, 2023.

[10]
A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu.
Theory of trotter error with commutator scaling.
Phys. Rev. X, 11(1):011020, 2021.

[11]
A. M. Childs and N. Wiebe.
Hamiltonian simulation using linear combinations of unitary
arXiv/1202.5822, 2012.

[12]
R. Cleve and C. Wang.
Efficient quantum algorithms for simulating Lindblad evolution.
In 44th International Colloquium on Automata, Languages, and
Programming, ICALP 2017, July 2017.

[13]
T. S. Cubitt.
Dissipative ground state preparation and the dissipative quantum
arXiv/2303.11962, 2023.

[14]
L. D'Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol.
From quantum chaos and eigenstate thermalization to statistical
mechanics and thermodynamics.
Advances in Physics, 65(3):239–362, 2016.

[15]
Y. Dong, L. Lin, and Y. Tong.
Ground-state preparation and energy estimation on early
fault-tolerant quantum computers via quantum eigenvalue transformation of
unitary matrices.
PRX Quantum, 3:040305, 2022.

[16]
Y. Ge, J. Tura, and J. I. Cirac.
Faster ground state preparation and high-precision ground energy
estimation with fewer qubits.
Journal of Mathematical Physics, 60(2):022202, 2019.

[17]
M. Gevrey.
Sur la nature analytique des solutions des équations aux
dérivées partielles. Premier mémoire.
Annales scientifiques de l'École Normale Supérieure, 3e
série, 35:129–190, 1918.

[18]
A. Gilyén, Y. Su, G. H. Low, and N. Wiebe.
Quantum singular value transformation and beyond: exponential
improvements for quantum matrix arithmetics.
In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory
of Computing, pages 193–204, 2019.

[19]
V. Gorini, A. Kossakowski, and E. C. G. Sudarshan.
Completely positive dynamical semigroups of N‐level systems.
Journal of Mathematical Physics, 17(5):821–825, 2008.

[20]
L. Hormander.
Analysis of linear partial differential operators I:
distribuction theory and fourier analysis.
New York: Springer-Verlag, 1990.

[21]
D. B. Kaplan, N. Klco, and A. Roggero.
Ground states via spectral combing on a quantum computer.
arxiv/1709.08250, 2017.

[22]
M. J. Kastoryano and F. G. Brandao.
Quantum gibbs samplers: The commuting case.
Communications in Mathematical Physics, 344:915–957, 2016.

[23]
J. Kempe, A. Kitaev, and O. Regev.
The complexity of the local Hamiltonian problem.
SIAM J. Comput., 35(5):1070–1097, 2006.

[24]
A. Y. Kitaev, A. Shen, and M. N. Vyalyi.
Classical and quantum computation.
American Mathematical Soc., 2002.

[25]
A. Kossakowski, A. Frigerio, V. Gorini, and M. Verri.
Quantum detailed balance and KMS condition.
Commun. Math. Phys., 57(2):97–110, 1977.

[26]
B. P. Lanyon, J. D. Whitfield, G. G. Gillett, M. E. Goggin, M. P. Almeida,
I. Kassal, J. D. Biamonte, M. Mohseni, B. J. Powell, M. Barbieri,
A. Aspuru-Guzik, and A. G. White.
Towards quantum chemistry on a quantum computer.
Nature Chemistry, 2(2):106–111, 2010.

[27]
S. Lee, J. Lee, H. Zhai, Y. Tong, A. M. Dalzell, A. Kumar, P. Helms, J. Gray,
Z.-H. Cui, W. Liu, et al.
Evaluating the evidence for exponential quantum advantage in
ground-state quantum chemistry.
Nature Comm., 14(1):1952, 2023.

[28]
X. Li and C. Wang.
Simulating Markovian open quantum systems using higher-order series
arXiv/2212.02051, 2022.

[29]
L. Lin and Y. Tong.
Near-optimal ground state preparation.
Quantum, 4:372, 2020.

[30]
L. Lin and Y. Tong.
Heisenberg-limited ground state energy estimation for early
fault-tolerant quantum computers.
PRX Quantum, 3:010318, 2022.

[31]
G. Lindblad.
On the generators of quantum dynamical semigroups.
Communications in Mathematical Physics, 48(2):119 – 130, 1976.

[32]
M. Metcalf, E. Stone, K. Klymko, A. F. Kemper, M. Sarovar, and W. A. de Jong.
Quantum Markov chain Monte Carlo with digital dissipative
dynamics on quantum computers.
Quantum Science and Technology, 7(2):025017, 2022.

[33]
X. Mi, A. A. Michailidis, S. Shabani, K. C. Miao, P. V. Klimov, J. Lloyd,
E. Rosenberg, R. Acharya, I. Aleiner, T. I. Andersen, M. Ansmann, F. Arute,
K. Arya, A. Asfaw, J. Atalaya, J. C. Bardin, A. Bengtsson, G. Bortoli,
A. Bourassa, J. Bovaird, L. Brill, M. Broughton, B. B. Buckley, D. A. Buell,
T. Burger, B. Burkett, N. Bushnell, Z. Chen, B. Chiaro, D. Chik, C. Chou,
J. Cogan, R. Collins, P. Conner, W. Courtney, A. L. Crook, B. Curtin, A. G.
Dau, D. M. Debroy, A. D. T. Barba, S. Demura, A. D. Paolo, I. K. Drozdov,
A. Dunsworth, C. Erickson, L. Faoro, E. Farhi, R. Fatemi, V. S. Ferreira,
L. F. B. E. Forati, A. G. Fowler, B. Foxen, E. Genois, W. Giang, C. Gidney,
D. Gilboa, M. Giustina, R. Gosula, J. A. Gross, S. Habegger, M. C. Hamilton,
M. Hansen, M. P. Harrigan, S. D. Harrington, P. Heu, M. R. Hoffmann, S. Hong,
T. Huang, A. Huff, W. J. Huggins, L. B. Ioffe, S. V. Isakov, J. Iveland,
E. Jeffrey, Z. Jiang, C. Jones, P. Juhas, D. Kafri, K. Kechedzhi, T. Khattar,
M. Khezri, M. Kieferova, S. Kim, A. Kitaev, A. R. Klots, A. N. Korotkov,
F. Kostritsa, J. M. Kreikebaum, D. Landhuis, P. Laptev, K. M. Lau, L. Laws,
J. Lee, K. W. Lee, Y. D. Lensky, B. J. Lester, A. T. Lill, W. Liu,
A. Locharla, F. D. Malone, O. Martin, J. R. McClean, M. McEwen, A. Mieszala,
S. Montazeri, A. Morvan, R. Movassagh, W. Mruczkiewicz, M. Neeley, C. Neill,
A. Nersisyan, M. Newman, J. H. Ng, A. Nguyen, M. Nguyen, M. Y. Niu, T. E.
OBrien, A. Opremcak, A. Petukhov, R. Potter, L. P. Pryadko, C. Quintana,
C. Rocque, N. C. Rubin, N. Saei, D. Sank, K. Sankaragomathi, K. J. Satzinger,
H. F. Schurkus, C. Schuster, M. J. Shearn, A. Shorter, N. Shutty, V. Shvarts,
J. Skruzny, W. C. Smith, R. Somma, G. Sterling, D. Strain, M. Szalay,
A. Torres, G. Vidal, B. Villalonga, C. V. Heidweiller, T. White, B. W. K.
Woo, C. Xing, Z. J. Yao, P. Yeh, J. Yoo, G. Young, A. Zalcman, Y. Zhang,
N. Zhu, N. Zobrist, H. Neven, R. Babbush, D. Bacon, S. Boixo, J. Hilton,
E. Lucero, A. Megrant, J. Kelly, Y. Chen, P. Roushan, V. Smelyanskiy, and
D. A. Abanin.
Stable quantum-correlated many body states via engineered
arxiv/2304.13878, 2023.

[34]
M. A. Nielsen and I. Chuang.
Quantum computation and quantum information.
Cambridge Univ. Pr., 2000.

[35]
R. Oliveira and B. M. Terhal.
The complexity of quantum spin systems on a two-dimensional square
arXiv preprint quant-ph/0504050, 2005.

[36]
P. J. J. O'Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean,
R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen,
Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero,
A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank,
A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven,
A. Aspuru-Guzik, and J. M. Martinis.
Scalable quantum simulation of molecular energies.
Phys. Rev. X, 6:031007, 2016.

[37]
T. E. O’Brien, B. Tarasinski, and B. M. Terhal.
Quantum phase estimation of multiple eigenvalues for small-scale
(noisy) experiments.
New Journal of Physics, 21(2):023022, 2019.

[38]
S. Polla, Y. Herasymenko, and T. E. O'Brien.
Quantum digital cooling.
Phys. Rev. A, 104:012414, Jul 2021.

[39]
P. Rall, C. Wang, and P. Wocjan.
Thermal state preparation via rounding promises.
arxiv/2210.01670, 2022.

[40]
B. Rost, L. D. Re, N. Earnest, A. F. Kemper, B. Jones, and J. K. Freericks.
Demonstrating robust simulation of driven-dissipative problems on
near-term quantum computers.
arXiv/2108.01183, 2021.

[41]
M. B. RUSKAI.
Beyond strong subadditivity? improved bounds on the contraction of
generalized relative entropy.
Reviews in Mathematical Physics, 06(05a):1147–1161, 1994.

[42]
O. Shtanko and R. Movassagh.
Preparing thermal states on noiseless and noisy programmable quantum
arXiv/2112.14688, 2023.

[43]
M. Srednicki.
The approach to thermal equilibrium in quantized chaotic systems.
Journal of Physics A: Mathematical and General, 32(7):1163,

[44]
K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, and F. Verstraete.
Quantum Metropolis sampling.
Nature, 471(7336):87–90, 2011.

[45]
L. N. Trefethen and J. A. C. Weideman.
The exponentially convergent trapezoidal rule.
SIAM Review, 56(3):385–458, 2014.

[46]
H. Trotter.
On the product of semi-groups of operators.
Proc. Amer. Math. Soc., 10:545, 1959.

[47]
L. Veis and J. Pittner.
Quantum computing applied to calculations of molecular energies: CH2
The Journal of Chemical Physics, 133(19):194106, 11 2010.

[48]
M.-H. Yung and A. Aspuru-Guzik.
A quantum–quantum Metropolis algorithm.
Proceedings of the National Academy of Sciences,
109(3):754–759, 2012.

Supplemental Material

§ NOTATION, FACTS, AND ORGANIZATION

In the supplementary material, we use capital letters for matrices and curly font for superoperators. Besides the usual $\Or$ notation, we use the following asymptotic notations: we write $f=\Omega(g)$ if $g=\Or(f)$; $f=\Theta(g)$ if $f=\Or(g)$ and $g=\Or(f)$; $f=\wt{\Or}(g)$ if $f=\Or(g\operatorname{polylog}(g))$. We use $\|\cdot\|$ to denote vector or matrix 2-norm: when $v$ is a vector we denote by $\|v\|$ its 2-norm, and when $A$ is matrix we denote by $\|A\|$ its operator norm (or the Schatten $\infty$-norm). Given any Hermitian matrix $H$ and a positive integer $N$, there exists a constant $C_N$ such that

\begin{equation}\label{eqn:matrix_fact}
\left\|\exp(-i Ht)-\sum^N_{n=0}\frac{(-iH)^n}{n!}t^n\right\|\leq C_N\|H\|^{N+1}t^{N+1}\,.
\end{equation}

The trace norm (or the Schatten $1$-norm) of a matrix $A$ is $\norm{A}_1=\Tr\left[\sqrt{A^{\dag}A}\right]$. A useful  inequality that we often use in our proof is that

\begin{equation}\label{eqn:trace_norm_inequ}
\|AB\|_1\leq \|A\|\|B\|_1\,.
\end{equation}

for any two matrices $A,B$ such that $AB$ is well defined.

Given a superoperator $\mathcal{L}$ that acts on operators, the induced $1$-norm is

\begin{equation}
\|\mathcal{L}\|_1:=\sup_{\|\rho\|_1\leq 1}\|\mathcal{L}(\rho)\|_1\,.
\end{equation}

According to GKLS theorem [19, 31], if $\mathcal{L}$ is a Lindbladian, which has the form

\[
\mathcal{L}(\rho)=-i [H, \rho]+\sum^M_{m=1}K_m \rho K^{\dag}_m-\frac12 \{K^{\dag}_mK_m,\rho\}\,,
\]

then $\exp(\mathcal{L}t)$ is a quantum channel, i.e., it is a completely positive trace-preserving (CPTP) map. It is also contractive under trace distance [41].
For any two density operators $\rho_1,\rho_2$, and any $t>0$,

\begin{equation}\label{eqn:lindblad_contractive}
\|\exp(\mathcal{L}t)\rho_1-\exp(\mathcal{L}t)\rho_2\|_1 \leq \|\rho_1-\rho_2\|_1.
\end{equation}

The following lemma provides an upper bound on the first-order Trotter error between two superoperators:
Given two superoperators $\mathcal{L}_1,\mathcal{L}_2$ such that $\exp(\mathcal{L
}_1 t)$, $\exp(\mathcal{L
}_2 t)$, and $\exp((\mathcal{L
}_1+\mathcal{L}_2) t)$ are quantum channels for all $t>0$. Then for all $t>0$,

\begin{equation}\label{eqn:operator_fact}
\left\|\exp((\mathcal{L
}_1+\mathcal{L}_2) t)-\exp(\mathcal{L
}_1 t)\exp(\mathcal{L
}_2 t)\right\|_1=\mathcal{O}\left(\|[\mathcal{L}_1,\mathcal{L}_2]\|_1t^2\right).
\end{equation}

This proof is similar to the proof of <cit.>. Define

\[
\mathcal{E}(t)=\exp(\mathcal{L
}_1 t)\exp(\mathcal{L
}_2 t)\,.
\]

Then, we obtain

\[
\frac{\ud \mathcal{E}(t)}{\ud t}=(\mathcal{L
}_1 t),\mathcal{L
}_2 t)\,,
\]

which implies

\[
\mathcal{E}(t)=\exp((\mathcal{L
}_1+\mathcal{L}_2) t)+\int^t_0\exp((\mathcal{L
}_1 \tau),\mathcal{L
}_2 \tau)\ud \tau\,.
\]

Because $\|[\exp(\mathcal{L
}_1 \tau),\mathcal{L
}_2]\|_1=\mathcal{O}\left([\mathcal{L}_1,\mathcal{L}_2]\tau\right)$ and $\|\exp((\mathcal{L
}_2 \tau)\|_1=1$, we obtain that

\[
\left\|\mathcal{E}(t)-\exp((\mathcal{L
}_1+\mathcal{L}_2) t)\right\|_1=\mathcal{O}\left([\mathcal{L}_1,\mathcal{L}_2]t^2\right)\,,
\]

which proves (<ref>).

Next, we bound the $1$-norm of $\mc{L}_K$ corresponding to the jump operator $K$ in <ref> of the Lindblad dynamics in <ref>:

Assume $\int^\infty_{-\infty}|f(t)|\ud t=\mathcal{O}(1)$, then

\begin{equation}\label{eqn:K_bound}
\|\mathcal{L}_K\|_1=\mathcal{O}(\|K\|^2)=\mathcal{O}(\|A\|^2)\,.
\end{equation}

To prove (<ref>), we first use (<ref>) to obtain $\|\mathcal{L}_K\|_1=\mathcal{O}(\|K\|^2)$. Then

\[
\|K\|\leq \|A\|\int^\infty_{-\infty}|f(s)|\ud s=\mathcal{O}(\|A\|).
\]

Finally, to characterize the convergence of the Lindblad dynamics, we define the mixing time as follows:

We define the continuous and discrete mixing time as follows.

* For any Lindbladian $\mc{L}$, the mixing time $\tmix$ is the smallest time for which

\begin{align*}
\norm{e^{\mc{L} \tmix }[\rho-\rho']}_{1} \le \frac{1}{2} \norm{\rho-\rho'}_{1} \quad \text{for any states}\quad \rho, \rho'.
\end{align*}

* For any quantum channel $\mathcal{N}$, the “discrete” mixing time $M_{\mathrm{mix}}$ is the smallest integer $M$ for which

\begin{align*}
\norm{\mathcal{N}^M[\rho-\rho']}_{1} \le \frac{1}{2} \norm{\rho-\rho'}_{1} \quad \text{for any states}\quad \rho, \rho'.
\end{align*}

In particular, if the Lindblad $\mc{L}$ has mixing time $\tmix$, any initial state must be $\epsilon$ close to the ground state after time $\tmix \cdot \log_2(2/\epsilon)$ since the trace distance between any two quantum states is at most $2$.

The rest of the supplementary material is organized as follows:

* Section <ref> introduces our algorithms for the approximate preparation of the ground state.

* Section <ref> analyzes the costs associated with these algorithms.

* Section <ref> discusses the convergence properties of the continuous-time Lindblad dynamics under the Eigenstate Thermalization Hypothesis (ETH) type ansatz.

* Section <ref> presents details of numerical results, including additional simulations with the TFIM and Hubbard models.

§ THE SIMULATION ALGORITHM

In this section, we first introduce the filter function needed to define the jump operator (Section <ref>) and present our algorithm for the approximate preparation of the ground state by simulating the continuous-time Lindblad dynamics (Section <ref>). Subsequently, we introduce the discrete-time Lindblad dynamics (Section <ref>). These algorithms offer the advantage of requiring only a single ancilla qubit and can be implemented without the need for constructing a block encoding of $K_s$.

§.§ Choice of filter function

As mentioned earlier, ensuring that the ground state remains a fixed point of the Lindblad dynamics requires $\hat{f}$ to satisfy the condition stated in <ref>.
Furthermore, as discussed in Key ideas, in order to numerically simulate $\exp(\mathcal{L}_K(t))$, we need to discretize and truncate the integral $\int^\infty_{-\infty} f(s)A(s)\ud s$. To ensure that this truncation does not introduce significant errors, additional assumptions must be imposed on $\hat{f}$ such that $f(t)$ decays rapidly as $|t|$ approaches infinity. In summary, we adopt the following assumption regarding $\hat{f}$:
[Filter function in the frequency domain]
Given parameter $S_\omega>2\Delta$, where $\Delta=\lambda_1-\lambda_0$. We assume the filter function $\hat f(\omega)\in \CC^{\infty}(\RR)$ satisfies

\begin{align}\label{eqn:f_hat_assumption}
0\le \hat f(\omega)\le 1 \quad \text{and}\quad \hat f(\omega) = \begin{cases}
\Omega(1) &\quad \omega = -\Delta\\
0 &\quad \omega \notin [-S_{\omega},0]
\end{cases}\quad \text{for all}\quad \omega\in\RR.
\end{align}

* Given any $N>0$, there exists a constant $C_N$ such that

\[
\left|\frac{\rd^N \hat{f}(\omega)}{\rd \omega^N}\right|\leq C_N\Delta^{-N}\,,\quad \text{for all}\quad \omega\in\mathbb{R}\,.
\]

The choice of $S_{\omega}$ significantly impacts the ease of discretizing the time integral (<ref>) and the mixing time of the Lindblad dynamics. Striking the right balance is crucial. Choosing a larger $S_\omega$ facilitates the transition of weight from high energy to low energy, leading to a reduction in the mixing time. However, it also introduces challenges in discretizing the time integral (<ref>), increasing the complexity of simulating the Lindblad dynamics. An example of such filter function that satisfies <ref> can be found in Section <ref>, where we set $S_\omega=\mathcal{O}(\|H\|)$.

We present two results on the Fourier transform of $\hat{f}(\omega)$

\begin{equation}\label{eqn:fomega_fourier}
f(s):=\frac{1}{2\pi}\int_{\RR} \hat{f}(\omega)e^{-i\omega s}\ud \omega
\end{equation}

This is slightly different from the standard convention of the Fourier transform, but it agrees with the convention used in Ref. [9].

Let $\hat{f}(\omega)$ satisfy <ref>. Then $f(s)$ in <ref> can be extended into an entire function denoted by $f(z)$. Given $N>0$,

\begin{equation}\label{eqn:f_super_decay}
\abs{f(s)}=\mathcal{O}\left(C_NS_\omega\Delta^{-N}(1+|s|)^{-N}\right)\, ,\quad \text{for all}\ s\in\RR,
\end{equation}

\begin{equation}\label{eqn:f_int_bound}
\sup_{|\omega|\leq 2\|H\|,\ |\zeta|\leq 1}\int^\infty_{\infty} \left|f(s+i\zeta)\exp(i\omega (s+i\zeta))\right|\ud s=\mathcal{O}\left(C_NS_\omega\Delta^{-N}\exp(2\|H\|+S_\omega)\right).
\end{equation}

Here $C_N$ is the constant from <ref>.

According to the Paley–Wiener–Schwartz theorem <cit.>, $f(s)$ can be extended to become an entire function in the complex plane. First, when $|z|\leq 1$, we have

\[
\left|f(z)\right|=\left|\frac{1}{2\pi}\int^0_{-S_\omega} \hat{f}(\omega)e^{-i\omega z}\ud \omega\right|\leq \frac{1}{2\pi}\int^0_{-S_\omega} \left|\hat{f}(\omega)e^{-i\omega z}\right|\ud \omega=\mathcal{O}\left(C_NS_\omega\exp(S_\omega |\mathrm{Im}(z)|)\right)\,.
\]

When $|z|>1$,

\[
\begin{aligned}
\left|f(z)\right|&=\left|\frac{1}{2\pi}\int^0_{-S_\omega} \hat{f}(\omega)e^{-i\omega z}\ud \omega\right|=\frac{1}{2\pi}\left|\int^0_{-S_\omega} \frac{\rd^N \hat{f}(\omega)}{\rd \omega^N}\frac{1}{\left(-iz\right)^N}e^{-i\omega z}\ud \omega\right|\\
&\leq \frac{1}{2\pi}\int^0_{-S_\omega} \frac{1}{|z|^N}\left|\frac{\rd^N \hat{f}(\omega)}{\rd \omega^N}\right|\left|e^{-i\omega z}\right|\ud \omega=\mathcal{O}\left(C_NS_\omega\Delta^{-N}|z|^{-N}\exp(S_\omega |\mathrm{Im}(z)|)\right)\,.
\end{aligned}
\]

For any $N\in \NN$, there exists a constant $C'_N>0$ such that

\begin{equation}\label{eqn:f_time_bound}
\abs{f(z)}=\mathcal{O}\left(C_NS_\omega\Delta^{-N}(1+|z|)^{-N} \exp(S_\omega \abs{\Im z})\right), \quad \text{for all}\ z\in\CC,
\end{equation}

which proves (<ref>) by choosing $z$ to be a real number.

Plugging (<ref>) into left-hand side (<ref>), we obtain that

\[
\int^\infty_{\infty} \left|f(s+i\zeta)\exp(i\omega (s+i\zeta))\right|\ud s\leq \mathcal{O}\left(C_NS_\omega\Delta^{-N}\exp(2\|H\|+S_\omega)\int^\infty_{\infty} (1+|s|)^{-N} \ud s\right)
\]

for any $|\omega|\leq 2\|H\|,\ |\zeta|\leq 1$. This proves (<ref>).

The conclusion of Lemma <ref> helps us control the error of the trapezoid rule:

Suppose $h$ is analytic in the strip $|\mathrm{Im}(z)|<1$. Suppose
\[\lim_{|s|\rightarrow\infty}\sup_{|\zeta|<1}|h(s+i\zeta)|\rightarrow0\,,\]

\[
\int^\infty_{-\infty}\left|h(s+i\zeta)\right|\ud s\leq M
\]

for any complex numbers $s$ and $|\zeta|<1$. Then, for any $\tau_s>0$, we have

\[
\left|\int^\infty_{-\infty}h(s)\ud s-\tau_s\sum^{\infty}_{n=-\infty}h(n\tau_s)\right|\leq \frac{2M}{\exp(2\pi/\tau_s)-1}\,.
\]

§.§ Simulating continuous-time Lindblad dynamics with one ancilla qubit

The Lindblad consists of the coherent part and the dissipative part. For simplicity, we use a first-order Trotter splitting

\begin{align}
e^{\mc{L}t} &= e^{(\mc{L}_H +\mc{L}_K)t} \\
& \approx (e^{\mc{L}_H\tau}e^{\mc{L}_K\tau} )^{t/\tau}\quad \text{for time step} \quad \tau.
\end{align}

This allows us to focus on simulating $e^{\mathcal{L}_K\tau}$ for a short time $\tau$.

§.§ Simulation of Lindblad dynamics

We define the dilated Hermitian jump operator using one ancilla qubit as

\begin{equation}
\wt{K}:=\begin{pmatrix}
0 & K^{\dag}\\
K & 0
\end{pmatrix}.
\end{equation}

Define the partial trace $\Tr_a\left(\sum^1_{i,j=0}\ket{i}\bra{j}\otimes \rho_{i,j}\right)=\sum^1_{i=0}\rho_{i,i}$, which traces out the ancilla qubit.  Our simulation approach is based on the following lemma:

\begin{align}
\sigma(t):=\Tr_a e^{-i \wt{K}\sqrt{t}} \left[\ket{0}\bra{0}\otimes\rho \right]e^{i \wt{K}\sqrt{t}}\label{eqn:simplied_lindblad_simulation}.
\end{align}

Then for a short time $t\ge 0$,

\begin{align}\label{eqn:inexact_lindblad}
\|\sigma(t)-\rho(t)\|_1=\mathcal{O}(\|K\|^4t^2)
\end{align}

where $\rho(t)$ is the solution to the Lindblad dynamics

\begin{equation}
\partial_t \rho(t)=\mc{L}_K[\rho(t)]\quad \text{with initial condition}\quad \rho(0) = \rho.
\end{equation}

According to <ref>, we first obtain

\[
\|\exp(-i\wt{K}\sqrt{t})-(1-i\wt{K}\sqrt{t}+\wt{K}^2t/2+\wt{K}^3t^{3/2}/6)\|=\mathcal{O}(\|\wt{K}\|^4t^2)=\mathcal{O}(\|K\|^4t^2)\,.
\]

Let $T_{3,K}=(1-i\wt{K}\sqrt{t}+\wt{K}^2t/2+\wt{K}^3t^{3/2}/6)$ denote the truncated Taylor expansion, and $1$ denotes the identity operator. We have

\begin{equation}\label{eqn:first_bound}
\left\|\sigma(t)-\Tr_a T_{3,K} \left[\ket{0}\bra{0}\otimes\rho \right]T_{3,K}^\dagger\right\|_1=\mathcal{O}(\|K\|^4t^2)\,.
\end{equation}

Because $\Tr_a (\ket{i}\bra{j}\otimes \rho)=\delta_{i,j}\rho$, We  note that terms with an odd power of $\Or(t^{1/2})$ vanish after applying the partial trace operation. This implies

\begin{equation}\label{eqn:second_bound}
\begin{aligned}
&\left\|\Tr_a T_{3,K} \left[\ket{0}\bra{0}\otimes\rho \right]T_{3,K}^\dagger-\left(1+K\rho K^\dagger t-\frac{1}{2}(K^\dagger K\rho t+\rho K^\dagger K t)\right)\right\|_1\\
=&\left\|\Tr_a T_{3,K} \left[\ket{0}\bra{0}\otimes\rho \right]T_{3,K}^\dagger-(1+\mathcal{L}_Kt)\rho\right\|_1\\=
\end{aligned}
\end{equation}

By combining equations (<ref>) and (<ref>), we arrive at the following expression:

\[
\|\sigma(t)-(1+\mathcal{L}_Kt)\rho\|_1=\mathcal{O}(\|K\|^4t^2)\,.
\]

Finally, since

\[
\|\rho(t)-(1+\mathcal{L}_Kt)\rho\|_1=\|\exp(\mathcal{L}_K t)-(1+\mathcal{L}_Kt)\rho\|_1=\mathcal{O}\left(\|\mathcal{L}_K\|^2_1t^2\right)=\mathcal{O}(\|K\|^4t^2)\,,
\]

we conclude the proof.

According to <ref>, we can approximate $e^{\mathcal{L}_K\tau}[\rho]$ by employing a dilated Hamiltonian simulation, tracing out the ancilla qubit, and resetting the ancilla qubit to the state $\ket{0}$.

§.§ Trotter-based simulation without block encoding the jump operator

Simulating <ref> still requires a block encoding for the dilation $\wt{K}$, which can require many ancilla qubits. Now we demonstrate that it is possible to use Trotter splitting to implement $\exp(-i\widetilde{K}\sqrt{t})$ directly using one ancilla qubit in total and without the block encoding of $K$ or $\widetilde{K}$.

Discretizing the time integral. As discussed before, to implement $ \exp \left(i \int_{-\infty}^\infty f(s) A(s) ds\right)$, we need to truncate and discretize the integral.

First, since $\hat{f}(\omega)$ is assumed to be compactly supported in <ref>, according to Appendix <ref> Lemma <ref>,  $f(s)$ decays super-polynomially as $\abs{s}\to \infty$. Thus, we can choose a suitable $S_s>0$, restrict the integration range to $[-S_s,S_s]$, and discretize this interval using a uniform grid $s_l=l\tau_s$ and $l=-M_s,\ldots,M_s$, where $\tau_s=S_s/M_s$. We can then approximate the integral using a trapezoidal rule:

\begin{equation}\label{eqn:K_s}
K_{s}:=  \sum_{l=-M_s}^{M_s} f(s_l)  e^{iH s_l}Ae^{-iH s_l} w_l,
\end{equation}

where $w_l=\tau_s/2$ for $l=\pm M_s$, and $w_l=\tau_s$ for $-M_s<l<M_s$.

Let $\hat{f}$ satisfy <ref>. Then, given any $N,\epsilon'>0$, if

\[
S_s=\Omega\left(\frac{1}{\Delta}\left(\frac{C_NS_\omega\|A\|}{\Delta\epsilon'}\right)^{\frac{1}{(N-1)}}\right),\quad \tau_s=\mathcal{O}\left(\left(\|H\|+S_\omega+\log\left(\frac{C_NS_\omega\norm{A}}{\Delta\epsilon'}\right)\right)^{-1}\right)\,,\]

we have

\begin{equation}\label{eqn:error_quadrature}
\|K-K_{s}\|=\mathcal{O}(\epsilon')\,.
\end{equation}

Here the constant $C_N$ comes from <ref>.

The proof is given in Appendix <ref>.
By taking $N\to \infty$ in <ref>, we obtain $S_s=\Or(\epsilon'^{-o(1)})$, the number of quadrature points $M_s=S_s/\tau_s=\widetilde{\Or}(\epsilon'^{-o(1)})$. Therefore the accuracy of the trapezoidal rule improves super polynomially with respect to the increase of $S_s,M_s$. We remark that by imposing stronger regularity conditions on $\hat{f}(\omega)$ (e.g., $\hat{f}(\omega)$ is compactly supported and is in the Gevrey class [17], which is stronger than $\hat{f}\in C^\infty(\RR)$), the aforementioned result can be enhanced to $S_s=M_s=\Or(\polylog(\epsilon'^{-1}))$. For simplicity, this work does not incorporate such an improvement.

According to Lemma <ref>, a small number of point $M_s$ already ensures a good approximation

\begin{equation}
\wt{K}\approx \wt{K}_s=
\begin{pmatrix}
0 & \left(\sum_{l=-M_s}^{M_s} f(s_l) A(s_l) w_l\right)^{\dag}\\
\sum_{l=-M_s}^{M_s} f(s_l) A(s_l) w_l & 0
\end{pmatrix}\\
=:\sum_{l=-M_s}^{M_s} \wt{H}_l.
\end{equation}

Since the Heisenberg evolution $A(s)=e^{iHs}Ae^{-iHs}$ is Hermitian, the term further factorizes as

\begin{equation}
\wt{H}_l=\begin{pmatrix}
0 & f^*(s_l) A(s_l) w_l\\
f(s_l) A(s_l) w_l& 0
\end{pmatrix}=\sigma_l\otimes A(s_l) \quad \text{where}\quad \sigma_l:=w_l(\sigma_x \Re f(s_l) + \sigma_y\Im f(s_{l}) )
\end{equation}

and $\sigma_x$ and $\sigma_y$ are Pauli matrices.

Second-order Trotter splitting.
After discretizing the time labels, the next step is to Trotterize the Hamiltonian evolution

\begin{align}
e^{-i\sqrt{\tau} \wt{K}_s} = e^{-i\sqrt{\tau} \sum_{l} \wt{H}_{l}}.
\end{align}

Specifically, we employ the second-order Trotter formula to balance between efficiency and accuracy. The second-order Trotter formula for $\exp(-i\sqrt{\tau} \widetilde{K})$ can be expressed as:

\begin{equation}
\rightprod_{l} e^{-i \frac{\sqrt{\tau}}{2} \wt{H}_l}\leftprod_{l} e^{-i \frac{\sqrt{\tau}}{2} \wt{H}_l}-e^{-i\sqrt{\tau} \wt{K}_s}=\tau^{3/2}\sum_{l_1,l_2,l_3}a_{l_1,l_2,l_3}\wt{H}_{l_1}\wt{H}_{l_1}\wt{H}_{l_3}+\Or(\tau^2),
\label{eq:eta3}
\end{equation}

where the coefficients $a_{l_1,l_2,l_3}$ can be calculated from Taylor expansion, and the left and right-ordered products are defined by

\begin{equation}\label{eqn:secondordertaylor_collect}
\leftprod_{l} e^{-i \sqrt{\tau} \wt{H}_{l}}:=e^{-i \sqrt{\tau} \wt{H}_{M_s}} \cdots e^{-i \sqrt{\tau} \wt{H}_{-M_s}} \quad \text{and} \quad
\rightprod_{l} e^{-i \sqrt{\tau} \wt{H}_l}:=e^{-i \sqrt{\tau} \wt{H}_{-M_s}} \cdots e^{-i \sqrt{\tau} \wt{H}_{M_s}}.
\end{equation}

The Trotter error bounds (<ref>) and Lemma <ref> together give an approximation scheme for simulating the Lindblad dynamics in (<ref>). For any initial state $\rho$, we have that

\begin{equation}
\begin{aligned}
&\mathrm{Tr}_a\left(\rightprod_{l} e^{-i \frac{\sqrt{\tau}}{2} \wt{H}_l}\leftprod_{l} e^{-i \frac{\sqrt{\tau}}{2} \wt{H}_l} \ket{0}\bra{0}\otimes\rho \rightprod_{l} e^{i \frac{\sqrt{\tau}}{2} \wt{H}_l} \leftprod_{l} e^{i \frac{\sqrt{\tau}}{2} \wt{H}_l}\right)\\
=&\mathrm{Tr}_a\left(e^{-i\sqrt{\tau} \wt{K}_s} [\ket{0}\bra{0}\otimes\rho ]e^{i\sqrt{\tau} \wt{K}_s}\right)+\Or(\tau^2)\\
=&e^{\mathcal{L}_{K_s}\tau}[\rho]+ \Or(\tau^2)\approx e^{\mathcal{L}_{K}\tau}[\rho]+ \Or(\tau^2)\,,
\end{aligned}
\end{equation}

where the second equality is obtained by using $\mathrm{Tr}_a\left(\wt{H}_{l_1}\wt{H}_{l_2}\wt{H}_{l_3}\ket{0}\bra{0}\otimes\rho\right)=0$ and $\mathrm{Tr}_a\left(\ket{0}\bra{0}\otimes\rho \wt{H}^\dagger_{l_3}\wt{H}^\dagger_{l_2}\wt{H}^\dagger_{l_1}\right)=0$. We also use Lemma <ref> in the last equality.

From the analysis above, we find that replacing the second-order formula with the first-order formula results in a local truncation error of $\Or(\tau)$. This substitution leads to an excessively large global error of $\Or(1)$.  While higher-order Trotter can further suppress the error in the first equality, the main error in the last equality comes from Lindbladian simulation by means of the method in <ref> and is not easily reducible. Therefore, the second-order Trotter seems to be adequate for the purposes of Lindbladian simulation[There exists another method to symmetrize the first-order formula to reduce the local truncation error to $\Or(\tau^2)$, which will not be discussed here.].

Canceling out back-and-forth Hamiltonian evolution.
Finally, to efficiently implement the products $e^{-i\frac{\sqrt{\tau}}{2} \widetilde{H}_l}$, we notice

\begin{align}
\exp\left(-i \frac{\sqrt{\tau}}{2} \sigma_l\otimes A(s_l)
\right)
= (I\otimes e^{iH s_{l}})\underset{=:\wt{A}_l(\sqrt{\tau})}{\underbrace{e^{-i \frac{\sqrt{\tau}}{2} \sigma_l\otimes A }}}(I\otimes e^{-iH s_{l}}).
\end{align}

Since $A$ is a simple operator (e.g., a local Pauli), the cost of $\wt{A}_l(\sqrt{\tau})$ is mostly negligible compared to that of the simulation of the system Hamiltonian.

What about the $e^{iH s_{l}}$ terms? A moment of thought reveals that we may rewrite the consecutive product in a form that efficiently cancels out the back-and-forth Hamiltonian evolution. We present the most abstract form to emphasize the simplicity of this observation.

Consider a time-order product at discretized times $s_l= l \tau_s $ for $l = - M_s,\cdots, M_s$. Then, for any Hamiltonian $H$ and a set of matrices $A_l$ depending on $l$,

\begin{align}
\rightprod_{l} A_{l}(s_{l}) &= e^{-iH S_s} \left(\rightprod_{l} A_l e^{iH \tau_s}\right) e^{iH (S_s+\tau_s)}\\
\leftprod_{l} A_{l}(s_{l}) &= e^{-iH (S_s+\tau_s)} \left(\leftprod_{l} e^{-iH \tau_s}A_l \right) e^{iH S_s}
\end{align}

where $A_l(s) = e^{i H s} A_l e^{-i Hs}$ denotes the Heisenberg evolution of $A_l$ with $H$.

The Hamiltonian evolution from two consecutive steps nearly cancels (for the right-ordered product for example)

\begin{align}
A_{l}(s_{l})A_{l}(s_{l+1}) = e^{i H s_l} A_{l} e^{i H \tau_s} A_{l} e^{-i H s_{l+1}} \quad \text{since}\quad e^{-iH s_l}e^{iH s_{l+1}}=e^{iH \tau_s}.
\end{align}

For our usage, our second-order formula consists of both the left and right products, which is

\begin{equation}
\begin{split}
&\rightprod_{l}   (I\otimes e^{iH s_l}) \wt{A}_l(\sqrt{\tau}) (I\otimes e^{-iH s_l})\leftprod_{l} (I\otimes e^{iH s_l}) \wt{A}_l(\sqrt{\tau})  (I\otimes e^{-iH s_l})\\
=&  (I\otimes e^{-iH S_s}) \underset{=:W(\sqrt{\tau})}{\underbrace{\left(\rightprod_{l} \wt{A}_l(\sqrt{\tau}) (I\otimes e^{iH \tau_s})\right)
\left(\leftprod_{l}  (I\otimes e^{-iH \tau_s})\wt{A}_l(\sqrt{\tau})\right)}} (I\otimes e^{iH S_s}).
\end{split}\label{eqn:W_eta}
\end{equation}

where $W(\sqrt{\tau})$ is a product of short-time Hamiltonian simulation.

Even nicer, the long-time simulation $(I\otimes e^{-iH S_s})$ from the previous time step exactly cancels with
$(I\otimes e^{iH S_s})$ from the subsequent time step. Therefore, we may remove both long-time evolution steps and define
the quantum channel

\begin{equation}
\mc{W}(\tau)[\rho]:=\Tr_a W(\sqrt{\tau}) \left[\ket{0}\bra{0}\otimes\rho\right] W^{\dag}(\sqrt{\tau}).
\end{equation}

The total simulation time of the system Hamiltonian for implementing the quantum channel $\mc{W}(\tau)$ now becomes

\begin{equation}\label{eqn:total_H_1}
\sum_{l=-M_s}^{M_s} \tau_s= \Or(S_s).
\end{equation}

In summary, the single ancilla simulation of the Lindblad dynamics takes the form

\begin{align}
&\rho_{m+1}=\mc{W}(\tau)[\rho_m]\quad\quad\quad\quad\quad\quad\quad\text{(purely dissipative)}\label{eqn:single_lindblad_shorttime}\\
\text{or}\quad &\rho_{m+1}=e^{\mathcal{L}_H\tau}\mc{W}(\tau)[\rho_m]\quad\quad\quad\quad\quad\text{(Trotterizing the coherent part)} \label{eqn:single_lindblad_shorttime_2}
\end{align}

The second line displays the option to include the coherence term via Trotter, which may introduce additional errors [The dissipative part itself already fixes the ground state. However, numerically, without this coherent term, we observe that the algorithm stops at other fixed points.]. The quantum circuits for $W(\sqrt{\tau})$ and (<ref>) are drawn in Figure <ref>.

[column sep=0.3cm]
$\ket{0}$     [2]A_-M_s(√(τ))        [2]A_-M_s+1(√(τ))   0em$\cdots$   [2]A_M_s-1(√(τ))      [2]A_M_s(√(τ))

[2]$\ket{\psi}$        [2]e^-iHτ_s       0em$\cdots$      [2]e^iHτ_s

[alternate]   [alternate]   [alternate]   [alternate]0em$\cdots$   [alternate]   [alternate]   [alternate]   [alternate]   [alternate]

Quantum circuits for the algorithm to simulate the continuous-time Lindblad dynamics. Upper: Quantum circuit used for $W(\sqrt{\tau})$; Lower: Quantum circuit used for continuous-time Lindblad simulation. The time step should satisfy $\tau=\Or(\epsilon)$. The measurement result of the ancilla qubit is discarded and the ancilla qubit is reset to $\ket{0}$ after each measurement. Here we assume $A$ is a local operator that only acts on one system qubit (or a small number of system qubits).

Finally, let $T= M_t \tau$ and allow $\tau$ to approach zero. We expect that $\rho_{M_t}$ from (<ref>) converges to $\exp(-\mathcal{L}_H S_s)[\rho(T)]=e^{iHS_s}\rho(T)e^{-iHS_s}$, where $\rho(t)$ is the solution of the modified Lindblad dynamics

\begin{equation}\label{eqn:modified_lindblad}
\partial_t \rho(t)=\mc{L}_{H}[\rho(t)]+\mc{L}_{K}[\rho(t)],\ \text{where}\ \rho(0)=\exp(\mathcal{L}_H S_s)[\rho_I].
\end{equation}

Compared to the original Lindblad dynamics in <ref>, the initial state is changed to $\exp(\mathcal{L}_H S_s)[\rho_I]=e^{-iHS_s}\rho_I e^{iHS_s}$. In particular, if $\rho_I$ commute with $H$ (e.g., $\rho_I$ is the density operator corresponding to an eigenstate of $H$, or the maximally mixed state), then $\exp(\mathcal{L}_H S_s)[\rho_I]=\rho_I$. At the end of the simulation, if $\rho(T)$ is the ground state $\rho_g$, then   $\exp(-\mathcal{L}_H S_s)[\rho_g]=\rho_g$. We, therefore, expect the behavior of the modified Lindblad dynamics to be very similar to that of the original dynamics.

§.§ Discrete-time Lindblad dynamics with one ancilla qubit

The short time propagator for simulating the continuous-time Lindblad dynamics in <ref>  is limited to first-order accuracy, which is mainly limited by the inexactness of the propagator according to <ref> (also discussed after <ref>). However, our goal is not to simulate the continuous-time Lindblad dynamics but to prepare the ground state. If we apply the dilated jump operator $\wt{K}$ to $\ket{0}\bra{0}\otimes \rho_g$, then according to <ref>,

\begin{equation}
\wt{K}\ket{0}\bra{0}\otimes \rho_g=\begin{pmatrix}
0\\ K\rho_g
\end{pmatrix}=\begin{pmatrix}
\end{pmatrix}.
\end{equation}

This implies

\[
\exp\left(-i\wt{K}t\right)[\ket{0}\bra{0}\otimes \rho_g]\exp\left(i\wt{K}t\right)=\ket{0}\bra{0}\otimes \rho_g\,,
\]

for any $t>0$. Thus, when the quadrature error and the Trotter error are properly controlled, we have

\begin{equation}
\exp(\mathcal{L}_H\tau)\mc{W}_a(\tau)[\rho_g]\approx \rho_g, \quad \text{for all}\  \tau>0.
\end{equation}

where $\mc{W}_a(\tau)[\rho_g]=\mathrm{Tr}_a\left(\exp\left(-i\wt{K}\sqrt{\tau}\right)[\ket{0}\bra{0}\otimes \rho_g]\exp\left(i\wt{K}\sqrt{\tau}\right)\right)$.

Therefore, the simulation scheme outlined in <ref> could conceivably be used for ground state preparation with an arbitrarily large time step $\tau$. By choosing a large time step, the quantum state $\rho_m$ may not necessarily approximate the exact dynamics $\rho(m\tau)$. However, we expect that after the mixing time $\tmix'=M_{\mathrm{mix}}\tau$, $\rho_{M_{\mathrm{mix}}}$ becomes a good approximation to $\rho_g$ because of the fixed point argument. Here, the “discrete" mixing time $M_{\mathrm{mix}}$ is given in <ref>.
If $\tmix'$ is not much larger than $\tmix$, we can significantly reduce the simulation cost thanks to the increase of $\tau$. This gives rise to the discrete-time Lindblad dynamics, which is defined as the continuous-time Lindblad simulation with a large time step.

When using a large $\tau$, additional attention must be paid to controlling the Trotter error. Specifically, $e^{-i \wt{K}\sqrt{\tau}}$ should be simulated as

\begin{equation}
e^{-i \wt{K}\sqrt{\tau}}=\left(e^{-i \wt{K}\sqrt{\tau}/r}\right)^r\approx \left(W(\sqrt{\tau}/r)\right)^r,
\end{equation}

where $r$ denotes the number of segments and should be properly chosen.

In summary, one single step of the discrete-time Lindblad dynamics is

\begin{equation}\label{eqn:accelerated_lindblad_shorttime}
\rho_{m+1}=\exp(\mathcal{L}_H\tau)\mc{W}(\tau,r)[\rho_m].
\end{equation}

\begin{equation}
\mc{W}(\tau,r)[\rho]=\Tr_a [W(\sqrt{\tau}/r)]^r \left(\ket{0}\bra{0}\otimes\rho\right)  [(W(\sqrt{\tau}/r))^{\dag}]^r\,.
\end{equation}

The quantum circuit for (<ref>) is drawn in Figure <ref>.

Since we also cancel the long-time simulation when formulating $W(\tau,r)$ in this case, if we increase $r$ and fix $T,\tau$, we expect that $\rho_{M_t}$ from
(<ref>) converges to the solution $\rho^a_{M_t}$ of the following discrete dynamics:

\begin{equation}\label{eqn:exact_acc_Lind}
\rho^a_{m+1}=\exp(\mathcal{L}_H\tau)\mathrm{Tr}_a\left(e^{-i \wt{K}\sqrt{\tau}}[\ket{0}\bra{0}\otimes\rho^a_m]e^{i \wt{K}\sqrt{\tau}}\right),\quad\text{where}\ \rho^a_0=\exp(\mathcal{L}_H S_s)[\rho_0]\,.
\end{equation}

Quantum circuit for discrete-time Lindblad dynamics simulation. The measurement result of the ancilla qubit is discarded and the ancilla qubit is reset to $\ket{0}$ after each measurement. The circuit for $W$ is given in <ref>. The time step $\tau$ can be chosen to be independent of the precision $\epsilon$, and the number of segments $r$ should be properly chosen to control the Trotter error.

§ ANALYSIS OF THE ALGORITHM

In this section, we investigate the complexity of simulating the Lindblad dynamics as described by equations (<ref>) and (<ref>) using the schemes (<ref>) and (<ref>). The rigorous version of Theorems <ref> and <ref> are presented in Theorems <ref> and <ref>, respectively.

§.§ Simulation cost of continuous-time Lindblad dynamics

Given a stopping time $T>0$ and the time step $\tau$, we assume the number of iterations $M_t=T/\tau$ is an integer. We now prove that $\rho_{M_t}$ from (<ref>) converges to $\rho(T)$ from (<ref>) with first order accuracy in $\tau$.

Given a stopping time $T>0$ and the accuracy $\epsilon>0$, we assume that $\hat{f}$ satisfies <ref>.
Then, to obtain $\left\|e^{iHS_s}\rho(T)e^{-iHS_s}-\rho_{M_t}\right\|_{1}\leq \epsilon$, we can choose

\[
S_s=\Theta\left(\frac{1}{\Delta}\left(\frac{C_NS_\omega\|A\|T}{\Delta\epsilon}\right)^{\frac{1}{(N-1)}}\right),\quad \tau_s=\Theta\left(\left(\|H\|+S_\omega+\log\left(\frac{C_NS_\omega\norm{A}T}{\Delta\epsilon}\right)\right)^{-1}\right),\  \tau=\Theta\left(\frac{\epsilon}{(\|A\|^4+\|A\|^2\|H\|)T}\right)\,,
\]

where $C_N$ is the constant that comes from <ref>.
In particular, the total Hamiltonian simulation time for $H$

\[T_{H,\mathrm{total}}=\Theta(T\tau^{-1}S_s

In addition, the total number of controlled-$A$ gates is

\[
\]

Note that the Hamiltonian simulation time is independent of the step size $\tau_s$ for discretizing the integral in forming the jump operator $K$.

The proof of <ref> is in Appendix <ref>.
We emphasize that the bottleneck leading to the first-order accuracy is not due to the Trotter splitting scheme used in <ref>, but the first-order accuracy of the Lindblad simulation method in <ref>.

§.§ Simulation cost of discrete-time Lindblad dynamics

The simulation cost of (<ref>) is shown in the following theorem:
Assume that $\hat{f}$ satisfies <ref>. Given a stopping time $T>0$, a time step $\tau=\mathcal{O}(\|A\|^{-2})$ such that $M_t=T/\tau\in \mathbb{N}$, we generate $\rho^a_{M_t}$ and $\rho_{M_t}$ using (<ref>) and (<ref>) correspondingly.

Given the accuracy $\epsilon>0$, to obtain $\|e^{iHS_s}\rho^a_{M_t}e^{-iHS_s}-\rho_{M_t}\|_1\leq \epsilon$, we can choose

\begin{equation}\label{eqn:accelerate_parameter_choice}
S_s=\Theta\left(\frac{1}{\Delta}\left(\frac{C_NS_\omega\|A\|T}{\Delta\epsilon}\right)^{\frac{1}{(N-1)}}\right),\quad \tau_s=\Theta\left(\left(\|H\|+S_\omega+\log\left(\frac{C_NS_\omega\norm{A}T}{\Delta\epsilon}\right)\right)^{-1}\right),\quad r=\Theta\left(\frac{\|A\|T^{1/2}}{\epsilon^{1/2}}\right)\,,
\end{equation}

where $C_N$ is a constant that only depends on $N,\hat{f}(\omega)$.
In particular, letting $N\to \infty$, the total Hamiltonian simulation time for $H$

\[T_{H,\mathrm{total}}=\Theta(T\tau^{-1}S_s

and the total number of controlled-$A$ evolution gates

\[
\]

We put the proof in Appendix <ref>.

Unlike the simulation of the continuous-time Lindblad dynamics in <ref>, the bottleneck of the second order accuracy in <ref> is due to the second order Trotter formula for the short time propagator $W(\sqrt{\tau})$. Replacing the second-order Trotter formula with a $p$-th order Trotter formula in defining $W(\sqrt{\tau})$, we can further improve the asymptotic scaling to be nearly linear in $T$. This is shown in <ref>.

Under the same assumptions of <ref>, but assume that the propagator $W(\sqrt{\tau}/r)$ is constructed using a $p$-th order Trotter method, we may choose

\[
\]

Then with the same choice of $S_s,\tau_s$ as in <ref>, and letting $N,p\to \infty$, the total Hamiltonian simulation time for $H$

\[T_{H,\mathrm{total}}=\Theta(T\tau^{-1}S_s

In addition, the total number of controlled-$A$ gates is

\[
\]

While <ref> partitions the coherent and dissipative components of the continuous-time Lindblad dynamics in a manner similar to a first-order Trotter method, this scheme still maintains the ground state. Therefore this first-order-like splitting may influence the mixing time, but may not contribute to the error of the ground state.

On the other hand, the accuracy of $W(\sqrt{\tau})$, which approximates $e^{-i\wt{K}\sqrt{\tau}}$ up to a change of frame, plays a vital role in the fixed point argument and needs to be implemented accurately. In practical scenarios (Sec. Numerics results and <ref>), we observe that employing a second-order Trotter method with $\tau=\Or(1)$ and $r=\Or(1)$ is often sufficient.

§.§ Proof of Lg

To prove <ref>, we first show the following proposition:
Under the assumptions of <ref>, we assume $\tau=\mathcal{O}(\|A\|^{-2})$, then

\begin{equation}\label{eqn:discretization_error}
\begin{aligned}
\left\|e^{iHS_s}\rho(T)e^{-iHS_s}-\rho_{M_t}\right\|_{1}=\mathcal{O}\left(\|K-K_{s}\|\|A\|T\right)+\mathcal{O}\left((\|A\|^4+\|A\|^2\|H\|)T\tau\right).
\end{aligned}
\end{equation}

We first notice that

\[
\left\|e^{iHS_s}\rho(T)e^{-iHS_s}-\rho_{M_t}\right\|_{1}= \left\|\rho(T)-e^{-iHS_s}\rho_{M_t}e^{iHS_s}\right\|_{1}\,.
\]

Define $\bar{\rho}_m=e^{-iHS_s}\rho_{m}e^{iHS_s}$. Then

\[
\bar{\rho}_{m+1}=e^{-iH\tau} \Tr_a \left(e^{-i H S_s} W(\sqrt{\tau}) e^{i H S_s}\right)\left(\ket{0}\bra{0}\otimes\bar{\rho}_m\right) \left(e^{-i H S_s} W(\sqrt{\tau}) e^{i H S_s}\right)^{\dag} e^{iH\tau}
\]

with $\bar{\rho}_0=e^{-iHS_s}\rho_{I}e^{iHS_s}$.

Because the Lindblad dynamics is contractive in trace distance according to <ref>, to prove (<ref>), it suffices to prove

\begin{equation}\label{eqn:one_step_error}
\begin{aligned}
\left\|\rho(\tau)-\bar{\rho}_1\right\|_{1}=\mathcal{O}\left(\|K-K_{s}\|\|A\|\tau\right)+\mathcal{O}\left((C^4_N\|A\|^4+\|H\|^2)\tau^2\right)
\end{aligned}
\end{equation}

\[
\widetilde{\rho}_1=e^{-iH\tau}\left(\Tr_a e^{-i\sqrt{\tau}\wt{K}} \ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}}\right)e^{iH\tau}\,,\]

\[\wt{\rho}_{1,s}=e^{-iH\tau}\left(\Tr_a e^{-i\sqrt{\tau}\wt{K}_{s}} \ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}_{s}}\right)e^{iH\tau}\,.\]

\begin{equation}\label{eqn:error_separation_2}
\left\|\rho(\tau)-\rho_1\right\|_{1}\leq \left\|\rho(\tau)-\widetilde{\rho}_1\right\|_{1}+\left\|\widetilde{\rho}_1-\wt{\rho}_{1,s}\right\|_{1}+\left\|\wt{\rho}_{1,s}-\bar{\rho}_1\right\|_{1}.
\end{equation}

The first term contains the Lindblad simulation error, the second term contains the error from the numerical integration, and the third term contains the error from Trotter splitting. According to Lemma <ref> and Lemma <ref> (<ref>),

we first bound the Lindblad simulation error:

\begin{equation}\label{eqn:Lindbald_simulation_error}
\begin{aligned}
\leq &\left\|\widetilde{\rho}_1-\exp(\mathcal{L}_H\tau)\exp(\mathcal{L}_K\tau)\rho(0)\right\|_{1}+\left\|\exp((\mathcal{L}_H+\mathcal{L}_K)\tau)\rho(0)-\exp(\mathcal{L}_H\tau)\exp(\mathcal{L}_K\tau)\rho(0)\right\|_{1}\\
\end{aligned}
\end{equation}

where we use <ref> in the first equality and <ref> in the last equality.

Next, in order to bound the error of the numerical integration, we observe that

\[e^{iH\tau}\widetilde{\rho}_1 e^{-iH\tau}=\Tr_a \left(e^{-i\sqrt{\tau}\wt{K}} \ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}}\right)\,,\]

\[
e^{iH\tau}\wt{\rho}_{1,s}e^{-iH\tau}=\Tr_a \left(e^{-i\sqrt{\tau}\wt{K}_s} \ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}_s}\right) \,.
\]

\[
\widetilde{\rho}_{1,a}(t)=e^{-it\wt{K}} \ket{0}\bra{0}\otimes\rho(0) e^{it\wt{K}},\quad \widetilde{\rho}_{1,s,a}(t)=e^{-it\wt{K}_s} \ket{0}\bra{0}\otimes\rho(0) e^{it\wt{K}_s}\,.
\]

\[
\left\|\widetilde{\rho}_{1,s,a}(t)-\widetilde{\rho}_{1,s,a}(0)\right\|_1=\mathcal{O}(\|\wt{K}_s\|t)\,,
\]

we have

\begin{equation}\label{eqn:numerical_integral_error}
\begin{aligned}
=&\left\|\mathrm{Tr}_a\left(\int^{\sqrt{\tau}}_0\exp\left(-i\wt{K}(\sqrt{\tau}-t)\right)[\wt{K}-\wt{K}_s,\widetilde{\rho}_{1,s,a}(t)]\exp\left(i\wt{K}(\sqrt{\tau}-t)\right)\ud t\right)\right\|_1\\
\leq& \left\|\mathrm{Tr}_a\left(\int^{\sqrt{\tau}}_0\exp\left(-i\wt{K}(\sqrt{\tau}-t)\right)[\wt{K}-\wt{K}_s,\widetilde{\rho}_{1,s,a}(t)-\widetilde{\rho}_{1,s,a}(0)]\exp\left(i\wt{K}(\sqrt{\tau}-t)\right)\ud t\right)\right\|_1\\
&+\left\|\mathrm{Tr}_a\left(\int^{\sqrt{\tau}}_0\exp\left(-i\wt{K}(\sqrt{\tau}-t)\right)[\wt{K}-\wt{K}_s,\widetilde{\rho}_{1,s,a}(0)]\exp\left(i\wt{K}(\sqrt{\tau}-t)\right)\ud t\right)\right\|_1\\
\leq& \left\|e^{-i\sqrt{\tau}\wt{K}} \ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}}-e^{-i\sqrt{\tau}\wt{K}_s} \ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}_s}\right\|_1\\
\end{aligned}\,,
\end{equation}

where we use the Duhamel principle in the second equality, $\mathrm{Tr}_a([\wt{K}-\wt{K}_s,\widetilde{\rho}_{1,s,a}(0)])=0$ in the second inequality, and $\|K_{s}\|=\mathcal{O
}(\|A\|)$ in the last equality.

Finally, we bound the error of Trotter splitting. Using (<ref>) to expand each $\exp\left(-i\tilde{H}_l\sqrt{\tau}/2\right)$ up to $N=3$, we obtain that:

\begin{equation}\label{eqn:trotter_error}
\begin{aligned}
\left\|\wt{\rho}_{1,s}-\bar{\rho}_1\right\|_{1}=&\mathcal{O}\left(\left(\sqrt{\tau}\sum_{l}\|H_l\|\right)^4\right)=\mathcal{O}\left(\tau^2\|A\|^4\left(\tau_s\sum_{l}|f(s_l)|\right)^4\right)=\mathcal{O}\left(\tau^2\|A\|^4\right)
\end{aligned}
\end{equation}

where we use <ref> (<ref>) with $N=1$ in the third equality. Here all $\tau^{3/2}$ terms disappear because of the partial trace. Plugging (<ref>), (<ref>), and (<ref>) into (<ref>), we prove (<ref>).

To bound $\|K-K_{s}\|$ in (<ref>), we need to use <ref>. We first give its proof there.

We separate the error into two parts:

\begin{equation}\label{eqn:error_separation}
\|K-K_{s}\|\leq \|K-K_\infty\|+\left\|K_\infty-K_{s}\right\|\,,
\end{equation}

where $K_\infty=\sum^\infty_{l=-\infty}f(s_l)e^{iH s_l}Ae^{-iHs_l} \tau_s$. For a given $N$, using <ref> (<ref>), the second part of error can be bounded as

\begin{equation}\label{eqn:dis_infty_M_s}
\left\|K_\infty-K_{s}\right\|=\mathcal{O}\left( C_NS_\omega\Delta^{-N}\|A\|\sum_{|l|\geq M_s} (1+\abs{s_l})^{-N}\tau_s\right)=\mathcal{O}\left( C_NS_\omega\Delta^{-N}\|A\|\int^\infty_{S_s}(1+s)^{-N}\ud s\right)= \mathcal{O}\left(\epsilon'\right).
\end{equation}

To bound the discretization error $\|K-K_\infty\|$, first define

\[
\hat{f}_{s,\infty}(\omega)=\sum^\infty_{l=-\infty}f(s_l)\exp(i\omega s_l) \tau_s.
\]

Use the fact that $\abs{\lambda_i-\lambda_j}\le 2\norm{H}$, we have

\[
\begin{split}
\|K-K_\infty\|\leq & \norm{A}\sup_{i,j}\abs{\hat{f}(\lambda_i-\lambda_j)-\hat{f}_s(\lambda_i-\lambda_j)}\\
\|A\|\sup_{|\omega|\leq 2\|H\|}\left|\int^\infty_{-\infty} f(s)\exp(i\omega s)\ud s-\sum^\infty_{l=-\infty}f(s_l)\exp(i\omega s_l) \tau_s\right|
\end{split}
\]

Fixing an arbitrary $\omega\in [-2\|H\|,2\|H\|]$, we define $h(z)=f(z)\exp(i\omega z)$. According to <ref> (<ref>), we obtain that $h(z)$ is an analytic function in $z$ and

\[
\sup_{|\omega|\leq 2\|H\|,\ |y|\leq 1}\int^\infty_{-\infty}\left|h(x+iy)\right|\ud x=\mathcal{O}\left(C_NS_\omega\Delta^{-N}\exp(2\|H\|+S_\omega)\right)
\]

We then use <ref> to obtain

\[
\left|\int^\infty_{-\infty} h(s)\ud s-\sum^\infty_{l=-\infty}h(s_l)\tau_s\right|=\mathcal{O}\left(C_NS_\omega\Delta^{-N}\exp(2\|H\|+S_\omega)\exp(-2\pi/\tau_s)\right)\,.
\]

This gives

\begin{equation}\label{eqn:dis_infty}
\|K-K_\infty\|=\mathcal{O}\left(C_NS_\omega\Delta^{-N}\|A\|\exp(2\|H\|+S_\omega)\exp(-2\pi/\tau_s)\right)=\mathcal{O}\left(\epsilon'\right)\,.
\end{equation}

where we use $\tau_s=\mathcal{O}\left(\left(\|H\|+S_\omega+\log\left(\frac{C_NS_\omega\norm{A}}{\Delta\epsilon'}\right)\right)^{-1}\right)$ in the last inequality.
Plugging (<ref>) and (<ref>) into (<ref>), we obtain $\|K-K_{s}\|=\mathcal{O}(\epsilon')$.

Now, we are ready to prove <ref>:

According to <ref> (<ref>) and <ref>, we set $\epsilon'=\frac{\epsilon}{\|A\|T}$ in $\eqref{eqn:error_quadrature}$ to obtain that

\[
\left\|\rho(T)-\rho_{M_t}\right\|_{1}=\frac{\epsilon}{2}+\mathcal{O}\left((\|A\|^4+\|A\|^2\|H\|)T\tau\right)\leq \epsilon\,,
\]

where we use $\tau=\mathcal{O}\left(\frac{\epsilon}{(\|A\|^4+\|A\|^2\|H\|)T}\right)$ in the inequality.

Next, we calculate the total Hamiltonian simulation time for $H$. According to (<ref>), each implementation of $W(\tau)$ needs to simulate the system Hamiltonian for time $\Theta(S_s)$. Thus, the total Hamiltonian simulation time:

\[
T_{H,\mathrm{total}}=\text{number of steps}\times (\tau+\Theta(S_s))=M_t(\tau+\Theta(S_s))=\Theta(T+TS_s
\tau^{-1})=\widetilde{\Theta}(T^{2+o(1)}\epsilon^{-1-o(1)})\,.
\]

Finally, to calculate the number of controlled-$A$ evolution gates, we notice that each trotter splitting step in $W(\tau)$ needs to implement $\Theta(1)$ controlled-$A$ evolution gates. This implies

\[
N_{A,\mathrm{gate}}=\text{number of steps}\times \Theta(S_s\tau^{-1}_s)=\widetilde{\Theta}(T^{2+o(1)}\epsilon^{-1-o(1)})\,.
\]

§.§ Proof of Lg and Lg

To prove <ref>, we first show the following proposition:
Under the conditions of <ref>, we further assume $r=\Omega(1)$, then

\begin{equation}\label{eqn:acc_discretization_error}
\begin{aligned}
\left\|e^{iHS_s}\rho^a_{M_t}e^{-iHS_s}-\rho_{M_t}\right\|_{1}=\mathcal{O}\left(\|K-K_{s}\|\|A\|T\right)+\mathcal{O}\left(\|A\|^4T\tau/r^2\right).
\end{aligned}
\end{equation}

For simplicity, we ignore the effect of $e^{iHS_s}$. Because the Lindblad dynamics is contractive in trace distance  according to <ref>, to prove (<ref>), it suffices to prove

\begin{equation}\label{eqn:one_step_error_acc}
\begin{aligned}
\left\|\rho^a_1-\rho_1\right\|_{1}=\mathcal{O}\left(\|K-K_{s}\|\|A\|\tau\right)+\mathcal{O}\left(\|A\|^4\tau^2/r^2\right).
\end{aligned}
\end{equation}

\[
\wt{\rho}_{1,s}=e^{-iH\tau}\left(\Tr_a e^{-i\sqrt{\tau}\wt{K}_s} \ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}_s}\right)e^{iH\tau}\,.\]

\begin{equation}\label{eqn:error_separation_2_acc}
\left\|\rho^a_1-\rho_1\right\|_{1}\leq \left\|\rho^a_1-\wt{\rho}_{1,s}\right\|_{1}+\left\|\wt{\rho}_{1,s}-\rho_1\right\|_{1}\,,
\end{equation}

where the first term contains the error from the numerical integration, and the third term contains the error from Trotter splitting. Similar to <ref> in the proof of Theorem <ref>, for the first term, we can bound is by:

\begin{equation}\label{eqn:numerical_integral_error_acc}
\left\|\rho^a_1-\wt{\rho}_{1,s}\right\|_{1}\leq \mathcal{O}(\|K-K_s\|\|A\|\tau)\,.
\end{equation}

To bound the second term, we notice

\begin{equation}\label{eqn:error_first_acc}
\begin{aligned}
\leq &\left\|\Tr_a \left(e^{-i\wt{K}_s\sqrt{\tau}} \ket{0}\bra{0}\otimes\rho(0) e^{i\wt{K}_s\sqrt{\tau}}-[W(\sqrt{\tau}/r)]^r \left(\ket{0}\bra{0}\otimes\rho(0)\right)  [(W(\sqrt{\tau}/r))^{\dag}]^r\right)\right\|_1.
\end{aligned}
\end{equation}

Using the second-order Trotter splitting formula and <ref> with $N=3$, we can rewrite $W(\sqrt{\tau}/r)$ as

\begin{equation}\label{eqn:W_eta_r}
\end{equation}

Here the derivation of the expression

\[
\]

is similar to that in <ref>. Moreover,

\[
\|E_2(\tau,r)\|=\mathcal{O}\left(\frac{\tau^2}{r^4}\left(\sum_{l}\|H_l\|\right)^4\right)=\mathcal{O}\left(\frac{\tau^2}{r^4}\|A\|^4\left(\tau_s\sum_{l}|f(s_l)|\right)^4\right)=\mathcal{O}\left(\frac{\|A\|^4\tau^2}{r^4}\right)\,.
\]

Here, in order to obtain the correct leading order error, we need to expand $W(\sqrt{\tau}/r)$ up to $\mathcal{O}(\tau^2)$.

Plugging (<ref>) into (<ref>), we have

\[
\begin{aligned}
\leq &\sum^r_{k=1}\left\|  \Tr_a \left(\frac{\tau^{3/2}}{r^3}e^{-i(k-1)\sqrt{\tau}\wt{K}_s/r}E_1e^{-i(r-k)\sqrt{\tau}\wt{K}_s/r}\ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}_s}\right)\right\|_1\\
+&\sum^r_{k=1}\left\|  \Tr_a \left(\frac{\tau^{3/2}}{r^3}e^{-i\sqrt{\tau}\wt{K}_s}\ket{0}\bra{0}\otimes\rho(0) e^{i(k-1)\sqrt{\tau}\wt{K}_s/r}E^\dagger_1e^{i(r-k)\sqrt{\tau}\wt{K}_s/r}\right)\right\|_1\\
\end{aligned}
\]

Here, we use the relationship $r=\Omega(1)$ and $\tau=\mathcal{O}(\|A\|^{-2})$ to incorporate all terms of the form $\left(\frac{|\tau|^{p/2}\|A\|^p}{r^p}\right)$ (for $p\ge 4$) into the asymptotic notation  $\mathcal{O}\left(\frac{\tau^2\|A\|^4}{r^4}\right)$.

Next, since $\Tr_a \left(E_1\ket{0}\bra{0}\otimes\rho(0) \right)=0$, we can expand $e^{i\sqrt{\tau}\wt{K}_s}$ to first order and obtain the refined estimate

\[
\left\|\Tr_a \left(\frac{\tau^{3/2}}{r^3}e^{-i(k-1)\sqrt{\tau}\wt{K}_s/r}E_1e^{-i(r-k)\sqrt{\tau}\wt{K}_s/r}\ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}_s}\right)\right\|_1=\mathcal{O}\left(\frac{\tau^2\|A\|^4}{r^3}\right)
\]

for all $1\leq k\leq r$. Here we use $\norm{E_1}=\Or(\norm{A}^3)$, and $\norm{K}=\Or(\norm{A})$.

\begin{equation}
\sum^r_{k=1}\left\|  \Tr_a \left(\frac{\tau^{3/2}}{r^3}e^{-i(k-1)\sqrt{\tau}\wt{K}_s/r}E_1e^{-i(r-k)\sqrt{\tau}\wt{K}_s/r}\ket{0}\bra{0}\otimes\rho(0) e^{i\sqrt{\tau}\wt{K}_s}\right)\right\|_1=\mathcal{O}\left(\frac{\tau^2\|A\|^4}{r^2}\right).
\end{equation}

Similarly using $\Tr_a \left(\ket{0}\bra{0}\otimes\rho(0)E^\dagger_1 \right)=0$, we have

\begin{equation}
\sum^r_{k=1}\left\|  \Tr_a \left(\frac{\tau^{3/2}}{r^3}e^{-i\sqrt{\tau}\wt{K}_s}\ket{0}\bra{0}\otimes\rho(0) e^{i(k-1)\sqrt{\tau}\wt{K}_s/r}E^\dagger_1e^{i(r-k)\sqrt{\tau}\wt{K}_s/r}\right)\right\|_1\\
\end{equation}

This gives

\[
\left\|\wt{\rho}_{1,s}-\rho_1\right\|_{1}=\mathcal{O}\left(\frac{\tau^2\|A\|^4}{r^2}\right)\,.
\]

Plugging this equality and (<ref>) into (<ref>), we prove (<ref>).

Now, we are ready to prove <ref>:

We note that <ref> also holds for this case. Thus, we apply <ref> (<ref>) and <ref> by setting $\epsilon'=\frac{\epsilon}{\|A\|T}$ in $\eqref{eqn:error_quadrature}$ to obtain that

\[
\left\|\rho(T)-\rho_{M_t}\right\|_{1}=\frac{\epsilon}{2}+\mathcal{O}\left(\|A\|^4T\tau/r^2\right)\leq \epsilon\,,
\]

where we use $r=\Theta(\|A\|T^{1/2}\epsilon^{-1/2})$ and $\tau=\mathcal{O}(\|A\|^{-2})$ in the inequality.

Next, we calculate the total Hamiltonian simulation time for $H$. According to (<ref>), the implementation of each $W(\tau/r)$ needs to simulate the system Hamiltonian for time $\Theta(S_s)$. Thus, the total Hamiltonian simulation time:

\[
T_{H,\mathrm{total}}=\text{number of steps}\times (\tau+r\Theta(S_s))=\Theta(T+rTS_s
\tau^{-1})=\widetilde{\Theta}(T^{3/2+o(1)}\epsilon^{-1/2-o(1)})\,.
\]

Finally, each step within the Trotter-splitting process of $W(\tau/r)$ requires the implementation of $\Theta(1)$ controlled-$A$ evolution gates.
This implies

\[
N_{A,\mathrm{gate}}=\text{number of steps}\times \Theta(S_s\tau^{-1}_sr)=\widetilde{\Theta}(T^{3/2+o(1)}\epsilon^{-1/2-o(1)})\,.
\]

We note that <ref> also holds for this case. In addition, when applying $p$-th order Trotter scheme, we obtain

\begin{equation}\label{eqn:acc_discretization_error_2}
\begin{aligned}
\left\|e^{iHS_s}\rho^a_{M_t}e^{-iHS_s}-\rho_{M_t}\right\|_{1}=\mathcal{O}\left(\|K-K_{s}\|\|A\|T\right)+\mathcal{O}\left(\|A\|^{p+2}T\tau^{p/2}/r^p\right)\,.
\end{aligned}
\end{equation}

Then, we apply (<ref>) and <ref> by setting $\epsilon'=\frac{\epsilon}{\|A\|T}$ in $\eqref{eqn:error_quadrature}$ to obtain that

\[
\left\|\rho(T)-\rho_{M_t}\right\|_{1}=\frac{\epsilon}{2}+\mathcal{O}\left(\|A\|^{p+2}T\tau^{p/2}/r^p\right)\leq \epsilon\,,
\]

where we use $r=\Theta(\|A\|^{2/p}T^{1/p}\epsilon^{-1/p})$ and $\tau=\mathcal{O}(\|A\|^{-2})$ in the inequality.

Similar to the previous proof, the total Hamiltonian simulation time is

\[
T_{H,\mathrm{total}}=\text{number of steps}\times (\tau+r\Theta(S_s))=\Theta(T+rTS_s
\tau^{-1})=\widetilde{\Theta}(T^{1+1/p}\epsilon^{-1/p})
\]

and the total number of controlled-$A$ gates is

\[
N_{A,\mathrm{gate}}=\text{number of steps}\times \Theta(S_s\tau^{-1}_sr)=\widetilde{\Theta}(T^{1+1/p}\epsilon^{-1/p})\,.
\]

This concludes the proof.

§ CONVERGENCE UNDER ETH TYPE ANSATZ

In this section, we will investigate the convergence of our algorithm to the ground state based on appropriate assumptions. Specifically, we focus on the continuous-time Lindblad dynamics (<ref>).

Firstly, we study the convergence property. Define $A_{i,j}=\left\langle \psi_i\right|A\ket{\psi_j}$ and $\hat{f}_{i,j}=\hat{f}(\lambda_i-\lambda_j)$.
In order to study the convergence of the continuous-time Lindblad dynamics (<ref>), we need to make some independent assumptions for $A$.

* (Random matrix elements) Assume for any $t\geq0$, $A$ is independently drawn from  random probability distribution $\Xi_A$ on the set of Hermitian matrices such that $A_{i,j}$ are independent and $\mathbb{E}(A_{i,j})=0$ when $i\neq j$. Denote $\sigma_{i,j}=\mathbb{E}(|A_{i,j}|^2)>0$.

* (Support of filter) $S_\omega>\min^{N-2}_{j=0}(\lambda_{j+1}-\lambda_j)$.

* (Diagonal initial state) $\rho(0)$ is a diagonal matrix in the basis of $\{\ket{\psi_i}\}^{N-1}_{i=0}$.

Strictly speaking, in the above assumption, $A$ should be a function of time and denoted as $A_t$. However, for the sake of simplicity and consistency with other notations, we will omit the subscript $t$.

Based on <ref>, we demonstrate the convergence in the following theorem:

* Under Assumption <ref>, $\rho^\star=\ket{\psi_0}\bra{\psi_0}$ is the unique fixed point of (<ref>) in the expectation sense. In particular,

\[
\lim_{t\rightarrow\infty} \mathbb{E}(\rho(t))=\ket{\psi_0}\bra{\psi_0}\,.
\]

We put the proof of Theorem <ref> in Appendix <ref>. It is important to highlight that while the expected operator $\mathbb{E}(\rho(t))$ shares a conceptual connection with $\rho(t)$, they are not identical. Nonetheless, this connection provides valuable insights into ergodicity, offering a promising perspective. For more details, please refer to Section <ref>.

Next, we study the convergence speed (or mixing time) of the continuous-time Lindblad dynamics (<ref>). According to the proof of Theorem <ref> in Appendix <ref>, we define $\textbf{p}(t)=(p_0(t),p_1(t),\dots,p_{N-1}(t))^\top$,
then the solution of (<ref>) satisfies

\begin{equation}
\frac{\ud\textbf{p}(t)}{\ud t}=\mathbf{T}\mathbf{p}(t), \quad \mathbb{E}(\rho(t))=\sum^{N-1}_{i=0}p_i(t)\ket{\psi_i}\left\langle \psi_i\right|.
\end{equation}

Here the transition matrix elements are

\begin{equation}\label{eqn:Transition_matrix}
\mathbf{T}_{j,i}=\hat{f}^2_{j,i}\sigma_{j,i},\quad i\neq j,\quad \mbox{and}\quad  \mathbf{T}_{i,i}=-\sum^{N-1}_{j\neq i}\hat{f}^2_{j,i}\sigma_{j,i}.
\end{equation}

We are now prepared to demonstrate that the solution of the continuous-time Lindblad dynamics (<ref>) rapidly converges to low-energy states, as proven by the following theorem:
Together with Assumption <ref>, we also assume that there exists a decreasing sequence $\{R_l\}^{L}_{l=1}$ with $L=\mathcal{O}(\mathrm{poly}(n))$ such that

* $R_1=N-1$, and $R_L=\mathcal{O}(\mathrm{poly}(n))$.

* For each $1\leq l\leq L-1$

\begin{equation}\label{eqn:weight_out}
\sum^{R_{l+1}}_{i=0}\hat{f}^2_{i,j}\sigma_{i,j}=\Omega(1/\mathrm{poly}(n)),\quad \text{for all}\  j\in (R_{l+1},R_l].
\end{equation}

Then, there exists $T^\star=\mathcal{O}(\mathrm{poly}(n))$ such that

\begin{equation}\label{eqn:convergence}
\mathbb{E}\left(\sum^{R_L}_{i=0}\bra{\psi_i}\rho(t)\ket{\psi_i}\right)=\sum^{R_L}_{i=0}p_i(t)=\Omega(1),\quad \text{for all}\  t>T^\star\,.
\end{equation}

for any initial condition $\textbf{p}(0)$ such that $p_i(0)\geq 0$ and $\sum^N_{i=1}p_i(0)=1$.

We would like to emphasize that condition (<ref>) is crucial in ensuring that the weight in $p_{R_{l+1}+1\leq i\leq R_{l}}$ can be efficiently transferred to $p_{i\leq R_{l}+1}$ within $\mathcal{O}(\mathrm{poly}(n))$ time. Since this transfer only needs to occur $L=\mathcal{O}(\mathrm{poly}(n))$ times, the total time required for the weight translation from $p_{R_{L}+1\leq i}$ to $p_{i\leq R_{L}}$ remains $\mathcal{O}(\mathrm{poly}(n))$. To provide a clear visual representation of the transition matrix's structure, we have included its graph in (<ref>).

The theorem above implies that, given appropriate assumptions, the continuous-time Lindblad dynamics (<ref>) can rapidly drive the quantum state towards low-energy eigenstates, resulting in an increased overlap with them. Now, if we assume that the mixing time of the sub-transition matrix $\textbf{T}_{(1:R_L,1:R_L)}\in \mathbb{R}^{R_L\times R_L}$ is $t^{\mathrm{sub}}_{\mathrm{mix}}=\mathcal{O}(1/\mathrm{poly}(n))$, then the mixing time for preparing the ground state is $t_{\mathrm{mix}}=\mathcal{O}(\mathrm{poly}(n))$. In other words, when $t>t_{\mathrm{mix}}$, the overlap with the ground state $\mathbb{E}\left(\bra{\psi_0}\rho(t)\ket{\psi_0}\right)=a_0(s)=\Omega(1)$.

Structure of the transition matrix $T$ in the eigenbasis of the Hamiltonian.

§.§ Proof of Lg

Since $\hat{f}(\omega)=0$ for $\omega\geq 0$, we have $\mc{L}_K[\ket{\psi_0}\bra{\psi_0}]=0$. Plugging this into (<ref>), we prove that
$\rho^\star=\ket{\psi_0}\bra{\psi_0}$ is a fixed point of (<ref>).

Now, we assume $A,\rho(0)$ satisfy Assumption <ref>. To track the evolution of $\mathbb{E}(\rho(t))$, we first observe,

\[
K\ket{\psi_i}=\sum^{N-1}_{j=0} \hat{f}_{j,i}A_{j,i}\ket{\psi_j},\quad K^\dagger\ket{\psi_i}=\sum^{N-1}_{j=0} \hat{f}_{i,j}A_{j,i}\ket{\psi_j}\,,
\]

where we use $A$ as a Hermitian matrix in the second equality. Taking the expectation on the randomness of $A$, we obtain

\[
\mathbb{E}\left(K\ket{\psi_i}\left\langle \psi_i\right|K^\dagger\right)=\sum^{N-1}_{j=0} \hat{f}^2_{j,i}\sigma_{j,i}\ket{\psi_j}\left\langle \psi_j\right|\,,
\]

\[
\begin{aligned}
\mathbb{E}\left(K^\dagger K\ket{\psi_i}\left\langle \psi_i\right|\right)&=\mathbb{E}\left(\sum^{N-1}_{j=0}\hat{f}_{j,i}A_{j,i}K^\dagger \ket{\psi_j}\left\langle \psi_i\right|\right)\\
&=\mathbb{E}\left(\sum^{N-1}_{j,k=0}\hat{f}_{j,i}A_{j,i}f_{j,k}A_{k,j}\ket{\psi_k}\left\langle \psi_i\right|\right)\\
&=\sum^{N-1}_{j=0} \hat{f}^2_{j,i}\sigma_{j,i}\ket{\psi_i}\left\langle \psi_i\right|\,.
\end{aligned}
\]

These two equalities imply that

\begin{equation}\label{eqn:lindbload_psi}
\mathbb{E}\left(\mathcal{L}_K\left(\ket{\psi_i}\left\langle \psi_i\right|\right)\right)=\sum^{N-1}_{j\neq i} \hat{f}^2_{j,i}\sigma_{j,i}\left(\ket{\psi_j}\left\langle \psi_j\right|-\ket{\psi_i}\left\langle \psi_i\right|\right)\,.
\end{equation}

Since $\rho(0)$ is a diagonal matrix in the basis of $\{\ket{\psi_i}\}^{N-1}_{i=0}$, plugging (<ref>) into (<ref>) and taking expectation on both sides, we find that $\mathbb{E}(\rho(t))$ is always a diagonal matrix in the basis of $\{\ket{\psi_i}\}^{N-1}_{i=0}$ and

\[
\frac{\ud\mathbb{E}(\rho(t))}{\ud t}=\mathbb{E}\left(\mathcal{L}_K\left(\rho(t)\right)\right)\,,
\]

where we use $[H,\mathbb{E}(\rho(t))]=0$ for any $t$. According to (<ref>), we can rewrite the above equation as

\[
\mathbb{E}(\rho(t))=\sum^{N-1}_{i=0}p_i(t)\ket{\psi_i}\left\langle \psi_i\right|\,
\]

where $p_i(t)$ solves

\begin{equation}\label{eqn:evolution_a_n}
\frac{\ud p_i(t)}{\ud t}=\sum^{N-1}_{j\neq i} \hat{f}^2_{i,j}\sigma_{i,j}p_{j}(t)-\sum^{N-1}_{j\neq i} \hat{f}^2_{j,i}\sigma_{j,i}p_i(t)\,.
\end{equation}

In particular,

\[
\frac{\ud a_0(t)}{\ud t}=\sum^{N-1}_{j=1} \hat{f}^2_{0,j}\sigma_{0,j}p_{j}(t)\,,
\]

which implies $p_0(t)=1,\ p_{i\geq 1}(t)=0$ is the unique fixed point. Since $\sigma_{i,j}>0$ for all $i,j$, we conclude that $\ket{\psi_0}\bra{\psi_0}$ is the unique fixed point of the Lindblad dynamics (<ref>), and the solution of (<ref>) converges to $\ket{\psi_0}\bra{\psi_0}$ as $t\rightarrow\infty$.

Next, we prove <ref>.

To show (<ref>), it suffices to prove that for any $1\leq l\leq L-1$, there exists $T_l=\mathcal{O}(\mathrm{poly}(n))$ such that

\begin{equation}\label{eqn:iteration}
\sum^{N-1}_{i=R_{l+1}+1} p_i(t)\leq \frac{1}{2}-\frac{1}{l+3},\quad \text{for all}\  t>T_l\,.
\end{equation}

Since $\hat{f}(x)=0$ when $x\geq 0$, we first notice that $\textbf{T}$ is an upper triangular matrix, which implies the weight can only move from high energy states to low energy states. Then, for $l=1$, using (<ref>), we obtain

\[
\frac{\ud \sum^{N-1}_{i=R_{2}+1} p_i(t)}{\ud t}\leq -\left(\min_{j\in (R_{2},R_1]}\sum^{R_{2}}_{i=0}\hat{f}^2_{i,j}\sigma_{i,j}\right)\left(\sum^{N-1}_{i=R_{2}+1} p_i(t)\right)=\mathcal{O}\left(-\frac{1}{\mathrm{poly}(n)}\left(\sum^{N-1}_{i=R_{2}+1} p_i(t)\right)\right)\,,
\]

Since $\sum^{N-1}_{i=R_{2}+1} p_i(0)\leq 1$, there exists $T_1=\mathcal{O}(\mathrm{poly}(n))$ such that

\begin{equation}\label{eqn:first_iteration}
\sum^{N-1}_{i=R_{2}+1} p_i(t)<\frac{1}{4},\quad \text{for all}\  t>T_1\,.
\end{equation}

Next, for $l=2$, according (<ref>), we obtain

\[
\frac{\ud \sum^{N-1}_{i=R_{3}+1} p_i(t)}{\ud t}\leq -\left(\min_{j\in (R_{3},R_2]}\sum^{R_{3}}_{i=0}\hat{f}^2_{i,j}\sigma_{i,j}\right)\left(\sum^{R_2}_{i=R_{3}+1} p_i(t)\right)=\mathcal{O}\left(-\frac{1}{\mathrm{poly}(n)}\left(\sum^{R_2}_{i=R_{3}+1} p_i(t)\right)\right)\,,
\]

where we use $\min_{j\in (R_{3},R_2]}\sum^{R_{3}}_{i=0}\hat{f}^2_{i,j}\sigma_{i,j}=\Omega(\frac{1}{\mathrm{poly}(n)})$ in the inequality. This implies that when $T>T_1$, we have

\[
\frac{\ud \sum^{N-1}_{i=R_{3}+1} p_i(t)}{\ud t}\leq \left\{
\begin{aligned}
&\mathcal{O}\left(-\frac{1}{20}\frac{1}{\mathrm{poly}(n)}\right),\quad \text{if}\quad \sum^{R_2}_{i=R_{3}+1} p_i(t)>\frac{1}{20}\\
&0,\quad \text{otherwise}\\
\end{aligned}
\right.\,.
\]

Because $\sum^{R_2}_{i=R_{3}+1} p_i(t)\leq \sum^{N-1}_{i=R_{3}+1} p_i(t)$, there exists $T_2=\max\left\{T_1,\mathcal{O}(\mathrm{poly}(n))\right\}$ such that

\[
\sum^{R_2}_{i=R_{3}+1} p_i(t)<\frac{1}{1+3}-\frac{1}{2+3}=\frac{1}{20},\quad \text{for all}\  t>T_2\,.
\]

Combining this with (<ref>), we can obtain that

\[
\sum^{N-1}_{i=R_{3}+1} p_i(t)<\frac{3}{10}=\frac{1}{2}-\frac{1}{2+3},\quad \text{for all}\  t>T_2\,.
\]

We can continue this process to any $l\leq L-1$. Specifically, for any $l$, we obtain that

\[
\frac{\ud \sum^{N-1}_{i=R_{l+1}+1} p_i(t)}{\ud t}\leq \left\{
\begin{aligned}
&\mathcal{O}\left(-\frac{1}{(l-1+3)(l+3)}\frac{1}{\mathrm{poly}(n)}\right),\quad \text{if}\quad \sum^{R_l}_{i=R_{l+1}+1} p_i(t)>\frac{1}{(l-1+3)(l+3)}\\
&0,\quad \text{otherwise}\\
\end{aligned}
\right.\,.
\]

Thus, there exists $T_l=\max\left\{T_{l-1},\mathcal{O}(\mathrm{poly}(n))\right\}$ such that

\[
\sum^{R_l}_{i=R_{l+1}+1} p_i(t)<\frac{1}{l-1+3}-\frac{1}{l+3}=\frac{1}{20}
\]

\[
\sum^{N-1}_{i=R_{l+1}+1} p_i(t)<\frac{1}{2}-\frac{1}{l+3}
\]

for any $t>T_l$. This proves (<ref>).

§.§ Concentration around deterministic dynamics

In this section, we demonstrate that, when the coupling operator $A$ satisfies <ref>, our continuous Lindblad dynamics simulation in <ref> approximates $\mathbb{E}\left(\exp(-\mathcal{L}_H S_s)[\rho(T)]\right)$ when $\tau\rightarrow0$, where $\rho(T)$ the solution of the modified Lindblad dynamics (<ref>).

For simplicity, we omit the phase shift $\exp(-\mathcal{L}_H S_s)$ and assume our simulation scheme can fit into the following classical setting: Given a stopping time $T>0$ and a small time step $\tau>0$ such that $M=T/\tau\in\mathbb{N}$, we approximate $\rho(n\tau)$ using $\rho_n$, which is defined by

\begin{equation}\label{eqn:simulation_update}
\rho_n=\mathcal{F}(\tau,A_n)\rho_{n-1},\quad \text{where}\ \rho_0=\rho(0)\,.
\end{equation}

Here, $\mathcal{F}(\tau,A_n)$ is a linear operator that depends on the random operator $A_n$ and $\{A_n\}^M_{n=1}$ are independently drawn from a probability distribution $\Xi_{A}$.

Next, we introduce the Frobenius norm of a matrix, which is defined as

\[
\|A\|_F:=\sqrt{\mathrm{Tr}\left(A^\dagger A\right)}\,.
\]

In addition, given a superoperator $\mathcal{L}$ that acts on matrices, the induced $F$-norm is

\[
\left\|\mathcal{L}\right\|_F:=\sup_{\|A\|_F\leq 1}\left\|\mathcal{L}(A)\right\|_F\,.
\]

Now, we are ready to introduce the approximation result, which is summarized in the following theorem:

Assuming our scheme possesses at least first-order accuracy, meaning that there exists a constant $C_\mathcal{F}$ such that

\[
\left\|\mathcal{F}(\tau,A)-\exp(\mathcal{L}_{A}\tau)\right\|_1\leq C_\mathcal{F}\tau^2
\]

almost surely for $A$ drawn from $\Xi_A$. Define $C_\Xi=\sup_{A\in\mathrm{supp}(\Xi_A)}\max\left\{\|\mathcal{L}_A\|_1,\|\mathcal{L}_A\|_F\right\}$.
If $\tau=\mathcal{O}\left(\min\left\{\frac{1}{T},\frac{1}{C_\Xi}\right\}\right)$, we have

\begin{equation}\label{eqn:half_order_converge}
\mathbb{E}\left\|\rho_M-\mathbb{E}(\rho(T))\right\|_F=\mathcal{O}\left(\exp(C_\Xi T)C_\Xi\sqrt{T\tau}+C_{\mathcal{F}}T\tau\right)
\end{equation}

The relationship between the $F$-norm and the trace norm in (<ref>) can be expressed as follows:

\[
\frac{1}{\sqrt{N}}\left\|\rho_M-\mathbb{E}(\rho(T))\right\|_1\leq \left\|\rho_M-\mathbb{E}(\rho(T))\right\|_F\leq \left\|\rho_M-\mathbb{E}(\rho(T))\right\|_1\,.
\]

It is important to highlight that these inequalities are sharp. Therefore, in the worst-case scenario, the upper bound of $\left\|\rho_M-\mathbb{E}(\rho(T))\right\|_1$ depends exponentially on the number of qubits $n$. It remains an intriguing question whether there exists potential for further improvement in the bound of $\left\|\rho_M-\mathbb{E}(\rho(T))\right\|_1$.

\[
\widetilde{\rho}_n=\left(1+\mathbb{E}(\mathcal{L}_{A_n})\tau\right)\widetilde{\rho}_{n-1},\quad \text{where}\ \widetilde{\rho}_0=\rho(0)\,.
\]

Here $\mathcal{L}_{A_n}$ is defined in (<ref>), and the subscript $A_n$ indicates its dependence on the random matrix $A_n$. We notice that

\[
\frac{\mathrm{d}\mathbb{E}(\rho(t))}{\mathrm{d}t}=\mathbb{E}(\mathcal{L}_A)\mathbb{E}(\rho(t))\,.
\]

Because $\exp(\mathbb{E}(\mathcal{L}_A) t)$ is a completely positive trace-preserving map for any $t\geq 0$,

\begin{equation}\label{eqn:d_small_1}
\left\|\mathbb{E}(\rho(T))-\widetilde{\rho}_M\right\|_F\leq \left\|\mathbb{E}(\rho(T))-\widetilde{\rho}_M\right\|_1=\mathcal{O}\left(C_\Xi^2T\tau\right)\,.
\end{equation}

Next, we rewrite (<ref>) as

\begin{equation}\label{eqn:simulation_update_new}
\begin{aligned}
\rho_n=&\left(1+\mathcal{L}_{A_n}\tau\right)\rho_{n-1}\\
\end{aligned}
\end{equation}

The second and third terms can be bounded by:

\begin{equation}\label{eqn:d_small_2}
\left\|\left(\mathcal{F}(\tau,A_n)-\exp(\mathcal{L}_{A_n}\tau)\right)\rho_{n-1}\right\|_1\leq C_\mathcal{F}\tau^2\,,
\end{equation}

\begin{equation}\label{eqn:d_small_3}
\left\|\left(\exp(\mathcal{L}_{A_n}\tau)-\left(1+\mathcal{L}_{A_n}\tau\right)\right)\rho_{n-1}\right\|_1=\mathcal{O}(C_\Xi^2\tau^2)\,.
\end{equation}

Plugging (<ref>), (<ref>) into (<ref>), when $\tau<\|\mathcal{L}_{A_n}\|_1^{-1}$, $(1+\mathcal{L}_{A_n}\tau)$ are completely positive trace-preserving maps, and we obtain

\begin{equation}\label{eqn:second_key}
\left\|\rho_M-\Pi^N_{n=1}\left(1+\mathcal{L}_{A_n}\tau\right)\rho_{0}\right\|_F\leq \left\|\rho_M-\Pi^N_{n=1}\left(1+\mathcal{L}_{A_n}\tau\right)\rho_{0}\right\|_1=\mathcal{O}((C_\mathcal{F}+C_\Xi^2)T\tau)\,.
\end{equation}

Finally, we borrow idea from <cit.> to bound $\mathbb{E}\left\|\Pi^N_{n=1}\left(1+\mathcal{L}_{A_n}\tau\right)\rho_{0}-\widetilde{\rho}_M\right\|^2_F$. First, we rewrite

\[
\begin{aligned}
\end{aligned}
\]

The key observation is that only (I) contains the random operator $\mathcal{L}_{A_N}$. Thus,

\[
\mathbb{E}\left(\mathrm{(I)}^\dagger\mathrm{(II)}\right)=\mathbb{E}\left(\mathbb{E}\left(\mathrm{(I)}^\dagger\mathrm{(II)}\middle|A_{1:n-1}\right)\right)=0\,.
\]

and $\mathbb{E}\left(\mathrm{(II)}^\dagger\mathrm{(I)}\right)=0$. These further imply

\[
\begin{aligned}
\leq \mathbb{E}\left\|\mathrm{(I)}\right\|^2_1+(1+M\tau)^2\mathbb{E}\left\|\left(\Pi^{N-1}_{n=1}\left(1+\mathcal{L}_{A_n}\tau\right)\rho_{0}-\Pi^{N-1}_{n=1}\left(1+\mathbb{E}\left(\mathcal{L}_{A_n}\right)\tau\right)\rho_{0}\right)\right\|^2_F\\
\le&4M^2\tau^2+\exp(2M\tau)\mathbb{E}\left\|\left(\Pi^{N-1}_{n=1}\left(1+\mathcal{L}_{A_n}\tau\right)\rho_{0}-\Pi^{N-1}_{n=1}\left(1+\mathbb{E}\left(\mathcal{L}_{A_n}\right)\tau\right)\rho_{0}\right)\right\|^2_F\,.
\end{aligned}
\]

Repeating the above calculation iteratively, we obtain

\[
\mathbb{E}\left\|\Pi^N_{n=1}\left(1+\mathcal{L}_{A_n}\tau\right)\rho_{0}-\widetilde{\rho}_M\right\|^2_F=\mathcal{O}\left(\exp(2C_\Xi T)C_\Xi^2T\tau\right)\,.
\]

\begin{equation}\label{eqn:d_small_4}
\mathbb{E}\left\|\rho_M-\widetilde{\rho}_M\right\|_F\leq \left(\mathbb{E}\left\|\Pi^N_{n=1}\left(1+\mathcal{L}_{A_n}\tau\right)\rho_{0}-\widetilde{\rho}_M\right\|^2_F\right)^{1/2}=\mathcal{O}\left(\exp(C_\Xi T)C_\Xi\sqrt{T\tau}\right)\,.
\end{equation}

Combining (<ref>), (<ref>), and (<ref>), we prove (<ref>).

§ DETAILS OF NUMERICAL EXPERIMENTS AND ADDITIONAL EXPERIMENTS

In our numerical simulations, we consider the following TFIM models with $L$ sites:

\begin{equation}\label{eqn:H_Ising}
H=-\left(\sum^{L-1}_{i=1} Z_{i}Z_{i+1}\right) -g\sum^L_{i=1} X_i\,,
\end{equation}

where $g$ is the coupling coefficient, $Z_i,X_i$ are Pauli operators for the $i$-th site and the dimension of $H$ is $2^L$. In the numerical test in Numerics results, we set $L=6$ and the coupling constant $g=1.2$. For the jump operator $K$, we utilize a local Hermitian operator $A = Z\otimes I_{2^{L-1}}$ and choose $\hat{f}(\omega)$ as follows:

\begin{equation}\label{eqn:f_omega}
\hat{f}(\omega) := \frac{1}{2}\left(\text{erf}\left(\frac{\omega+a}{\delta_a}\right)-\text{erf}\left(\frac{\omega+b}{\delta_b}\right)\right),
\end{equation}

where $\text{erf}(\omega) = \frac{2}{\sqrt{\pi}}\int_0^\omega e^{-x^2} \mathrm{d}x$ denotes the error function. The parameters $a$ and $\delta_a$ are chosen to be of the order $S_\omega$, while $b$ and $\delta_b$ are of the order $\Delta$. The inverse Fourier transform $f(s) = \int_{\mathbb{R}} f(\Omega)e^{-i\omega s}\mathrm{d}\omega$ is given by

\begin{equation}\label{eqn:F}
f(s) = \frac{e^{-\frac{\delta_a^2 s^2}{4}}e^{ias}-e^{-\frac{\delta_b^2 s^2}{4}}e^{ibs}}{2\pi is},.
\end{equation}

It is worth noting that $\lim_{t\to 0} f(s) = \frac{a-b}{2\pi}$ is well-defined, and $f(s)$ is a smooth complex function that is approximately supported on the interval $[-S_s, S_s]$, where $S_s = \Theta(\delta_b^{-1})$. Although (<ref>) does not strictly satisfy <ref>, as $\hat{f}(\omega)$ is approximately supported in $[-2a,0]$ and $\hat{f}(0)\approx 0$, our numerical tests have shown that this choice of $f$ still yields satisfactory results. The graph of $f$ can be found in Figure <ref>.

Illustration of $\hat{f}(\omega)$ and its Fourier transform $f(s)$ following <ref> and <ref>.

For simulating the continuous/discrete-time Lindblad dynamics, we always choose $a=2.5\|H\|$, $\delta_a=0.5\|H\|$, $b=\delta_b=\Delta$, $S_s=5/\delta_b$, and $\tau_s=\pi/(2a)$.

The code used for generating the numerical results is available on Github (<https://github.com/lin-lin/oneancillaground>).

§.§ TFIM-4 model

To simulate the Lindblad dynamics with TFIM-4 model, we choose $L=4$, $g=1.2$ in (<ref>), and the Hermitian operator $A=Z\otimes I_{2^{L-1}}$. In our numerical experiment, we set $\tau=1$ and $r=1$ for the discrete-time Lindblad dynamics simulation, and $\tau=0.1$ for continuous-time Lindblad dynamics simulation. The stopping times are set to $T=80$ for TFIM-4. Each Lindblad dynamics simulation starts from the initial state with zero overlaps between the ground state and is repeated 100 times, and we calculate the average energy and overlap with the ground state. The findings are depicted in Figure <ref>. Our observations indicate that both Lindblad dynamics exhibit efficient convergence to the ground state starting from the initial state with zero overlaps. Furthermore, the discrete-time Lindblad dynamics ($\tau=1$) can achieve a one-order reduction in total Hamiltonian simulation time compared to the approximate continuous-time Lindblad dynamics ($\tau=0.1$).

[Lindblad simulation time vs energy]
[Lindblad simulation time vs overlap]

[Hamiltonian simulation time vs energy]
[Hamiltonian simulation time vs overlap]
Continuous vs discrete-time Lindblad dynamics for TFIM-4.

§.§ Hubbard model

Consider the one-dimensional Hubbard model defined on $L$ spinful sites with open boundary conditions

\[
\]

Here $c_{j,\sigma}(c^\dagger_{j,\sigma})$ denotes the fermionic annihilation (creation) operator on the site $j$ with spin $\sigma$. $\left\langle\cdot,\cdot\right\rangle$ denotes sites that are adjacent to each other. $n_{j,\sigma}=c^\dagger_{j,\sigma}c_{j,\sigma}$ is the number operator.

We choose $L=4$, $t=1$, $U=4$, and $A=X\otimes I_{2^{L-1}}$. We set $\tau=0.5$ and $r=2$ for discrete-time Lindblad dynamics simulation and $\tau=0.025$ for continuous-time Lindblad dynamics simulation. The stopping times are set to $T=100$, and each Lindblad dynamics simulation is repeated 100 times starting from an initial state with zero overlap between the ground state. The results are presented in Figures <ref>. In our observations, we find that in both dynamics, the energy decreases to $\lambda_0$ and the overlap with the ground state increases to $1$. We observe that both the continuous-time and discrete-time Lindblad dynamics exhibit a fast convergence rate over time and achieve an accurate ground state construction.

[Lindblad simulation time vs energy]
[Lindblad simulation time vs overlap]

[Hamiltonian simulation time vs energy]
[Hamiltonian simulation time vs overlap]
Continuous vs discrete-time Lindblad dynamics for Hubbard-4.