Large time-step discretisation of adiabatic quantum dynamics

TL;DR

This paper demonstrates that large, uniform time steps can be used in discretizing adiabatic quantum dynamics, significantly reducing complexity while maintaining accuracy, especially with boundary cancellation conditions.

Contribution

It introduces a novel analysis showing large, uniform time steps are feasible in adiabatic quantum simulation, improving efficiency without sacrificing accuracy.

Findings

Large time steps reduce simulation complexity.

Uniform time steps are effective regardless of convergence order.

Boundary cancellation exponentially suppresses diabatic errors.

Abstract

Adiabatic quantum computing is a general framework for preparing eigenstates of Hamiltonians on quantum devices. However, its digital implementation requires an efficient Hamiltonian simulation subroutine, which may introduce extra computational overhead or complicated quantum control logic. In this work, we show that the time step sizes in time discretization can be much larger than expected, and the overall complexity is greatly reduced. Remarkably, regardless of the general convergence order of the numerical method, we can choose a uniform time step size independent of tolerated error and evolution time for sufficiently accurate simulation. Furthermore, with the boundary cancellation condition where the continuous diabatic errors are exponentially suppressed, we provide strong evidence on an exponential convergence of even first-order Trotter with uniform time step size. We apply our analysis to the example of adiabatic unstructured search and show several preferable features of the Trotterized adiabatic approach: it can match the Grover lower bound, it does not require a priori knowledge on the number of marked states, and its performance can be asymptotically comparable with that of the quantum approximate optimization algorithm.

Figures

Tables

Table 1. Table 1: Summary of the complexity estimates on discretizing  Eq.   1 . Query complexity is measured by the number of time propagation steps required in the evolution. We assume the gap condition only for the time-dependent Hamiltonian H ​ ( t ) H(t) , and that both | f ′ ​ ( s ) | |f^{\prime}(s)| and | f ′′ ​ ( s ) | |f^{\prime\prime}(s)| are uniformly bounded. Here α = ‖ H 0 ‖ + ‖ H 1 ‖ \alpha=\|H_{0}\|+\|H_{1}\| , α ~ p = ∑ γ 0 , ⋯ , γ p ∈ { 0 , 1 } ‖ [ H γ p , ⋯ , [ H γ 1 , H γ 0 ] ] ‖ \widetilde{\alpha}_{p}=\sum_{\gamma_{0},\cdots,\gamma_{p}\in\left\{0,1\right\}}\|[H_{\gamma_{p}},\cdots,[H_{\gamma_{1}},H_{\gamma_{0}}]]\| , T T is the evolution time, ϵ \epsilon is the tolerated error, and Δ ∗ \Delta_{*} represents the minimum gap of H ​ ( s ) H(s) .

Method	Error bound	Time step size	Scaling
First-order exponential integrator	Standard	$𝒪\left(ϵ{\alpha }^{-1}\right)$	$𝒪\left({\alpha }^3{\Delta }_{\ast }^{-3}{ϵ}^{-2}\right)$
First-order product formula	Standard	$𝒪\left(ϵ{\alpha }^{-2}{T}^{-1}\right)$	$𝒪\left({\alpha }^6{\Delta }_{\ast }^{-6}{ϵ}^{-3}\right)$
$p$-th order product formula	Standard	$𝒪\left({ϵ}^{1/p}{\alpha }^{-1-1/p}{T}^{-1/p}\right)$	$𝒪\left({\alpha }^{3+3/p}{\Delta }_{\ast }^{-3-3/p}{ϵ}^{-1-2/p}\right)$
First-order exponential integrator	Corollary 4	${\alpha }^{-1}$	$𝒪\left({\alpha }^3{\Delta }_{\ast }^{-3}{ϵ}^{-1}\right)$
First-order product formula	Corollary 8	$𝒪\left(\min \{{\alpha }^{-1},{\Delta }_{\ast }^{1/2}{\tilde{\alpha }}_2^{-1/2}\}\right)$	$𝒪\left({\alpha }^3{\Delta }_{\ast }^{-3}{ϵ}^{-1}\max \{1,{\tilde{\alpha }}_2^{1/2}{\Delta }_{\ast }^{-1/2}{\alpha }^{-1}\}\right)$
$p$-th order simplified product formula	Corollary 9	$𝒪\left(\min \{{\alpha }^{-1},{\Delta }_{\ast }^{1/p}{\tilde{\alpha }}_p^{-1/p}\}\right)$	$𝒪\left({\alpha }^3{\Delta }_{\ast }^{-3}{ϵ}^{-1}\max \{1,{\tilde{\alpha }}_p^{1/p}{\Delta }_{\ast }^{-1/p}{\alpha }^{-1}\}\right)$

Table 2. Table 2: Spectral gaps of H ​ ( s ) H(s) and W ​ ( s ) W(s) in the first example. Here the gap reported is the minimal gap over equi-distant points with step size 10 − 4 10^{-4} .

$ϵ$	${10}^{-1}$	$5\times {10}^{-2}$	$2\times {10}^{-2}$	${10}^{-2}$	$5\times {10}^{-3}$	$2\times {10}^{-3}$	${10}^{-3}$	$5\times {10}^{-4}$	$2\times {10}^{-4}$	${10}^{-4}$	0
Gap of $H\left(s\right)$	$5.1\times {10}^{-2}$	$2.3\times {10}^{-2}$	$7.9\times {10}^{-3}$	$3.0\times {10}^{-2}$	$5.6\times {10}^{-4}$	$8.9\times {10}^{-4}$	$1.4\times {10}^{-3}$	$1.6\times {10}^{-3}$	$1.8\times {10}^{-3}$	$1.8\times {10}^{-3}$	$1.9\times {10}^{-3}$
Gap of $W\left(s\right)$	$5.2\times {10}^{-2}$	$2.5\times {10}^{-2}$	$9.7\times {10}^{-3}$	$4.8\times {10}^{-2}$	$2.6\times {10}^{-3}$	$9.5\times {10}^{-4}$	$4.8\times {10}^{-4}$	$2.4\times {10}^{-4}$	$1.0\times {10}^{-4}$	$5.2\times {10}^{-5}$	$1.1\times {10}^{-16}$

Table 3. Table 3: Spectral gaps of H ​ ( s ) H(s) and W ​ ( s ) W(s) in the second example. Here the gap reported is the minimal gap over equi-distant points with step size 10 − 4 10^{-4} .

$ϵ$	${10}^{-1}$	$5\times {10}^{-2}$	$2\times {10}^{-2}$	${10}^{-2}$	$5\times {10}^{-3}$	$2\times {10}^{-3}$	${10}^{-3}$	$5\times {10}^{-4}$	$2\times {10}^{-4}$	${10}^{-4}$	0
Gap of $H\left(s\right)$	$5.1\times {10}^{-2}$	$2.5\times {10}^{-2}$	$9.6\times {10}^{-3}$	$4.7\times {10}^{-3}$	$2.4\times {10}^{-3}$	$9.4\times {10}^{-4}$	$4.7\times {10}^{-4}$	$2.3\times {10}^{-4}$	$1.0\times {10}^{-4}$	$5.1\times {10}^{-5}$	$3.3\times {10}^{-16}$
Gap of $W\left(s\right)$	$5.3\times {10}^{-2}$	$2.6\times {10}^{-2}$	$1.1\times {10}^{-2}$	$6.6\times {10}^{-3}$	$4.2\times {10}^{-3}$	$2.8\times {10}^{-3}$	$2.4\times {10}^{-3}$	$2.1\times {10}^{-3}$	$2.0\times {10}^{-3}$	$1.9\times {10}^{-3}$	$1.9\times {10}^{-3}$

Equations

$$i{\partial }_t|\psi (t)⟩=H(t/T)|\psi (t)⟩,0\le t\le T,$$

$$H(s)=(1-f(s))H_0+f(s)H_1,$$

$$\mathcal{T}{e}^{-i{\int }_{jh}^{(j+1)h}H(t/T)dt}=U_{\text{num}}((j+1)h,jh)+\mathcal{O}(h^{p+1}),$$

$$\mathcal{T}{e}^{-i{\int }_0^{T}H(t/T)dt}\approx \prod _{j=0}^{T_d-1}{U}_{\text{num}}((j+1)h,jh)$$

$$U_{\exp }(t+h,t)=\exp \left(-ihH(t/T)\right).$$

$$U_{\text{pf},1}(t+h,t)=\exp \left(-ihf(t/T)H_1\right)\exp \left(-ih(1-f(t/T))H_0\right),$$

$$U_{\text{pf},2}(t+h,t)=\exp \left(-ih/2(1-f((t+h/2)/T))H_0\right)\exp \left(-ihf((t+h/2)/T)H_1\right)\exp \left(-ih/2(1-f((t+h/2))/T)H_0\right).$$

$$U_{\text{pf},p}(t+h,t)=\prod _{k=0}^{K_p}\exp \left(-i{\beta }_{p,k}hf((t+{\delta }_{p,k}h)/T)H_1\right)\exp \left(-i{\alpha }_{p,k}h(1-f((t+{\gamma }_{p,k}h)/T))H_0\right)$$

$$U_{\text{spf},p}(t+h,t)=\prod _{k=0}^{K_p}\exp \left(-i{\beta }_{p,k}hf(t/T)H_1\right)\exp \left(-i{\alpha }_{p,k}h(1-f(t/T))H_0\right).$$

$$U(n/T_d)|\psi ⟩=\prod _{j=0}^{n-1}W(j/T_d)|\psi ⟩,1\le n\le T_d.$$

$$\Vert D^{(k)}W(s)\Vert \le \frac{c_k(s)}{T_d^k}$$

$${\sigma }_P^{(k)}(s)\supseteq \bigcup _{l=0}^k{\sigma }_P(s+l/T_d),{\sigma }_Q^{(k)}(s)\supseteq \bigcup _{l=0}^k{\sigma }_Q(s+l/T_d).$$

$$U_A(n/T_d)|\psi ⟩=\prod _{j=0}^{n-1}{W}_A(j/T_d)|\psi ⟩,1\le n\le T_d.$$

$$U_A(n/T_d)P(0)=P(n/T_d)U_A(n/T_d)$$

$$\Vert U(n/T_d)-U_A(n/T_d)\Vert \le \frac{C}{T_d}\left(\frac{{\hat{c}}_1(0)}{{\check{\Delta }}_2(0{)}^2}+\frac{{\hat{c}}_1(n/T_d)}{{\check{\Delta }}_2(n/T_d{)}^2}+\sum _{j=0}^{n-1}\frac{{\hat{c}}_1(j/T_d{)}^2}{T_d{\check{\Delta }}_2(j/T_d{)}^3}+\sum _{j=0}^{n-1}\frac{{\hat{c}}_2(j/T_d)}{T_d{\check{\Delta }}_2(j/T_d{)}^2}\right).$$

$$\mathcal{O}\left(\frac{h}{T}\left(\frac{h(\Vert H_0\Vert +\Vert H_1\Vert )}{{\tilde{\Delta }}_{\ast }^2}+\frac{h^2(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{{\tilde{\Delta }}_{\ast }^2}+\frac{h^2(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{{\tilde{\Delta }}_{\ast }^3}\right)\right).$$

$$W(j/T_d)=e^{-ihH(j/T_d)}=e^{-ihH(jh/T)}=U_{\exp }((j+1)h,jh).$$

$$c_1(s)=\mathcal{O}(h(\Vert H_0\Vert +\Vert H_1\Vert )),c_2(s)=\mathcal{O}(h(\Vert H_0\Vert +\Vert H_1\Vert )+h^2(\Vert H_0\Vert +\Vert H_1\Vert {)}^2).$$

$$T=\mathcal{O}\left(\frac{(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{{\Delta }_{\ast }^3ϵ}\right),$$

$$h=\frac{1}{\Vert H_0\Vert +\Vert H_1\Vert },$$

$$T_d=\mathcal{O}\left(\frac{(\Vert H_0\Vert +\Vert H_1\Vert {)}^3}{{\Delta }_{\ast }^3ϵ}\right).$$

$$T=C\frac{(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{{\Delta }_{\ast }^3ϵ}$$

$$T_d=T/h=C\frac{(\Vert H_0\Vert +\Vert H_1\Vert {)}^3}{{\Delta }_{\ast }^3ϵ}\ge C\frac{(\Vert H_0\Vert +\Vert H_1\Vert {)}^3}{{\Delta }_{\ast }^3}.$$

$$T_d\ge \frac{C}{4}\frac{\Vert H_0\Vert +\Vert H_1\Vert }{{\Delta }_{\ast }}=\frac{C}{4h{\Delta }_{\ast }}.$$

$$\mathcal{O}\left(\frac{h}{T}\left(\frac{h(\Vert H_0\Vert +\Vert H_1\Vert )}{h^2{\Delta }_{\ast }^2}+\frac{h^2(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{h^2{\Delta }_{\ast }^2}+\frac{h^2(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{h^3{\Delta }_{\ast }^3}\right)\right)$$

$$\le \mathcal{O}\left(\frac{1}{T}\left(\frac{\Vert H_0\Vert +\Vert H_1\Vert }{{\Delta }_{\ast }^2}+\frac{h(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{{\Delta }_{\ast }^2}+\frac{(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{{\Delta }_{\ast }^3}\right)\right)$$

$$\le \mathcal{O}\left(\frac{(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{T{\Delta }_{\ast }^3}\right).$$

$$T=\mathcal{O}\left(\frac{(\Vert H_0\Vert +\Vert H_1\Vert {)}^2}{{\Delta }_{\ast }^3ϵ}\right).$$

$$h=\frac{1}{\Vert H_0\Vert +\Vert H_1\Vert }$$

$$T_d\sim \frac{(\Vert H_0\Vert +\Vert H_1\Vert {)}^3}{{\Delta }_{\ast }^3ϵ}=\frac{1}{ϵ}{\left(\frac{(\Vert H_0\Vert +\Vert H_1\Vert )}{{\Delta }_{\ast }}\right)}^2\frac{\Vert H_0\Vert +\Vert H_1\Vert }{{\Delta }_{\ast }},$$

$${}_a⟨0|V_j|0{⟩}_a=\frac{1}{{\eta }_j}{H}_j.$$

Full text

1. I Introduction

2. II Summary of the results

3. II.1 Setup and numerical methods

4. II.2 Main results

5. II.2.1 Improved complexity

6. II.2.2 Improved complexity with boundary cancellation condition

7. II.2.3 Advantage of Trotterized evolution over continuous AQC

8. II.2.4 Discrete adiabatic Grover search and implication on QAOA

9. III Improved time step size and complexity estimate

10. III.1 Key idea

11. III.2 A review of the discrete adiabatic theorem

12. III.3 First-order exponential integrator

13. III.3.1 General case with gap assumption on the walk operator

14. III.3.2 General case without gap assumption on the walk operator

15. III.3.3 Implementation

16. III.4 Product formula

17. III.4.1 General case with gap assumption on the walk operator

18. III.4.2 General case without gap assumption on the walk operator

19. IV Improved time step size and complexity estimate with boundary cancellation

20. IV.1 A review of high-order discrete adiabatic theorem

21. IV.2 Our results

22. V Gap conditions of product formula and discrete adiabatic evolution for gapless Hamiltonian

23. VI Optimized gap dependence with application to Grover search and comparison with QAOA

24. VII Conclusions

25. A Continuous adiabatic theorems

26. B Multistep gap condition

27. C Bounds for finite differences

28. D Numerical methods and their standard complexity estimates

29. D.1 Standard time discretization error bounds

30. D.2 Standard complexity estimates

31. E Complexity estimate of the high-order simplified product formulae without gap condition on the walk operator

32. F High-order discrete adiabatic theorem

33. F.1 A slight variant

34. F.2 Missing steps in the proof of the high-order discrete adiabatic theorem and a numerical validation

35. G Spectral gap of the walk operator in adiabatic Grover search

36. G.1 Direct computations

37. G.2 Proof of Lemma 14

38. H Proof of Theorem 15

39. H.1 $p=1$

40. H.2 $1<p<2$

41. I Proof of Theorem 16

I Introduction

Preparing eigenstates of Hamiltonians is a fundamental task in quantum physics, quantum chemistry, and quantum information. Classical algorithms have a cost that scales polynomially with the dimension of the Hilbert space, which motivates the exploration of quantum algorithms for potential speedup. An approach to solving the eigenstate preparation problem is adiabatic quantum computing (AQC) Albash and Lidar (2018). In the framework of AQC, we consider the dynamics of a time-dependent Hamiltonian given by

[TABLE]

where $H(s)=(1-f(s))H_{0}+f(s)H_{1}$ is an interpolating Hamiltonian that smoothly evolves from an initial Hamiltonian $H_{0}$ to a target Hamiltonian $H_{1}$, and $f(s)$ is called the scheduling function such that $f(0)=0$ and $f(1)=1$. Typically, $H_{0}$ is designed to be simple with easily preparable eigenstates. Then, under the gap condition that the desired eigenpath of $H(s)$ is separated from the rest of the spectrum, the dynamics of Eq. 1 will approximately drive the system from the eigenstate of $H_{0}$ to that of $H_{1}$, as $T\rightarrow\infty$.

While AQC is a promising technique for solving general eigenstate preparation problems, its implementation on digital quantum devices requires an efficient subroutine for simulating Eq. 1, which may introduce non-optimal asymptotic scalings. Commonly used methods include the first-order exponential integrator Hochbruck and Ostermann (2010); Schiffer et al. (2022) and the first-order product formula Farhi et al. (2000); van Dam et al. (2001); Childs et al. (2021); An et al. (2021). Both methods, despite their simplicity, have only first-order convergence, and thus make the number of discrete evolution steps super-linear in the evolution time $T$ (for first-order methods, it is specifically quadratic in $T$). This can be expensive since $T$ is usually large to minimize diabatic excitation in AQC dynamics. Alternatively, high-order time-dependent Hamiltonian simulation algorithms, such as the high-order Trotter-Suzuki formula Suzuki (1993); Childs et al. (2021) and the truncated Dyson series method Kieferová et al. (2019); Low and Wiebe (2018), can reduce the number of evolution steps to almost linear in $T$. However, these methods may still introduce extra computational complication. The number of exponentially many local operators in the product formulae may grow rapidly as the convergence order increases, and the truncated Dyson series method requires complicated quantum control logic Su et al. (2021); Rajput et al. (2022).

In this work, we develop general methodology of analyzing numerical discretization methods for digitally simulating near adiabatic dynamics and obtain greatly improved error bounds. According to our analysis, the time step sizes in discretizing near-adiabatic dynamics can be much larger than expected, and the overall complexity is greatly reduced. Specifically, we consider the first-order exponential integrator and the product formula for time discretization. We show that, as long as the time step size is smaller than a uniform threshold independent of tolerated error $\epsilon$ and evolution time $T$, time discretization errors will be non-dominant and the digital simulation can achieve sufficient accuracy. Furthermore, under an additional assumption called the boundary cancellation condition, we show that in practice even a first-order product formula with uniform time step size can achieve exponential convergence, which achieves an exponential speedup in precision over previous results.

Our analysis is based on the discrete adiabatic theorem, which are discrete analog of the continuous adiabatic theorems and guarantee that sequentially applying a set of slowly varying gapped unitary walk operators can approximately drive the system from the eigenstate of the initial walk operator to the connected eigenstate of the final walk operator. In our analysis, we treat the numerical integrators as the set of discrete walk operators and show that the discretized dynamics indeed follows an adiabatic path with the same initial and final eigenstates as the continuous one but different in the middle. This not only allows us to obtain aforementioned improved error bounds, but also shows that the discretized evolution can potentially solve eigenstate preparation problems which are beyond the reach of the original continuous formulation. For example, if the Hamiltonian $H(s)$ becomes gapless, then the continuous AQC based on Eq. 1 cannot follow the desired eigenpath. However, we find that, even for a gapless Hamiltonian, it is possible that the first-order Trotterized operator with large time step size satisfies the gap condition, so the corresponding discretized evolution may still succeed in preparing the desired eigenstates.

As an application, we use our analysis to design a discrete adiabatic algorithm for the unstructured search problem with multiple marked items. Our algorithm has several preferable features: it can match the Grover lower bound, it does not require a priori knowledge on the number of marked states, and it can achieve a logarithmic scaling in the inverse precision. Additionally, we contribute to the active discussions on the relation between AQC and quantum approximate optimization algorithm (QAOA) Farhi et al. (2014), by showing that our discrete-adiabatic-based algorithm gives a set of QAOA angles which can achieve the optimal query complexity $\mathcal{O}(\sqrt{N})$ for finding the single marked state among $N$ states.

There have been a few recent efforts on improving the analysis of discretizing AQC. For instance, Yi (2021) demonstrates the robustness of digital adiabatic simulation using first-order Trotter formula with large time step sizes by applying continuous adiabatic theorems to an effective Hamiltonian. Additionally, a more recent work Kovalsky et al. (2023) proves a similar result by relating the first-order Trotter error estimate to its second-order counterpart. Our results improve over prior work in the following ways:

Our analysis is applicable to the first-order exponential integrator and product formulae of any order. As a comparison, previous works Yi (2021); Kovalsky et al. (2023) only focus on the first-order Trotter formula, and their analysis heavily relies on its specific structure. 2. 2.

Under the boundary cancellation condition, our analysis can give an exponentially better complexity. Specifically, our analysis explains why first-order methods with large time step sizes are very likely to have an exponential time convergence in this case, while previous works Yi (2021); Kovalsky et al. (2023) still only show a linear convergence as it is unclear how their approaches can be naturally combined with the additional boundary cancellation condition. 3. 3.

Our analysis uncovers the intrinsic property of time discretization of the AQC as following another (discrete) adiabatic path. This is not only conceptually different from the standard numerical analysis which views time discretization only as introducing numerical errors and driving the dynamics away from the adiabatic path, but also the key for revealing new findings. For example, with the discrete adiabatic point of view, we are able to show that Trotterized adiabatic evolution may still succeed even for a gapless Hamiltonian, but all the previous works Yi (2021); Kovalsky et al. (2023) only apply to gapped Hamiltonians.

The rest of this paper is organized as follows. We first summarize our main results in Section II. In Section III, we show our improved complexity estimates based on the discrete adiabatic theorem in the general case, and Section IV is devoted to the case with the boundary cancellation condition. Section V discusses the possibility of using the first-order product formula with large step for gapless Hamiltonians. In Section VI, we present our results for robust adiabatic Grover search and compare it with QAOA. Section VII summarizes this work and discusses potential further work.

II Summary of the results

II.1 Setup and numerical methods

We first formally restate the state preparation problem.

Problem 1** (State preparation).**

Let $H_{1}$ be a Hamiltonian and $|\phi\rangle$ be a target eigenstate of $H_{1}$. For $\varepsilon>0$, the goal of the state preparation problem is to prepare a quantum state $|\widetilde{\phi}\rangle$ such that $\||\widetilde{\phi}\rangle\langle\widetilde{\phi}|-|\phi\rangle\langle\phi|\|\leq\varepsilon$.

In the framework of AQC, given another simple Hamiltonian $H_{0}$ with known eigenstate, we consider the dynamics of Eq. 1 with interpolating Hamiltonian

[TABLE]

and the scheduling function $f(s)$. Under the gap condition, starting from the eigenstate of $H_{0}$, the solution of Eq. 1 at the final time will be a good approximation of the corresponding eigenstate of $H_{1}$.

The AQC dynamics of Eq. 1 needs further time discretization to be implemented on a digital quantum device. The most commonly used strategy for numerical propagation is to divide the time interval into $T_{d}$ equidistant sub-intervals and approximate the exact local evolution operator by a local unitary numerical propagator. Specifically, a $p$-th order local unitary evolution operator $U_{\text{num}}$ is constructed such that for any $0\leq j\leq T_{d}-1$, it approximates the time-ordering operator from $jh$ to $(j+1)h$ as follows

[TABLE]

and the global exact evolution operator can be approximated as

[TABLE]

where $T_{d}$ is the number of discrete time steps, and $h=T/T_{d}$ is the step size.

In practice, there are two common classes of numerical methods, namely the first-order exponential integrator and the product formulae. They differ in the level of time discretisation. If we only approximate the time-ordering operator by a time-independent Hamiltonian simulation problem, we obtain the first-order exponential integrator, which is given by

[TABLE]

In this method, each propagation step only involves one matrix exponential but the exponent is still an interpolation of the Hamiltonians $H_{0}$ and $H_{1}$. If we further split the evolution of $H(t/T)$ into exponentials of $H_{0}$ and $H_{1}$ separately, we obtain the product formulae. The first- and second-order product formulae are given as

[TABLE]

and

[TABLE]

Arbitrary $p$-th order formulae can be constructed via several approaches, such as Trotter-Suzuki recursion Suzuki (1993); Wiebe et al. (2010), Yoshida method Yoshida (1990) and their improved versions Morales et al. (2025); Blanes et al. (2024). The resulting $p$-th order formula can be written in the general form as

[TABLE]

for some real parameters $\alpha_{p,k},\beta_{p,k},\gamma_{p,k},\delta_{p,k}$. We will also consider a simplified higher-order product formula

[TABLE]

Equation (9) is obtained by replacing all the discrete time points Eq. 8 with the fixed initial time. If the function $f$ is a constant function (which is not the case for the adiabatic dynamics), then Eq. 9 represents a $p$-th order time-independent product formula.

II.2 Main results

II.2.1 Improved complexity

In this work, we propose a general framework for analyzing algorithms to simulate near adiabatic dynamics and obtain greatly improved error bounds in Section III. Our results and a comparison with previous results are summarized in Table 1.

We first state our results for the first-order exponential integrator. In Section III.3, we demonstrate that this first-order method can utilize a time step size as large as $\alpha^{-1}$ without introducing significant discretization errors, where $\alpha=\|H_{0}\|+\|H_{1}\|$. In comparison, the standard error bound typically requires a time step size of $\mathcal{O}(\epsilon\alpha^{-1})$, where $\epsilon$ is the tolerated level of errors (see Appendix D). Therefore, our analysis greatly increases the allowed time step size for sufficiently accurate simulation by fully excluding the error dependence in the time step size. Since the overall evolution time $T$ should be $\mathcal{O}(\alpha^{2}\Delta_{*}^{-3}\epsilon^{-1})$ to mitigate diabatic errors Jansen et al. (2007), where $\Delta_{*}$ is the minimum gap of $H(s)$, our result reduces the number of the total propagation steps from $\mathcal{O}(\alpha^{3}\Delta_{*}^{-3}\epsilon^{-2})$ to $\mathcal{O}(\alpha^{3}\Delta_{*}^{-3}\epsilon^{-1})$.

The main theoretical tool in our analysis is the discrete adiabatic theorem, which can be viewed as a discrete analog of the continuous adiabatic theorem. The theorem was first established in Dranov et al. (1998), and an improved version with explicit gap dependence was introduced in the recent work Costa et al. (2022). In the discrete adiabatic theorem, we are provided with a set of slowly varying unitary operators, along with a promise of the gap condition of these unitary operators. Sequentially applying these operators allows us to steer the system from the eigenstate of the initial unitary to the eigenstate of the final unitary, with a diabatic error bounded by $\mathcal{O}(1/T_{d})$, where $T_{d}$ is the number of the unitary operators. In our analysis, we treat the numerical integrators as a set of discrete walk operators. This allows us to directly apply the discrete adiabatic theorems to bound the error between the actual state and the exact eigenstates, making our approach tighter and more straightforward compared to standard analyses where errors are typically bounded by the sum of continuous diabatic errors and time discretization errors. The gap condition of the numerical integrator $e^{-ihH(t/T)}$ can be guaranteed from the gap condition of $H(s)$, as long as the step size $h\leq\alpha^{-1}$. The changes between adjacent numerical integrators are determined by the scheduling function $f(t/T)$, which is already sufficiently slow with large evolution time $T$ and suppresses the discrete adiabatic errors. This is why further reduction of the time step size is unnecessary.

Our analysis based on the discrete adiabatic theorem also applies to the high-order product formula. However, unlike the first-order exponential integrator case, now the gap condition of the Trotter numerical integrators can be very different from the gap condition of the Hamiltonian $H(s)$. For example, even for the first-order Trotter method $e^{-ihf(s)H_{1}}e^{-ih(1-f(s))H_{0}}$, it is not clear how its gap is related to the gap condition of the Hamiltonian $H(s)$ when the step size $h$ is large. It becomes even more complicated when we consider higher-order Trotter in Eq. 8, because the exponentials are the time evolution of Hamiltonian $H(s)$ evaluated at different times. Therefore, in Section III.4, we instead consider the simplified higher-order product formula in Eq. 9 with all the Hamiltonians evaluated at the same time. Our strategy is to reduce the step size $h$ such that the numerical operator $U_{\text{spf},p}$ is close to the walk operator $e^{-ihH(t/T)}$ by the distance at most $c\Delta_{*}$ for a small constant $c$, then the gap of $H(s)$ would imply the gap of $U_{\text{spf},p}$. Such a step size can be determined using the time-independent Trotter error bound in Childs et al. (2021). An important observation is that the sole purpose of reducing the time step size is to ensure the gap condition, so the time step size depends only on the norms and gaps of the Hamiltonians, but is independent of the overall evolution time $T$ and the tolerated error $\epsilon$.

Combining this with the discrete adiabatic theorem, we show in Section III.4 that, in a $p$-th simplified product formula, a sufficiently accurate simulation may be obtained with time step size $h=\mathcal{O}(\min\{\alpha^{-1},\Delta_{*}^{1/p}\widetilde{\alpha}_{p}^{-1/p}\})$, where $\widetilde{\alpha}_{p}=\sum_{\gamma_{0},\cdots,\gamma_{p}\in\left\{0,1\right\}}\|[H_{\gamma_{p}},\cdots,[H_{\gamma_{1}},H_{\gamma_{0}}]]\|$ is the scaling of nested commutators of depth $p$. Again, our analysis shows that the time step size can be independent of $T$ and $\epsilon$, while the best previous analysis Wiebe et al. (2010) gives the estimate of $h=\mathcal{O}(\epsilon^{1/p}\alpha^{-1-1/p}T^{-1/p})$ which still depends on $\epsilon$ and $T$. As a result, our results still improve the standard complexity estimate by suggesting a larger time step size and smaller number of propagation steps for accurate simulation. Another interesting remark is that our analysis implies the advantage of higher-order methods from a completely different perspective – the advantage is because higher-order methods can better match the gap condition with larger time step sizes, instead of better $\epsilon$ and $T$ dependence as suggested by the standard numerical analysis.

II.2.2 Improved complexity with boundary cancellation condition

We have shown that the overall complexity of simulating near adiabatic dynamics using the first-order exponential integrator or product formula is $\mathcal{O}(1/\epsilon)$, omitting dependence on other parameters. This linear dependence is due to the scaling of the adiabatic error in general case, namely $\mathcal{O}(1/T)$ in the continuous adiabatic theorem and $\mathcal{O}(1/T_{d})$ in the discrete adiabatic theorem. In Section IV, we show how the overall complexity can be further improved by taking advantage of high-order adiabatic theorems, under an additional assumption which is known as the boundary cancellation condition.

The boundary cancellation condition assumes the scheduling function $f(s)$ has vanishing derivatives of any order at the boundary $s=0,1$. Then, in the continuous case, the diabatic error can be suppressed by the high-order adiabatic theorems to a super-polynomial scaling $\mathcal{O}(1/T^{k})$ for any positive integer $k$ Jansen et al. (2007), or even an exponential convergence $\mathcal{O}(e^{-cT})$ Nenciu (1993); Ge et al. (2016). In the discrete case, the work Dranov et al. (1998) also claims a high-order discrete adiabatic theorem, which has the error scaling $\mathcal{O}(1/T^{k})$ for any positive integer $k$. Again, in the time discretization of Eq. 1, we view the numerical integrators as a set of discrete walk operators. Then the high-order discrete adiabatic theorem directly guarantees that, if the Hamiltonian has constant spectral norm and gap, then it suffices to use an $\mathcal{O}(1)$ time step size and the overall number of steps is $\mathcal{O}(1/\epsilon^{o(1)})$. Such a high-order convergence holds for any numerical methods, even including first-order ones like the first-order exponential integrator and the first-order product formula.

However, we remark that the improved scaling $\mathcal{O}(1/\epsilon^{o(1)})$ we find here is not mathematically rigorously proved yet. This is because our analysis is based on the high-order discrete adiabatic theorem claimed in Dranov et al. (1998), but we find a missing step in its proof and this missing step seemingly cannot be fixed in a simple way (see Section F.2 for details). Nevertheless, in Section F.2 we show numerical evidence that supports the correctness of the high-order discrete adiabatic theorem, so our improved complexity would hold in practice.

II.2.3 Advantage of Trotterized evolution over continuous AQC

So far, we have given improved complexity estimates of discretized AQC for gapped Hamiltonian $H(s)$. Notice that in order to apply the discrete adiabatic theorem, what we need is the gap condition of the numerical integrators rather than the Hamiltonian $H(s)$. In the first-order product formula with large time step size, these two types of gap conditions can be very different. Our strategy so far is to moderately reduce the time step size $h$ such that the Trotter operator $e^{-ihf(s)H_{1}}e^{-ih(1-f(s))H_{0}}$ is sufficiently close to the walk operator $e^{-ihH(s)}$, whose spectral gap is closely related to that of the Hamiltonian $H(s)$.

It seems that the difference of the gap conditions between the interpolating Hamiltonian $H(s)$ and the unitary numerical operators only introduces extra computational overheads, as reduced time step size is required to match the gap conditions. However, in Section V, we discuss an interesting phenomenon that such difference of the gap conditions can actually be advantageous when the interpolating Hamiltonian has an exponentially small spectral gap or becomes gapless. In such cases, Trotterized evolution with large time step size may significantly outperform the continuous AQC and even becomes the only way for adiabatic state preparation. More specifically, we consider the case where the gap of $H(s)$ is exponentially small or just becomes [math]. Then, continuous AQC requires an exponentially large evolution time or may even fail Avron and Elgart (1999). However, numerical discretization with a large time step size, such as the first-order Trotter method with time step size $1$ (that is, $e^{-if(s)H_{1}}e^{-i(1-f(s))H_{0}}$), may still have a large gap. Thus, discrete adiabatic theorems suggest that the large-step-size evolution may be able to find the desired eigenstates with much lower complexity. Notice that in this case, reducing the time step size further could make the situation worse, as the discrete evolution would be forced to approximate the continuous AQC that does not work. Therefore, for Hamiltonians with exponentially small gaps, their large-step-size discretization may surprisingly yield an exponential speedup compared to continuous AQC. We provide numerical evidence for this phenomenon using a toy model in Section V.

II.2.4 Discrete adiabatic Grover search and implication on QAOA

All of our aforementioned results have super-linear spectral gap dependence, as they hold generically and in the worst-case scenario. However, in continuous AQC, a linear gap dependence can be achieved by carefully designing an interpolation function that takes into account simultaneous gap information Roland and Cerf (2002); Jansen et al. (2007); Albash and Lidar (2018); An and Lin (2022). We observe that the linear gap dependence can be preserved with the same interpolating function when discretizing the dynamics using large time step size.

To illustrate this, in Section VI, we consider the example of adiabatic Grover search Roland and Cerf (2002), where we aim to find a uniform superposition of $M$ target states in an $N$-dimensional Hilbert space. We consider two possible interpolation functions: one generalized from Roland and Cerf (2002); Dalzell et al. (2017) that changes speed proportional to the gap size, and the second, a smooth connection of two identical glue functions, used in An and Lin (2022) to construct adiabatic quantum linear system solvers. The latter simultaneously slows down when gap closes and satisfies the boundary cancellation condition. We show that the AQC dynamics with either interpolation function can be discretized using a time step size $1$ without incurring extra discretization errors, and the overall number of steps using the second interpolation function can be $\widetilde{\mathcal{O}}(\sqrt{N/M})$ in dimension, almost recovering the quadratic Grover speedup, and $\mathcal{O}(1/\epsilon^{o(1)})$ in precision assuming the correctness of high-order discrete adiabatic theorems. Our second scheduling function provides a super-polynomial improvement in precision compared to previous work on adiabatic quantum search Dalzell et al. (2017). A remarkable feature of our algorithm is that it is robust against the estimations of $M$ in two aspects. First, constructing the interpolation functions does not require an estimate of $M$ (and for the second interpolation function, not even information about $N$). Additionally, increasing the number of propagation steps improves accuracy, similar to fixed-point iteration Yoder et al. (2014); Dalzell et al. (2017), instead of “overcooking” the problem as in the original Grover algorithm Grover (1996).

Another powerful approach for eigenstate preparation is the Quantum Approximate Optimization Algorithm (QAOA) Farhi et al. (2014). In QAOA, an ansatz is considered that alternatively evolves an input state with the target Hamiltonian $H_{1}$ and a mixing Hamiltonian $H_{0}$. The goal is to optimize the evolution times, also called QAOA angles, such that the final output state approximates the desired eigenstate. QAOA is closely related to AQC and especially discretized AQC, since for gapped Hamiltonian the Trotterization of AQC directly gives a candidate of the QAOA angles, yet QAOA is believed to have the potential to outperform AQC by resulting in a shorter overall evolution time and having better gap dependence. Even for the simple example of the Grover search problem, it was previously believed that the performance and analysis of QAOA should involve a very different process from estimating the gap dependence in the AQC approach Jiang et al. (2017). However, our improved AQC results demonstrate the existence of a Trotterized AQC with time step size $1$, which can readily solve the Grover search problem with a cost of $\mathcal{O}(\sqrt{N})$ in the case of a single marked state. Therefore, the performance of discretized AQC can match that of QAOA for this task.

III Improved time step size and complexity estimate

III.1 Key idea

We first explain the key components in our improved analysis. The first observation is that, unlike the standard analysis which uses the triangle inequality and bounds both the continuous adiabatic error and the time discretization error, we directly view the numerical propagation operator $U_{\text{num}}((j+1)h,jh)$ as the discrete adiabatic walk operator. Then, as long as $U_{\text{num}}$ satisfies the assumptions in the discrete adiabatic theorems and the eigenstate of the last numerical propagation operator coincides with our target eigenstate, the overall approximation error can be directly bounded using the discrete adiabatic theorems. This is a more straightforward and tighter way to bound the overall error. The second observation is that numerical operators with relatively large time step size suffice to fulfill the assumption of slow changes in the discrete adiabatic theorems. This is because the Hamiltonian $H(t/T)$ varies on the rescaled time $s=t/T$, which already controls the change of the walk operators at a low level and thus smaller time step sizes are not necessary.

A further result from our analysis is a better time discretization error bound for the first-order exponential integrator and the product formula. We may use the triangle inequality in a different way to bound the numerical discretization error by the sum of the continuous and discrete adiabatic errors, which are both well bounded even with large time step sizes.

III.2 A review of the discrete adiabatic theorem

Discrete adiabatic evolution is a discrete analog of continuous adiabatic evolution. In continuous adiabatic evolution (see Appendix A for a detailed discussion), we are given a time-dependent Schrödinger equation. If the time-dependent Hamiltonian varies slowly and there is a gap in its spectrum, then propagating the Schrödinger equation approximately drives the state from an eigenstate of the initial Hamiltonian to an eigenstate of the final Hamiltonian. In the discrete setting, we are given a set of unitary walk operators and sequentially apply those on an input quantum state. If two adjacent walk operators are close and satisfy the gap condition, then such discrete evolution will approximately drive the quantum state from the eigenstate of the first walk operator to the eigenstate of the final walk operator.

The discrete adiabatic theorem bounds the distance between the output of the discrete evolution and the ideal eigenspace of the final walk operator. It was first proved in Dranov et al. (1998), and a recent work Costa et al. (2022) improved the linearly convergent discrete adiabatic theorem to an explicit form with clear gap dependence. Here, we briefly review the result in Costa et al. (2022).

Let $W(s)$ be a set of parameterized unitary walk operators for $s\in[0,1]$. For an input quantum state $|\psi\rangle$ and a positive integer $T_{d}$, consider the discrete evolution $U(n/T_{d})$ defined as

[TABLE]

Suppose that the discrete evolution satisfies the following conditions.

For $k=1,2$, there exist real-valued functions $c_{k}(s)$ such that

[TABLE]

for $0\leq s\leq 1-k/T_{d}$, where $D^{(k)}W(s)$ represents the $k$-th order finite difference of $W(s)$ with step size $1/T_{d}$. 2. 2.

The eigenvalues of $W(s)$ can be separated into two groups $\sigma_{P}(s)$ and $\sigma_{Q}(s)$ such that for any $k=0,1,2$ and $0\leq s\leq 1-k/T_{d}$, the angular distance between the arcs $\sigma_{P}^{(k)}(s)$ and $\sigma_{Q}^{(k)}(s)$ is bounded from below by $\Delta_{k}(s)>0$, where

[TABLE] 3. 3.

The input state $|\psi\rangle$ is within the eigenspace corresponding to $\sigma_{P}(0)$.

The first of these conditions means that the rate of variation of the walk operator is bounded, and is a similar condition to the bound on the first and second time derivatives in the continuous case. The second condition is ensuring that there is a gap. For generality it is allowed to be a gap between eigenspaces, rather than just a ground state and excited states. There is a further subtlety in the discrete case that the gap is considered between successive steps as well; that is, the eigenspace $\sigma_{P}$ at one step cannot overlap with $\sigma_{Q}$ at the next step. This is accounted for by considering the union of the eigenspaces at three successive steps as in Eq. (12). However, we remark that if $T_{d}$ is large enough, then this multistep gap condition can be guaranteed by a uniformly bounded spectral gap of each $W(s)$. Intuitively, this is because the eigenvalues of $W(s)$ do not change much from one step to another, and thus different arcs of the spectrum at different time steps will not overlap. We discuss this with more technical details in Appendix B. The third condition just specifies $\sigma_{P}$ as the eigenspace that the state should remain within during the discrete adiabatic evolution.

For reference we also need to define an ideal adiabatic evolution operator $U_{A}(n/T_{d})$, which drives the input state to the ideal final eigenspace without any leakage. The operator $U_{A}(n/T_{d})$ is defined through a set of ideal adiabatic walk operators $W_{A}(s)$ such that

[TABLE]

The ideal adiabatic walk operators $W_{A}(s)$ include a counter-diabatic correction of the original walk operators $W(s)$. For technical simplicity, we omit the explicit definition of $W_{A}(s)$ here and refer to Dranov et al. (1998); Costa et al. (2022). An important proposition involving $U_{A}(n/T_{d})$ is that for all $1\leq n\leq T_{d}$,

[TABLE]

where $P(s)$ denotes the spectral projection onto the eigenspace $\sigma_{P}(s)$. The left-hand side corresponds to projecting onto the desired initial eigenspace, then evolving for $n$ steps, which should give the eigenspace at step $n$ if the evolution is adiabatic. That is indeed what the right-hand side shows, because it has the evolution for $n$ steps followed by the projection onto the eigenspace, so acting with this operator on any state must yield a state in the correct eigenspace. Therefore $U_{A}(n/T_{d})$ exactly preserves the simultaneous eigenspace of $W(n/T_{d})$ corresponding to $\sigma_{P}(n/T_{d})$.

Now we are ready to state the first discrete adiabatic theorem proved in Costa et al. (2022).

Lemma 2**.**

Suppose that $T_{d}\geq\sup_{s\in[0,1]}(4\hat{c}_{1}(s)/\check{\Delta}_{2}(s))$, where the hat (check) notations of a function at $s$ represent the maximum (minimum) of its function values at neighboring points $\{s-1/T_{d},s,s+1/T_{d}\}$ whenever well-defined. Then for any $1\leq n\leq T_{d}$, there exists an absolute constant $C$ such that

[TABLE]

Because $U_{A}(n/T_{d})$ is exact adiabatic evolution, it would map an initial state in the initial eigenspace to one in the desired final eigenspace. The bound on the difference of operators therefore gives a bound on how much the state can leak outside the eigenspace under the true evolution $U(n/T_{d})$.

III.3 First-order exponential integrator

III.3.1 General case with gap assumption on the walk operator

Now we consider walk operators that are given by Hamiltonian evolution for time $h$, $W(s)=e^{-ihH(s)}$, where $H(s)$ is varied in an adiabatic way as in Eq. (2). We can then use this form to determine $c_{k}(s)$ and provide a theorem for the error in terms of the properties of the Hamiltonians $H_{0},H_{1}$. For this theorem the gap is still given for the walk operator, though it will be implicitly dependent on $h$.

In general cases, we can obtain a linear convergence in the evolution time. We first prove the following result for arbitrary time step size $h$ under the multistep gap assumption for the walk operator.

Theorem 3**.**

Consider 1 using digital adiabatic simulation with first-order exponential integrator $W(s)=e^{-ihH(s)}$, where $H(s)$ is given as Eq. 2. Suppose that both $|f^{\prime}(s)|$ and $|f^{\prime\prime}(s)|$ are uniformly bounded over $[0,1]$. Let $\widetilde{\Delta}(s)$ denote the multistep spectral gap of $U_{\exp}(sT+h,sT)$, and $\widetilde{\Delta}_{*}=\min_{s\in[0,1-h/T]}\widetilde{\Delta}(s)$. Then, the overall error between the actual and the ideal evolution can be bounded by

[TABLE]

Proof.

Let $T_{d}=T/h$. For $s\in[0,1]$ and any $0\leq j\leq T_{d}-1$, we have

[TABLE]

We show in Lemma 22 that, for this walk operator $W(s)$, we can choose its corresponding $c_{1}(s)$ and $c_{2}(s)$ (defined in Eq. 11) as

[TABLE]

Then we can apply Lemma 2 to obtain the desired error bound. ∎

III.3.2 General case without gap assumption on the walk operator

Theorem 3 assumes the multistep gap condition for the first-order exponential integrator. While this may be established independently in some applications, a more common assumption is the (fixed-time) gap condition for the Hamiltonian $H(s)$. When the time step size is small enough, we may relate the multistep gap condition for $W(s)$ to the fixed-time gap condition for $H(s)$ and show the following result.

Corollary 4**.**

Consider 1 using digital adiabatic simulation with the first-order exponential integrator Eq. 5, $W(s)=e^{-ihH(s)}$. Suppose that both $|f^{\prime}(s)|$ and $|f^{\prime\prime}(s)|$ are uniformly bounded over $[0,1]$. Let $\Delta(s)$ denote the spectral gap of $H(s)$, and $\Delta_{*}=\min_{s\in[0,1]}\Delta(s)$. Then, in order to bound the overall error by $\epsilon\in(0,1)$, it suffices to choose

[TABLE]

and

[TABLE]

and the overall number of steps becomes

[TABLE]

Proof.

Following the notations in the proof of Theorem 3, we use $T_{d}=T/h$ to denote the number of steps and define the walk operator $W(s)=e^{-ihH(s)}$. Notice that when $h=1/(\|H_{0}\|+\|H_{1}\|)$, we always have $\|hH(s)\|\leq 1$, so the gap of $W(s)$ is exactly the same as that of $hH(s)$, which is bounded from below by $h\Delta_{*}$. Suppose that our choice of $T$ is

[TABLE]

for a constant $C>0$. Then we have

[TABLE]

Notice that $\Delta_{*}\leq 2\|H(s)\|\leq 2(\|H_{0}\|+\|H_{1}\|)$, so we have

[TABLE]

According to Lemma 22, we have $c_{1}(s)=\mathcal{O}(h(\|H_{0}\|+\|H_{1}\|))=\mathcal{O}(1)$, so we can choose a sufficiently large constant $C$ such that $T_{d}\geq\frac{2\pi}{h\Delta_{*}}\sup c_{1}(s)$. Then Lemma 20 ensures that the multistep gap of $W(s)$ is bounded from below by $h\Delta_{*}/2$.

Plugging this back into Theorem 3 yields the error bound

[TABLE]

In order to further bound the error by $\epsilon$, it suffices to choose

[TABLE]

Since the error bound is independent of the time step $h$, we directly choose the largest possible time step size

[TABLE]

which ensures the fixed-time gap condition of $W(s)$. Finally, the condition $T/h>2\widetilde{c}(\|H_{0}\|+\|H_{1}\|)/\Delta_{*}$ is naturally satisfied, since $\Delta_{*}\leq 2\|H(s)\|\leq 2(\|H_{0}\|+\|H_{1}\|)$ and

[TABLE]

from which we can also infer that $T_{d}\gtrsim\frac{\|H_{0}\|+\|H_{1}\|}{\Delta_{*}}$ and thus satisfies the condition for $T/h$. ∎

III.3.3 Implementation

So far we have shown an improved estimate on the number of steps in the first-order exponential integrator. However, it is still not straightforward to implement the first-order exponential integrator $e^{-ihH(s)}$ from the input models of $H_{0}$ and $H_{1}$. Here we briefly describe how to implement $e^{-ihH(s)}$. Suppose that we are given unitaries $V_{0}$ and $V_{1}$ such that, for $j=0,1$, $V_{j}$ is a $(\eta_{j},n_{a},0)$-block-encoding of $H_{j}$, i.e.,

[TABLE]

Here $\eta_{j}$ represents the block-encoding factor such that $\eta_{j}\geq\|H_{j}\|$. The number of the ancilla qubits for both block-encodings is the same as $n_{a}$ without loss of generality (if the number of the ancilla qubits are not the same, then we may define $n_{a}$ to be the larger one and supplement the other unitary with extra trivial ancilla qubits). We further assume that the block-encoding factors are asymptotically optimal in the sense that $\eta_{j}=\mathcal{O}(\|H_{j}\|)$ as well.

For every fixed $s\in[0,1]$, we first block encode $H(s)=(1-f(s))H_{0}+f(s)H_{1}$. This can be done using the linear combination of unitaries technique supplemented with an extra ancilla qubit. Let the rotation be

[TABLE]

and the select oracle be

[TABLE]

Then, according to (Gilyén et al., 2019a, Lemma 52), the operator $(\mathrm{H}\otimes I_{n_{a}}\otimes I_{n})\text{SEL}(R(s)\otimes I_{n_{a}}\otimes I_{n})$ is a $(1,n_{a}+1,0)$-block-encoding of

[TABLE]

To implement $e^{-ihH(s)}$, we may use the qubitization Low and Chuang (2019) or the quantum singular value transformation technique Gilyén et al. (2019a). In particular, write

[TABLE]

so we need to evolve the rescaled Hamiltonian $\frac{1}{\sqrt{2}\sqrt{(1-f(s))^{2}\eta_{0}^{2}+f(s)^{2}\eta_{1}^{2}}}H(s)$ to time $h\sqrt{2}\sqrt{(1-f(s))^{2}\eta_{0}^{2}+f(s)^{2}\eta_{1}^{2}}$, which is still $\mathcal{O}(1)$. This is because $h\sqrt{2}\sqrt{(1-f(s))^{2}\eta_{0}^{2}+f(s)^{2}\eta_{1}^{2}}\leq\sqrt{2}h(\eta_{0}+\eta_{1})$, which is $\mathcal{O}(1)$ due to the choice of $h$ in Corollary 4. Then, by (Gilyén et al., 2019a, Corollary 60), we may obtain a $(1,n_{a}+3,\epsilon^{\prime})$-block-encoding of $e^{-ihH(s)}$, using the unitaries $V_{0}$ and $V_{1}$ a total number of times $\mathcal{O}(\log(1/\epsilon^{\prime}))$. Finally, the overall numerical propagator can be obtained by multiplying all the block-encoded unitaries, which, by (Gilyén et al., 2019b, Corollary 55), yields a $(1,n_{a}+3,4(T/h)^{2}\epsilon^{\prime})$-block-encoding of it. To bound the overall error by $\epsilon$, it suffices to choose $\epsilon^{\prime}=\epsilon h^{2}/T^{2}$, and this contributes to an extra factor of $\log((\|H_{0}\|+\|H_{1}\|)/(\epsilon\Delta_{*}))$ in the complexity of the gate-based implementation.

III.4 Product formula

III.4.1 General case with gap assumption on the walk operator

We now consider the product formula

[TABLE]

and its simplified version

[TABLE]

We start with the simplest case where we assume the multistep gap condition for the product walk operator.

Theorem 5**.**

Consider 1 using digital adiabatic simulation with the product formula $U_{\text{pf},p}$ or the simplified product formula $U_{\text{spf},p}$. Suppose that both $|f^{\prime}(s)|$ and $|f^{\prime\prime}(s)|$ are uniformly bounded over $[0,1]$. Let $\widetilde{\Delta}(s)$ denote the multistep spectral gap of $U_{\text{pf},p}(sT+h,sT)$ or $U_{\text{spf},p}(sT+h,sT)$, and $\widetilde{\Delta}_{*}=\min_{s\in[0,1-h/T]}\widetilde{\Delta}(s)$. Then, the overall error between the actual and the ideal evolution can be bounded by

[TABLE]

Proof.

For $s\in[0,1-h/T]$ and any $0\leq j\leq T_{d}-1$, the walk operator $W(j/T_{d})$ is $U_{\text{pf},p}((j+1)h,jh)$ or $U_{\text{pfs},p}((j+1)h,jh)$. We show in Lemma 23 and Lemma 24 that, for the walk operator $W(s)$ in either case, we can choose $c_{1}(s)$ and $c_{2}(s)$ (defined in Eq. 11) as

[TABLE]

Then we can apply Lemma 2 to obtain the desired error bound. ∎

III.4.2 General case without gap assumption on the walk operator

Theorem 5 assumes the gap condition for the product walk operator, which is not natural when the time step size $h$ is large. When the gap condition is not assumed, we can indeed decrease the time step size to make the walk operator closer to the evolution $e^{-ihH(s)}$, of which the gap is closely related to that of $H(s)$. Here we only consider the simplified product formula Eq. 9.

We first prove the result for first- and second-order product formulae. An interesting phenomenon is that the first- and second-order formulae share the same eigenvalues as in the following lemma.

Lemma 6**.**

Consider the simplified first-order product formula $U_{\text{spf},1}(t+h,t)=e^{-ihf(t/T)H_{1}}e^{-ih(1-f(t/T))H_{0}}$ and the simplified second-order product formula $U_{\text{spf},2}(t+h,t)=e^{-ih(1-f(t/T))H_{0}/2}e^{-ihf(t/T)H_{1}}e^{-ih(1-f(t/T))H_{0}/2}$. The eigenvalues of $U_{\text{spf},1}(t+h,t)$ and $U_{\text{spf},2}(t+h,t)$ are exactly the same.

Proof.

Notice that $U_{\text{spf},2}$ is just a rotation of $U_{\text{spf},1}$ according to a unitary, i.e.,

[TABLE]

and such a rotation does not change eigenvalues. ∎

Now we can estimate the gap of the product walk operators by reducing the time step size and using perturbation theory.

Lemma 7**.**

For $s\in[0,1-h/T]$ and a scheduling function $f(s)$, let $H(s)=(1-f(s))H_{0}+f(s)H_{1}$ be the interpolating Hamiltonian, $U_{\text{spf},1}$ and $U_{\text{spf},2}$ be the simplified time-dependent product formulae as $U_{\text{spf},1}(sT+h,sT)=e^{-ihf(s)H_{1}}e^{-ih(1-f(s))H_{0}}$ and $U_{\text{spf},2}(sT+h,sT)=e^{-ih(1-f(s))H_{0}/2}e^{-ihf(s)H_{1}}e^{-ih(1-f(s))H_{0}/2}$. Let $\Delta_{H}(s)$ and $\Delta_{U,p}(s)$ denote the spectral gap of $H(s)$ and $U_{\text{spf},p}$, respectively. Then, for any $h\leq 1/(\|H_{0}\|+\|H_{1}\|)$ and $p=1,2$, we have

[TABLE]

Proof.

According to to Lemma 6, it suffices to consider only the second-order formula.

Notice that the spectral gaps of $e^{-ihH(s)}$ and $hH(s)$, though defined in slightly different ways, are of the same value for all $h\leq 1/(\|H_{0}\|+\|H_{1}\|)$. So it suffices to investigate the spectral gaps of $U_{\text{spf},2}(sT+h,sT)$ and $e^{-ihH(s)}$. For every fixed $s$, we may bound the difference between $U_{\text{spf},2}(sT+h,sT)$ and $e^{-ihH(s)}$ by using the error bound for the second-order time-independent Trotter method because the discrete time points are fixed in different exponentials. For example, by (Childs et al., 2021, Proposition 16),

[TABLE]

Equation (40) implies that the difference between $U_{\text{spf},2}(sT+h,sT)$ and $e^{-iH(s)}$ is small if we reduce the time step size $h$, and further implies that the eigenvalues of $U_{\text{spf},2}(sT+h,sT)$ and $e^{-iH(s)}$ are close according to Bhatia and Davis (1984); Elsner and He (1993). Specifically, according to the eigenvalue perturbation theorem of Refs. Bhatia and Davis (1984); Elsner and He (1993), there exists an ordering $\lambda_{j}(s)$ and $\widetilde{\lambda}_{j}(s)$ of the eigenvalues of $U_{\text{spf},p}(sT+h,sT)$ and $e^{-iH(s)}$, respectively, such that

[TABLE]

Now let $\mathbb{S}^{1}$ denote the unit circle, and $|z_{1}-z_{2}|_{\mathbb{S}^{1}}$ denote the angular distance between $z_{1},z_{2}\in\mathbb{S}^{1}$. Then we have

[TABLE]

Because the arcsine has a monotonically increasing first derivative, the inequality in Eq. (III.4.2) implies

[TABLE]

Let $\lambda_{0}(s)$ and $\lambda_{1}(s)$ denote the two eigenvalues that determine the gap. Then the gap $U_{\text{spf},2}(sT+h,sT)$ can be bounded by

[TABLE]

and

[TABLE]

where we use the fact again that the gaps of the Hamiltonian and the walk operator are the same. ∎

According to Lemma 7, we can reduce the time step size $h$ to make the product walk operator satisfy the gap condition as well. An important observation is that the time step size will not depend on the tolerated error $\epsilon$, because the only reason to reduce the time step size is for the gap condition, which is independent of the numerical error.

Now we are ready to prove the following complexity estimate for the product formula without assuming the gap condition of the product walk operator.

Corollary 8**.**

Consider 1 using digital adiabatic simulation with the simplified first- and second-order product formulae as in Eq. 9. Suppose that both $|f^{\prime}(s)|$ and $|f^{\prime\prime}(s)|$ are uniformly bounded over $[0,1]$. Let $\Delta(s)$ denote the spectral gap of $H(s)$, and $\Delta_{*}=\min_{s\in[0,1]}\Delta(s)$. Then, in order to bound the overall error between the actual and the ideal evolution by $\epsilon\in(0,1)$, it suffices to choose

[TABLE]

and

[TABLE]

and the overall number of steps becomes

[TABLE]

Proof.

Let $T_{d}=T/h$, and we assume $T_{d}$ is an integer. Define the walk operator

[TABLE]

We first bound the fixed-time gap of $W(s)$. According to the choice of $h$, we always have

[TABLE]

Using Lemma 7, the gap of $W(s)$ is bounded from below by the gap of $\frac{1}{2}h\Delta(s)$, which has a further lower bound $\frac{1}{2}h\Delta_{*}$.

For the multistep gap, we can use Lemma 20, and we need to verify that $T_{d}$ is sufficiently large. Specifically, suppose that our choice of $T$ is

[TABLE]

for a constant $C>0$. Then we have

[TABLE]

Notice that $\Delta_{*}\leq 2\|H(s)\|\leq 2(\|H_{0}\|+\|H_{1}\|)$, so we have

[TABLE]

According to Lemma 24, we have $c_{1}(s)=\mathcal{O}(h(\|H_{0}\|+\|H_{1}\|))=\mathcal{O}(1)$, so we can choose a sufficiently large constant $C$ such that $T_{d}\geq\frac{2\pi}{h\Delta_{*}/2}\sup c_{1}(s)$. Then Lemma 20 ensures that the multistep gap of $W(s)$ is bounded from below by $h\Delta_{*}/4$.

Therefore, the error bound in Theorem 5 becomes

[TABLE]

and the desired choices of $T$ and $h$ directly follow from this error bound and the conditions for the multistep gap. ∎

For high-order product formulae, we can estimate the complexity by the same strategy as the first- and second-order formulae and replacing the second nested commutators by higher-order nested commutators due to the high-order Trotter error bound in (Childs et al., 2021, Theorem 11). Now we state the result in the following corollary, and give its proof in Appendix E.

Corollary 9**.**

Consider 1 using digital adiabatic simulation with simplified $p$-th order product formulae as in Eq. 9. Suppose that both $|f^{\prime}(s)|$ and $|f^{\prime\prime}(s)|$ are uniformly bounded over $[0,1]$. Let $\Delta(s)$ denote the spectral gap of $H(s)$, and $\Delta_{*}=\min_{s\in[0,1]}\Delta(s)$. Then, for any $p\geq 1$,

the overall error between the actual and the ideal evolution can be bounded by

[TABLE] 2. 2.

in order to bound the overall error by $\epsilon$, it suffices to choose

[TABLE]

and

[TABLE]

and the overall number of steps becomes

[TABLE]

IV Improved time step size and complexity estimate with boundary cancellation

IV.1 A review of high-order discrete adiabatic theorem

The high-order discrete adiabatic theorem aims to establish a super-polynomial convergence in $1/T$ under further assumptions. In particular, we suppose the discrete evolution satisfies the following conditions.

For $k\geq 1$, there exist real-valued functions $c_{k}(s)$ such that $\|D^{(k)}W(s)\|\leq c_{k}(s)/T^{k}$ for $0\leq s\leq 1-k/T$, where $D^{(k)}W(s)$ represents the $k$-th order finite difference of $W(s)$ with step size $1/T$. 2. 2.

The eigenvalues of $W(s)$ can be separated into two groups $\sigma_{P}(s)$ and $\sigma_{Q}(s)$ such that each continuous eigenpath $\lambda(s)$ of $W(s)$ consistently belongs to one of the two groups and the distance between $\sigma_{P}(s)$ and $\sigma_{Q}(s)$ is bounded from below by $\Delta(s)>0$. 3. 3.

The input state $|\psi\rangle$ is within the eigenspace corresponding to $\sigma_{P}(0)$.

Under these conditions, Ref. Dranov et al. (1998) claims an upper bound of the discrete diabatic error, achieving $\mathcal{O}(1/T^{k})$ for any positive integer $k$. Unfortunately, after a careful examination, we find a missing step in the proof provided in Ref. Dranov et al. (1998), so the claimed high-order discrete adiabatic theorem is not rigorously supported. Nevertheless, we also perform a numerical test, verifying that the conclusion on the high-order convergence is numerically correct. Therefore, to be precise, here we present the result in Dranov et al. (1998) as a conjecture of which we have strong confidence on the correctness, and refer to Appendix F for more details on the missing step in the proof and our numerical validation.

Conjecture 10**.**

Suppose that $W(s)$ satisfies the boundary cancellation condition in the sense that $W^{(k)}(0)=W^{(k)}(1)=0$ for all $k\geq 1$. Then, for any integer $T$ and any positive integer $k$, we have

[TABLE]

where $U(s)$ is the discrete evolution defined in Eq. 10, $P(s)$ is the spectral projection of $W(s)$ onto $\sigma_{P}(s)$, and $C_{k}$ is a constant independent of $T$.

Similar to the continuous case, 10 bounds the leakage rate of the actual state out of the desired eigenspace. Compared with Lemma 2, 10 shows an error bound with super-polynomial convergence in $T$ under an additional assumption analogous to the continuous boundary cancellation condition. However, 10 does not give an explicit dependence on the spectral gap. It does not show how the pre-constant depends on the order $k$ either, which prevents us from obtaining an exponentially small error because there might be an exponential dependence on $k$ hidden in the constant $C_{k}$.

IV.2 Our results

We now focus on the scenario with the additional boundary cancellation condition and establish the super-polynomial convergence result. Due to the lack of a discrete adiabatic theorem with super-polynomial convergence and explicit gap dependence at the same time, it is hard to simultaneously establish both the gap dependence and the precision dependence. Notice that if we want to track only the gap dependence, we can assume the tolerated error to be at a constant level and directly use Corollary 4 or Corollary 9 in the general case, which gives an $\mathcal{O}(\Delta_{*}^{-3})$ dependence. Throughout this section, we assume that the spectral norm of the Hamiltonians $H_{0}$ and $H_{1}$ are bounded by $1$ and the minimum spectral gap to be bounded from below by a constant, and we focus on the dependence on $\epsilon$.

For the first-order exponential integrator and product formulae, we can prove the following two theorems, which show that for any method it suffices to choose an $\mathcal{O}(1)$ time step size and the number of steps can be $\mathcal{O}(1/\epsilon^{o(1)})$. Remarkably, such results also hold even for first-order numerical methods.

Theorem 11**.**

Consider 1 using exponential propagators Eq. 5. Suppose that 10 is true, both $\|H_{0}\|$ and $\|H_{1}\|$ are bounded by $1$, $H(s)$ satisfies the gap condition with constant-level minimum gap, and the smooth interpolation function has $f^{(k)}(0)=f^{(k)}(1)=0$ for all $k\geq 1$. Then, in order to prepare an $\epsilon$-approximation of the eigenstate of $H_{1}$, it suffices to choose $h=1$, and the evolution time and overall number of steps as

[TABLE]

for any $k\geq 1$.

Proof.

Let $W(s)=e^{-iH(s)}$ where $H(s)=(1-f(s))H_{0}+f(s)H_{1}$. We only need to verify that the walk operator $W(s)$ satisfies all the conditions of 10. Since the function $f(s)$ is smooth over $[0,1]$, the operator $W(s)$ is also smooth. The walk operator $W(s)$ also satisfies the gap condition since $H(s)$ is assumed to satisfy the gap condition and $\|H(s)\|\leq 1$. Finally, the boundary cancellation condition of $f(s)$ implies the boundary cancellation of $W(s)$. Therefore, our claim directly follows from 10. ∎

Theorem 12**.**

Consider 1 using simplified product formula as in Eq. 9. Suppose that 10 is true, both $\|H_{0}\|$ and $\|H_{1}\|$ are bounded by $1$, $H(s)$ satisfies the gap condition with constant-level minimum gap, and the smooth interpolation function $f(s)$ has $f^{(k)}(0)=f^{(k)}(1)=0$ for all $k\geq 1$. Then, in order to prepare an $\epsilon$-approximation of the eigenstate of $H_{1}$, it suffices to choose $h=\mathcal{O}(1)$, and the evolution time and overall number of steps as

[TABLE]

for any $k\geq 1$.

Proof.

Define the walk operator

[TABLE]

Since the function $f(s)$ is smooth over $[0,1]$, the operator $W(s)$ is also smooth. Using Lemma 27, the walk operator $W(s)$ also satisfies the gap condition as $H(s)$ does by choosing a sufficiently small (but independent of $\epsilon$) time step size $h$. Finally, according to the chain rule, each term in the derivatives of $W(s)$ contains the derivatives of $f(s)$ as a multiplicative factor, which implies that $W(s)$ also satisfies the boundary cancellation condition. Therefore, our claim directly follows from 10. ∎

As a side product, we can obtain an improved error bound on the numerical discretization errors as well. The idea is to use the triangle inequality to bound the numerical discretization error by the sum of the continuous and discrete adiabatic errors. Here we state and prove an error bound for the first-order product formula. Similar results for other numerical methods can be established in the same way.

Corollary 13**.**

Let $|\psi(T)\rangle$ be the solution of Eq. 1 at the final time, and $|\widetilde{\phi}\rangle$ be the numerical solution obtained by the first-order product formula with step size $\mathcal{O}(1)$. Suppose that 10 is true, both $\|H_{0}\|$ and $\|H_{1}\|$ are bounded by $1$, the Hamiltonian $(1-s)H_{0}+sH_{1}$ satisfies the gap condition, the spectrum of interest only consists of one simple eigenvalue, and the smooth interpolation function $f(s)$ has $f^{(k)}(0)=f^{(k)}(1)=0$ for any $k\geq 1$. Then, for any $k\geq 1$, we have

[TABLE]

Proof.

First, according to the continuous high-order adiabatic theorem, which we detail in Lemma 19 in Appendix A, there exists an eigenstate $|u\rangle$ of $H_{1}$ corresponding to the eigenvalue of interest such that for any $k$,

[TABLE]

Similarly, for the numerical solution, according to Theorem 12, there exists an eigenstate $|v\rangle$ of $e^{-ihH_{1}}$ corresponding to the eigenvalue of interest such that for any $k$,

[TABLE]

Notice that $|u\rangle$ and $|v\rangle$, if not exactly the same, only differ by a phase factor, since we assume that $\|H_{1}\|\leq 1$ and the spectrum of interest only consists of one simple eigenvalue. Therefore, the triangle inequality gives

[TABLE]

∎

We compare Corollary 13 with the standard Trotter error bound in two aspects. First, standard Trotter error bounds are typically in terms of operator norm, and thus simultaneously bound the fidelity error and the phase error. However, our result in Corollary 13 only bounds the difference in the density matrix, and there is nothing we can conclude from Corollary 13 about the error within the phase. Nevertheless, in the scenario of near adiabatic evolution, the phase factor is usually not of importance and our bound on the density matrices is usually sufficient. Second, in terms of the density matrices, our error bound significantly improves the standard ones. Our Trotter error bound for near adiabatic evolution is $\mathcal{O}(T^{-k})$ for an arbitrarily large $k$, while the standard error bound of first-order Trotter with time step size $\mathcal{O}(1)$ is only $\mathcal{O}(1)$. We remark that such an improvement requires additional conditions, including the regularity conditions and the Hamiltonian gap conditions specified in Corollary 13.

V Gap conditions of product formula and discrete adiabatic evolution for gapless Hamiltonian

Our analysis requires the gap condition of the local numerical propagator due to the usage of the discrete adiabatic theorems. For product formula, if only the gap in the Hamiltonian $H(s)$ is assumed rather than the product walk operator, we can reduce the time step size $h$ to make sure that the product walk operator also has a spectral gap (as shown in Lemma 27 and Corollary 9).

However, when the time step size is large (that is, as large as $\mathcal{O}(1)$), the relation between the gap of the product walk operator and the gap of the Hamiltonian is not straightforward yet since the perturbation theory does not directly apply. We show by toy numerical examples that, in the worst case, even if $\|H_{0}\|$ and $\|H_{1}\|$ are bounded, the gap condition of the Hamiltonian does not necessarily imply that of the first-order Trotter operator, and vice versa.

The first example is that the gap of $H(s)$ remains constant while the gap of $W(s)$ closes. We consider a $4$-dimensional system, and $H_{0}$ and $H_{1}$ are constructed as follows. Let $Q$ be the eigenbasis of the matrix

[TABLE]

$D=\text{diag}(-0.5,-0.5+\epsilon,0.2,0.6)$ for a small parameter $\epsilon\geq 0$, and define $\widetilde{W}=Qe^{-iD}Q^{\dagger}$. Then we choose $H_{1}=\text{diag}(-1,-0.6,0,1)$, $H_{0}=2i\log(e^{iH_{1}/2}\widetilde{W})$, and the scheduling function is linear (i.e., $f(s)=s$). By definition, $W(0.5)$ is exactly $\widetilde{W}$, so we can control the gap of $W(s)$ by tuning $\epsilon$. In particular, the gap vanishes when we choose $\epsilon=0$. As an illustration, the spectrum of $H(s)$ and $W(s)$ are shown in Fig. 1. Notice that the gap of $W(s)$ is measured using the distance in $\mathbb{S}^{1}$ (i.e., the spectral gap of $i\log(W(s))$). At a first glance of Fig. 1, the spectrum of $H(s)$ and $i\log(W(s))$ are surprisingly very similar. However, a careful computation yields small differences. We numerically compute the spectral gap between ground and first excited states, and results with different $\epsilon$ are shown in Table 2. As $\epsilon$ decreases, the gap of $W(s)$ decreases, but the gap of $H(s)$ remains at a constant level (although as small as $\sim 10^{-3}$) when $\epsilon$ is small enough. In the extreme case where $\epsilon=0$, the gap of $W(s)$ becomes [math] as expected (numerical gap is $\sim 10^{-16}$, reaching machine epsilon), while the gap of $H(s)$ is still $\sim 10^{-3}$.

Similarly, we can construct a second example where the gap of $W(s)$ remains constant while the gap of $H(s)$ closes. Let $Q$ and $D=\text{diag}(-0.5,-0.5+\epsilon,0.2,0.6)$ be the same as in the first example. Define $\widetilde{H}=QDQ^{\dagger}$, $H_{1}=\text{diag}(-1,-0.6,0,1)$, $H_{0}=2\widetilde{H}-H_{1}$, and the scheduling function $f(s)=s$. Again by definition, $H(0.5)$ is exactly $\widetilde{H}$ and thus can be gapless when $\epsilon=0$. Table 3 shows the spectral gaps with different $\epsilon$, and we can observe that, as $\epsilon$ decreases, the gap of $H(s)$ vanishes while the gap of $W(s)$ remains at a constant level.

The second example where the gap of $H(s)$ closes while the gap of $W(s)$ of the first-order Trotter operator remains constant suggests the possibility of efficient adiabatic optimization for systems with exponentially small gap or even gapless systems. Specifically, for gapless systems, we can try first-order Trotter with time step size $1$. Thanks to the discrete adiabatic theorems, such an approach will work and approximate the desired eigenstate if the product walk operators satisfy the gap condition, which is possible even for gapless systems (and we have observed this in our second example). Meanwhile, decreasing the time step sizes might worsen the computation because it makes the numerical evolution closer to the continuous dynamics which is known to be away from the desired eigenstate.

We illustrate this phenomenon using our second toy model again. Here we directly choose $\epsilon=0$, so the Hamiltonian $H(s)$ is gapless and the walk operator $W(s)$ of first-order Trotter with time step size $1$ still satisfies the gap condition.  Fig. 2 shows the spectrums of $H(s)$ and $W(s)$, and confirms that $H(s)$ is gapless and $W(s)$ is gapped (though the gap is indeed very small). We choose time $T=10^{k}$ for $k=3,4,5,6,7$, and implement first-order Trotter method with time step sizes $h=2^{-k}$ for $k=0,1,2,3,4,5$. In Fig. 3, we plot the fidelity $|\innerproduct{\widetilde{\phi}}{e_{j}}|$ for $j=0,1$ where $|\widetilde{\phi}\rangle$ is our final actual state and $|e_{j}\rangle$ is the $j$-th eigenstate of the final Hamiltonian $H_{1}$. When the time step size is fixed to be $1$, the actual state will better approximate the ground state of $H_{1}$ as $T$ becomes larger. This can be explained by the discrete adiabatic theorems and the analysis we establish in this work. However, if we decrease the time step size, then the numerical solution will go away from the ground state and tend to the first excited state instead. This is because, as the time step size decreases, the numerical solution will be closer to the solution of the continuous dynamics governed by $H(s)$, which is gapless and does not even approximately follow the ground eigenpath. The observation that it goes to the first excited state can be explained according to the adiabatic theorem without gap condition Avron and Elgart (1999) since the ground state of the initial Hamiltonian $H_{0}$ and the first excited state of the final Hamiltonian $H_{1}$ is smoothly connected. In summary, this toy example shows that for gapless Hamiltonian, continuous adiabatic evolution does not approximate the ground state but the Trotterization with $\mathcal{O}(1)$ time step size may still work as long as the product walk operator satisfies the gap condition.

VI Optimized gap dependence with application to Grover search and comparison with QAOA

We have shown that the time step size in discretizing near-adiabatic dynamics can be larger as previously expected by viewing the numerical integrator as the set of discrete adiabatic walk operators. In our analysis, we mainly focus on the improved scalings in $\epsilon$, and we simply bound the time-dependent spectral gap by the minimal gap, resulting in a generally super-linear dependence on the gap. However, the gap dependence can be improved to be linear if we are given sufficient information of the simultaneous gap and carefully choose the scheduling function $f(s)$ according to the size of the gap. This has been observed in several applications of the continuous adiabatic evolution as well as using the first-order exponential integrator to construct optimal linear system solver. Here we demonstrate by the adiabatic Grover search example that the linear gap dependence can also be achieved in first-order Trotter method as well.

We consider the unstructured search problem in $N$ dimensional Hilbert space. Let $\left\{|j\rangle\right\}_{j=0}^{N-1}$ denote the computational basis, and $\mathcal{M}\subset\left\{|j\rangle\right\}_{j=0}^{N-1}$ is the set of the marked states with $|\mathcal{M}|=M$. We assume $N\geq 2M$. Starting with the uniform superposition of the computational basis, our goal is to find the uniform superposition of all marked states. Here we do not assume a priori knowledge on $M$ since estimating the bound of $M$ can be costly – instead we are interested in designing robust algorithm that does not use the information of $M$ in implementing the algorithm and automatically becomes more accurate when $M$ is larger.

Adiabatic formulation Albash and Lidar (2018) takes the initial Hamiltonian to be

[TABLE]

where $|u\rangle=\frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}|j\rangle$ is the uniform superposition, and the final Hamiltonian to be

[TABLE]

Notice that both $H_{0}$ and $H_{1}$ only effectively act on a two-dimensional subspace. The orthonormal basis of this subspace is $\left\{|e_{0}\rangle,|e_{1}\rangle\right\}$ where

[TABLE]

and the corresponding matrix representations of $H_{0}$ and $H_{1}$ are

[TABLE]

Furthermore, $|u\rangle$ is the ground state of $H_{0}$ and our target state $|e_{0}\rangle$ is the ground state of $H_{1}$. Let $H(s)=(1-f(s))H_{0}+f(s)H_{1}$ where $f(s)$ is the scheduling function such that $f(0)=0$ and $f(1)=1$. We may solve the time evolution governed by $H(s)$ with the initial state $|u\rangle$ to approximate the target state. The corresponding digital simulation algorithm, using first-order Trotter method with time step size $1$ (so $T=T_{d}$ in this case), goes as

[TABLE]

For time-dependent Hamiltonian $H(s)=(1-f(s))H_{0}+f(s)H_{1}$, straightforward computations yield that two relevant eigenvalues which determine the diabatic errors are Albash and Lidar (2018)

[TABLE]

The gap of $H(s)$ is

[TABLE]

However, in the discrete evolution, the gap of the walk operator $W(s)=\exp\left(-if(s)H_{1}\right)\exp\left(-i(1-f(s))H_{0}\right)$ is what determines the scaling of $T$. Since $W(s)$ is represented by a $2$-dimensional matrix, it is possible to explicitly compute its eigenvalues and spectral gap (see Appendix G). However, the explicit expression of the gap turns out to be complicated for us to estimate the adiabatic error using Lemma 2. To facilitate further computations, we alternatively relate the gap of $W(s)$ to that of $H(s)$ in the following lemma, for which the proof can be found in Appendix G.

Lemma 14**.**

Consider the adiabatic Grover search problem. Let $H(s)$ be the Hamiltonian with the gap $\Delta_{H}(s;N,M,f)$ specified in Eq. 74, and $\Delta_{W}(s;N,M,f)$ denote the gap of the walk operator $W(s)$. Then

[TABLE]

The minimum gap appears when $s=f^{-1}(1/2)$ and $\sim\sqrt{M/N}$. As a result, a linear interpolation ($f(s)=s$) leads to the non-optimal scaling $T\geq N/M$ which is at least quadratically worse than Grover’s algorithm. To recover the Grover speedup, we need to choose a carefully designed scheduling function such that it slows down when the gap closes.

We consider two options. The first one is to use the scheduling function designed for $M=1$, so the scheduling function is independent of general values of $M$. The intuitive reason why this may work is that the case $M=1$ is the hardest one with smallest gap near $s=f^{-1}(1/2)$, so the resulting scheduling function is also sufficiently slow near the smallest gap for general $M>1$. One example is defined through the differential equation

[TABLE]

where the normalization constant $d_{N,p}$ is chosen such that $f(1)=1$ and has the expression

[TABLE]

Such a scheduling function was originally proposed in Roland and Cerf (2002) with $p=2$ and generalized in Jansen et al. (2007) with $1<p<2$ for continuous AQC. The work Dalzell et al. (2017) considers the time discretization error with $p=2$ and shows that the time step size can be as large as $1$ via efficient Hamiltonian approach. For $1\leq p<2$, we give the error bound in the following theorem, and the proof can be found in Appendix H.

Theorem 15**.**

Consider solving the Grover search problem in $N$ dimension with $M$ marked states using AQC discretized by first-order product formula with time step size $1$. Suppose that $M$ is unknown.

Suppose that $T\geq\mathcal{O}(\log(N))$ and the scheduling function is defined in Eq. 76 with $p=1$. Then the overall error is bounded by

[TABLE] 2. 2.

Suppose that $T\geq\mathcal{O}(N^{\frac{p-1}{2}})$ and the scheduling function is defined in Eq. 76 with $1<p<2$. Then the overall error is bounded by

[TABLE]

Theorem 15 implies that, to achieve an $\epsilon$ error, the smallest number of Trotter steps required is $\mathcal{O}\left(\sqrt{\frac{N}{M}}\frac{\log(N)}{\epsilon}\right)$ for $p=1$, and $\mathcal{O}\left(\frac{\sqrt{N}}{M^{1-p/2}}\frac{1}{\epsilon}\right)$ for $1<p<2$. The discrete adiabatic evolution with $p=1$ is nearly optimal (up to a logarithmic factor) in both $N$ and $M$, and that with $1<p<2$ is optimal in $N$ but not optimal in $M$. Notice that these critical points are only for comparison with other algorithms and not for designing algorithms since we do not know $M$ a priori – instead we will use sufficiently large $T$ and the error becomes automatically smaller if $T$ is too large.

Now we consider another scheduling function satisfying boundary cancellation condition. We also want it to slow down at $s=1/2$ to preserve the Grover speedup. We start with the function

[TABLE]

where

[TABLE]

One can show that $g(0)=0$, $g(1)=1$, and $g^{(k)}(0)=g^{(k)}(1)=0$ for all $k\geq 1$. Consider

[TABLE]

which is illustrated in Fig. 4.

Intuitively, we symmetrically connect two $g$’s and the resulting function $f(s)$ is smooth, satisfies the boundary cancellation and becomes slow in the middle. With this scheduling function Eq. 82, we can bound the error as follows, and the proof can be found in Appendix I.

Theorem 16**.**

Consider solving the Grover search problem in $N$ dimension with $M$ marked states using AQC with the scheduling function Eq. 82, discretized by first-order product formula with time step size $1$. Suppose that 10 is true, and $T\geq\mathcal{O}(\sqrt{N/M})$. Then, for any integer $k\geq 1$, the overall error can be bounded by

[TABLE]

Here $C$ is an absolute constant, and $C_{k}(N,M)$ is a constant that depend on $k,N,M$.

Theorem 16 can be understood from two angles. On the one hand, for constant level of error $\epsilon$, we use the first error bound and the number of Trotter steps required scales $\widetilde{\mathcal{O}}(\sqrt{N/M})$, achieving the Grover scaling up to logarithmic factors. On the other hand, for fixed $N$ and $M$, we use the second error bound to conclude that the convergence in $T$ is super-polynomial and $T$ scales $\widetilde{\mathcal{O}}(1/\epsilon^{o(1)})$. Notice that the scheduling function Eq. 82 is independent of both $M$ and $N$.

Finally we compare our results with QAOA. Our results imply an approach to construct the angles in QAOA via adiabatic evolution for Grover search problem. Specifically, QAOA for Grover search takes the ansatz

[TABLE]

where the parameters $(\beta_{j},\gamma_{j})$ are optimized to minimize $\innerproduct{\psi}{H_{1}|\psi}$. According to our discrete adiabatic result, the first-order Trotter method with time step 1 outputs

[TABLE]

and $|\widetilde{\phi}\rangle$ is a good approximation of the ground state of $H_{1}$. This implies a set of near-optimal QAOA parameters

[TABLE]

where $f(s)$ is the scheduling function in Eq. 76 with $1<p<2$. When $M=1$, the depth and time complexity of QAOA are both $p=T\sim\sqrt{N}$ which already matches the Grover lower bound. Therefore, our result shows that discrete adiabatic evolution can asymptotically match the performance of optimized QAOA for Grover search problem, and the speedup of QAOA is at most a constant factor.

VII Conclusions

In this work, we establish improved error bounds and complexity analysis for the first-order exponential integrator and product formula when applied to simulating near-adiabatic quantum dynamics. The key step in our analysis is to view the numerical integrator as the set of slowly varying unitary walk operators, and to directly bound the distance between the actual output state and the ideal eigenstate. Compared to the standard error analysis which combines the continuous adiabatic theorem and the bounds for time discretization error in the general setting, our method provides a much more direct point of view and thus yields improvements. In the first-order exponential integrator or product formula of any order, for accurate simulation, it suffices to choose a uniform time step size independent of the error $\epsilon$ and the evolution time $T$. This greatly improves the previous estimate which scales as $\epsilon^{1/p}/T^{1/p}$ for a $p$-th order product formula. Remarkably, when the Hamiltonian is bounded by $1$, a first-order method unit time step size would suffice. Furthermore, under the boundary cancellation condition, even first-order methods with uniform time step size can achieve exponential convergence.

The key technical tool in our analysis is the discrete adiabatic theorems. In particular, under the boundary cancellation condition, we use a high-order discrete adiabatic theorem to show exponential convergence for first-order methods. However, there are two open questions in the high-order discrete adiabatic theorem. First, as discussed in Section F.2, although the conclusion of the high-order discrete adiabatic theorem is numerically correct, the proof presented in Dranov et al. (1998) seems not valid due to a missing step. It is our future work to fix this gap and establish a rigorous proof for the high-order discrete adiabatic theorem. Another drawback of the current high-order discrete adiabatic theorem is that it only analyzes the scaling of errors in terms of the evolution time $T$. It would be interesting and practically useful to enhance the result with explicit dependence on the spectral gap as well, similar to what has been done for first-order discrete adiabatic theorem in Costa et al. (2022).

To apply the discrete adiabatic theorems, we require the numerical integrator to satisfy the gap condition as well. When analyzing high-order product formula, we match the gap condition for the numerical integrator with that of the time-dependent Hamiltonian $H(s)$, as long as the time step size is smaller than a threshold which only depends on the Hamiltonian $H(s)$ itself. However, when the time step size is large, there is no clear relation between the gap of the numerical integrator and the gap of the Hamiltonian, as shown in Section V via toy models. Remarkably, it is possible that, even if the Hamiltonian $H(s)$ has small gap or becomes gapless, its trotterization with large time step may regain the gap condition, indicating that adiabatic evolution with suitable discretization may achieve exponential speedup over the standard AQC framework. It is very interesting to find if there is any application of practical interest that has such feature in its gap condition.

We apply our analysis to the example of adiabatic Grover search problem. For this problem, we find that the QAOA approach does not yield better asymptotic scaling compared to the Trotterized adiabatic approach. Relation and comparison between AQC and QAOA for eigenstate preparation and optimization is a long-standing topic, and a thorough and comprehensive understanding may significantly advance the theoretical study on QAOA or even more general variational quantum algorithms. Our work here is the first step to bridge the connection, in the scenario where the Hamiltonian is bounded and well gapped. The next step would be to further explore AQC (and especially discrete AQC) and QAOA for gapped Hamiltonian but with large spectral norm, and then even for gapless Hamiltonian.

Acknowledgements

The authors thank Lucas Brady, Jacob Bringewatt, Andrew Childs, Alexey Gorshkov, Michael Jarret, and Lin Lin for helpful discussions. DA acknowledges funding from Innovation Program for Quantum Science and Technology via Project 2024ZD0301900, the support by The Fundamental Research Funds for the Central Universities, Peking University, and the support by the Department of Defense through the Hartree Postdoctoral Fellowship at QuICS. DWB is supported by Australian Research Council Discovery Projects DP210101367 and DP220101602.

Appendix A Continuous adiabatic theorems

Here we briefly review the previous results about continuous adiabatic theorems, which aim to bound the distance between the dynamics in Eq. 1 and the ideal eigenspace.

For a general interpolation function $f(s)$, the continuous adiabatic error can be bounded linearly in the inverse evolution time. Here we present a result from Ref. Jansen et al. (2007). Let $U(t)$ denote the exact evolution operator of Eq. 1, and $P(s)$ denote the spectral projection onto the eigenspace of $H(s)$ of interest. Suppose that $\Delta(s)$ is the spectral gap of the Hamiltonian $H(s)$.

Lemma 17** (Continuous adiabatic theorem).**

Suppose that $H(s)$ is second-order continuously differentiable. Let $\Delta(s)$ denote the spectral gap of $H(s)$. Then for any $t\in[0,T]$ and $s=t/T$, we have

[TABLE]

for a constant $C>0$ which only depends on the number of the eigenstates contained in $P(s)$.

When the Hamiltonian $H(s)$ satisfies the boundary cancellation condition, namely that the derivatives of $H(s)$ of any order vanish at the boundary, convergence of the final adiabatic error can be improved to be exponential in time. Such a high-order adiabatic theorem was rigorously proved in Nenciu (1993) without explicit gap dependence and in Ge et al. (2016) with cubic inverse gap dependence. Based on their approaches, a recent work An (2021) improved the gap dependence to be quadratic. Here we present the main result in An (2021).

Definition 18** (Gevrey class).**

A function $g(s)$ defined on $[0,1]$ is in the Gevrey class $G^{\alpha}$ for $\alpha>0$ if there exist constants $C,D>0$ such that for all $k\geq 0$,

[TABLE]

Here $\|\cdot\|$ represents the absolute value if $g$ is scalar-valued and the spectral norm if $g$ is operator-valued.

Lemma 19** (Continuous adiabatic theorem with boundary cancellation (An, 2021, Theorem 7)).**

Suppose that $H(s)$ is in the Gevrey class $G^{\alpha}$ for $\alpha>0$ and assume that $H^{(k)}(0)=H^{(k)}(1)=0$ for all $k\geq 1$. Let $\Delta(s)$ denote a lower bound of the spectral gap of $H(s)$ and let $\Delta_{*}=\min_{s\in[0,1]}\Delta(s)$. Then we have

[TABLE]

where $c$ and $c^{\prime}$ are positive constants only depending on $C$, $D$ (specified according to Definition 18) and $\alpha$.

Both Lemma 17 and Lemma 19 bound the leakage of the dynamics out of the desired eigenspace. We remark that, up to a constant factor $2$, it can also serve as the upper bound of the distance between the dynamics and an ideal eigenstate $P(s)|\psi(t)\rangle/\|P(s)|\psi(t)\rangle\|$ using a linear algebra result that

[TABLE]

Appendix B Multistep gap condition

The discrete adiabatic theorem in Lemma 2 (Section III.2) needs the multistep gap condition, i.e., the gap is assumed between successive steps, which is a further subtlety. Here we show that, if the number of steps $T_{d}$ is large enough and each $W(s)$ has a uniformly bounded spectral gap, then the multistep gap condition can be guaranteed.

Lemma 20**.**

Let $W(s)$ be the walk operator, $T_{d}$ be the number of discrete walk steps, and $c_{1}(s)$ be the function such that $\|DW(s)\|\leq c_{1}(s)/T_{d}$. Suppose that the gap between two groups $\sigma_{P}(s)$ and $\sigma_{Q}(s)$ of the eigenvalues of $W(s)$ is bounded from below by $\Delta_{*}>0$. Then, as long as $T_{d}\geq(2\pi/\Delta_{*})\sup_{\widetilde{s}\in\{s,s+1/T_{d}\}\cap[0,1]}c_{1}(\widetilde{s})$, the multistep gap between the arcs $\sigma_{P}^{(2)}(s)$ and $\sigma_{Q}^{(2)}(s)$ is bounded from below by $\Delta_{*}/2$.

Proof.

Let $\lambda(s)$ and $\mu(s)$ denote the two eigenvalues of $W(s)$ that determine the gap, i.e., at the corresponding boundary of the arcs. For any $s^{\prime},s^{\prime\prime}\in\{s,s+1/T_{d},s+2/T_{d}\}\cap[0,1]$, according to Bhatia and Davis (1984) and the definition of $c_{1}(s)$, we have

[TABLE]

Let $\mathbb{S}^{1}$ denote the unit circle, and define $|z_{1}-z_{2}|_{\mathbb{S}^{1}}$ to be the angular distance between $z_{1},z_{2}\in\mathbb{S}^{1}$. Using $|z_{1}-z_{2}|_{\mathbb{S}^{1}}\leq\frac{\pi}{2}|z_{1}-z_{2}|$, we have

[TABLE]

We have required that $T_{d}\geq\frac{2\pi}{\Delta_{*}}\sup_{\widetilde{s}\in\{s,s+1/T_{d}\}\cap[0,1]}c_{1}(\widetilde{s})$, which gives

[TABLE]

Therefore, by the triangle inequality, we have

[TABLE]

which means that the multistep gap is at least $\Delta_{*}/2$. ∎

Appendix C Bounds for finite differences

To apply the discrete adiabatic theorem stated in Lemma 2, we need to estimate the scalings of $c_{1}(s)$ and $c_{2}(s)$. Here we show how to bound $c_{1}(s)$ and $c_{2}(s)$ by the derivatives of the walk operator $W(s)$, and derive explicit estimates for the first-order exponential integrator and the product formula.

We first show that for each $k$, $c_{k}(s)$ is closely related to the $k$-th order derivative of $W(s)$.

Lemma 21**.**

Let $W(s)$ be a set of parameterized unitary walk operators for $s\in[0,1]$, and $D^{(k)}$ represents the $k$-th order finite difference of $W(s)$ with step size $1/T_{d}$ for an integer $T_{d}$. For any positive integer $k$, if $W(s)$ is $k$-th order continuously differentiable, then we have

[TABLE]

where

[TABLE]

Proof.

For $k=1$, using the fundamental theorem of calculus, we have

[TABLE]

More generally,

[TABLE]

This is the matrix form of a standard result for scalar finite differences (Milne-Thomson, 1933, p. 10), and is easily found by induction with

[TABLE]

This expression gives, by the triangle inequality,

[TABLE]

Therefore $\|D^{(k)}W(s)\|$ is upper bounded with $c_{k}$ as in Eq. (96). ∎

Now we give explicit choices of $c_{1}(s)$ and $c_{2}(s)$ for the first-order exponential integrator and the product formula.

Lemma 22**.**

Let $W(s)=U_{\exp}(sT+h,sT)$ be the first-order exponential integrator defined in Eq. 5, where $h=T/T_{d}$ is the time step size. Suppose that both $|f^{\prime}(s)|$ and $|f^{\prime\prime}(s)|$ are uniformly bounded over $[0,1]$. Then we can choose

[TABLE]

Proof.

First, using the formula that for any operator $X(s)$,

[TABLE]

we can compute that

[TABLE]

and, according to the chain rule,

[TABLE]

Therefore, using the bound $\|H_{1}-H_{0}\|\leq(\|H_{0}\|+\|H_{1}\|)$, we have

[TABLE]

and

[TABLE]

Then we can choose $c_{1}(s)$ and $c_{2}(s)$ as stated according to Lemma 21 and Eqs. 105 and 106. ∎

Lemma 23**.**

Let

[TABLE]

be the product formula defined in Eq. 8, where $h=T/T_{d}$ is the time step size. Suppose that both $|f^{\prime}(s)|$ and $|f^{\prime\prime}(s)|$ are uniformly bounded over $[0,1]$. Then we can choose

[TABLE]

Proof.

Let

[TABLE]

and

[TABLE]

so that $W(s)=\prod_{k=0}^{K_{p}}V_{k}(s)$. By the product rule, we have

[TABLE]

and

[TABLE]

Because the exponent of each operator $V_{k,j}$, $j=0,1$, involves a single operator, they commute with their derivatives, so we can use

[TABLE]

so

[TABLE]

where $\varphi_{k}(s)=(sT+\gamma_{p,k}h)/T$ and $\phi_{k}(s)=(sT+\delta_{p,k}h)/T$. So we have $\|V_{k,0}^{\prime}(s)\|=\mathcal{O}(h\|H_{0}\|)$ and $\|V_{k,1}^{\prime}(s)\|=\mathcal{O}(h\|H_{1}\|)$, and

[TABLE]

Then, according to Lemma 21, we can choose $c_{1}(s)=\max\|W^{\prime}(s)\|=\mathcal{O}\left(h(\|H_{0}\|+\|H_{1}\|)\right)$.

Now we estimate the second-order derivative of $W(s)$. Using the product rule again, we can compute

[TABLE]

Then

[TABLE]

Notice that

[TABLE]

so we have $\|V_{k,0}^{\prime\prime}(s)\|=\mathcal{O}(h\|H_{0}\|+h^{2}\|H_{0}\|^{2})$ and $\|V_{k,1}^{\prime\prime}(s)\|=\mathcal{O}(h\|H_{1}\|+h^{2}\|H_{1}\|^{2})$. Together with $\|V_{k,0}^{\prime}(s)\|=\mathcal{O}(h\|H_{0}\|)$ and $\|V_{k,1}^{\prime}(s)\|=\mathcal{O}(h\|H_{1}\|)$, we have

[TABLE]

Then we can choose $c_{2}(s)=\mathcal{O}\left(h(\|H_{0}\|+\|H_{1}\|)+h^{2}(\|H_{0}\|+\|H_{1}\|)^{2}\right)$ according to Lemma 21. ∎

For the simplified product formula, the same bound also holds because the simplified product formula is just a special case of the general product formula.

Lemma 24**.**

Let

[TABLE]

be the simplified product formula defined in Eq. 9, where $h=T/T_{d}$ is the time step size. Suppose that both $|f^{\prime}(s)|$ and $|f^{\prime\prime}(s)|$ are uniformly bounded over $[0,1]$. Then we can choose

[TABLE]

Appendix D Numerical methods and their standard complexity estimates

In the main text, we focus on two types of numerical methods for discretizing Eq. 1, namely the first-order exponential integrator Eq. 5 and the product formula Eq. 8. Here we present their standard time discretization error bounds and complexity estimates.

D.1 Standard time discretization error bounds

For the first-order exponential integrator, its standard error bound is given as follows.

Lemma 25**.**

Let $U_{\exp}(t+h,t)$ be the first-order exponential integrator defined in Eq. 5. Suppose that $|f^{\prime}(s)|$ is uniformly bounded over $[0,1]$. Then we have $\left\|U_{\exp}(t+h,t)-\mathcal{T}\exp\left(\int_{t}^{t+h}H(\tau/T)d\tau\right)\right\|\leq C_{1,H}h^{2}$, where

[TABLE]

Proof.

For a fixed time $t$ and a step size $h$, let $\tau\in[0,h]$. The exact evolution operator from $t$ to $t+\tau$ is

[TABLE]

and the first-order exponential integrator can be written as

[TABLE]

These two operators satisfy the differential equations, respectively,

[TABLE]

and

[TABLE]

According to the variation of parameters formula (as known as Duhamel’s principle), we may write $U_{\exp}$ as a perturbation of $U$, that is,

[TABLE]

In the second line we use the fundamental theorem of calculus, and we are getting an extra $1/T$ scaling because the Hamiltonian is changing slowly on the scale of $T$ (technically, we are taking the derivative with respect to $\tau^{\prime}$ for the function $H(t/T+\tau^{\prime}/T)$, in which the argument $\tau^{\prime}$ is divided by $T$). Since $\|H^{(1)}(s)\|=|f(s)|\|H_{1}-H_{0}\|\leq\mathcal{O}(\|H_{0}\|+\|H_{1}\|)$, we have

[TABLE]

∎

We show an error bound for the Trotter-Suzuki formula in the next result, which is a direct corollary of the general analysis in Wiebe et al. (2010). We also expect a similar error bound for general product formulae beyond Trotter-Suzuki.

Lemma 26**.**

Let $U_{\text{pf},p}(t+h,t)$ denote the $p$-th order time-dependent product formula obtained by Sukuzi recursion. Suppose that the scheduling function $f(s)$ is smooth with uniformly bounded derivatives of any order, and $\|H_{0}\|+\|H_{1}\|\geq 1$. Then, for any $p\geq 1$ and $T\geq 1$, we have

[TABLE]

where $C_{p}$ is a constant that depends on $p$ and the derivatives of $f(s)$.

Proof.

According to (Wiebe et al., 2010, Theorem 3), there exists a $p$-dependent constant $\widetilde{C}_{p}$ such that

[TABLE]

Here $\Lambda_{p}$ is a constant such that

[TABLE]

By $H(s)=(1-f(s))H_{0}+f(s)H_{1}$ and the triangle inequality, it suffices to choose

[TABLE]

Let $F_{m}$ be the upper bound of $|f^{(m)}(s)|$ for all $s\in[0,1]$, then we have

[TABLE]

So we can choose $\Lambda_{p}=\max_{0\leq m\leq p}F_{m}^{1/(m+1)}\left(\left\|H_{0}\right\|+\left\|H_{1}\right\|\right)$, and the claimed error bound can be obtained by defining $C_{p}=\max_{0\leq m\leq p}F_{m}^{(p+1)/(m+1)}\widetilde{C}_{p}$. ∎

D.2 Standard complexity estimates

Now we show the standard approach to analyze the complexity of using discretized AQC to prepare an eigenstate. Specifically, let $|\phi\rangle$ denote the exact target eigenstate and

[TABLE]

where $U_{\text{num}}$ is a $p$-th order local numerical propagator in the sense that the local truncation error can be bounded by

[TABLE]

where $C_{p,H}$ is a constant that might depend on $p$ and certain norms of $H$, $H_{0}$, $H_{1}$, but is independent of the time step size $h$. Notice that in the first-order exponential integrator, this $C$ parameter also depends on $T$. We would like to study how many local numerical propagators, namely $M$, scale if we want $\||\widetilde{\phi}\rangle\langle\widetilde{\phi}|-|\phi\rangle\langle\phi|\|\leq\epsilon$ for some $0<\epsilon<1$.

A standard approach for estimating the final eigenstate error $\varepsilon_{\text{Eig}}$ is by using $\varepsilon_{\text{AQC}}$ as the upper bound from the continuous adiabatic theorem plus the error from the time discretization $\varepsilon_{\text{Disc}}$ :

[TABLE]

To bound the overall approximation error by $\mathcal{O}(\epsilon)$, it suffices to bound both the continuous adiabatic error and the time discretization error by $\mathcal{O}(\epsilon)$. In order to bound the time discretization error, which is $C_{p,H}Th^{p}$ for a $p$-th order method, we can choose

[TABLE]

Therefore the overall number of queries to $U_{\text{num}}$ can be bounded by

[TABLE]

and we only need to further study the scaling of the evolution time $T$, which needs to be chosen large enough such that the adiabatic error is bounded as well. For the general second-order continuously differentiable interpolation function $f(s)$ which is independent of the Hamiltonians, Lemma 17 implies that the adiabatic error is bounded by

[TABLE]

then the evolution time should be chosen as

[TABLE]

The overall number of steps is thus

[TABLE]

Now we consider the boundary cancellation condition. Here we further assume that both $\|H_{0}\|$ and $\|H_{1}\|$ are bounded by $1$. We can use Lemma 19 to bound the adiabatic error, which implies the choice of evolution time to be

[TABLE]

and thus the overall number of steps to be

[TABLE]

Appendix E Complexity estimate of the high-order simplified product formulae without gap condition on the walk operator

Here we give a proof of Corollary 9. The first step is to establish a connection between the gap of the Trotter walk operator and the Hamiltonian. We show this connection in the following lemma.

Lemma 27**.**

For $s\in[0,1-h/T]$ and a scheduling function $f(s)$, let $H(s)=(1-f(s))H_{0}+f(s)H_{1}$ be the interpolating Hamiltonian and $U_{\text{spf},p}(sT+h,sT)$ be the simplified product formula as in Eq. 9. Let $\Delta_{H}(s)$ and $\Delta_{U,p}(s)$ denote the spectral gap of $H(s)$ and $U_{\text{spf},p}$, respectively. Then, for any $h\leq 1/(\|H_{0}\|+\|H_{1}\|)$ and $p\geq 1$, there exist $p$-dependent constants $c_{p}$ such that

[TABLE]

Proof.

Notice that the spectral gaps of $e^{-ihH(s)}$ and $hH(s)$ are of the same value for all $h\leq 1/(\|H_{0}\|+\|H_{1}\|)$. So it suffices to investigate the spectral gaps of $U_{\text{spf},p}(sT+h,sT)$ and $e^{-ihH(s)}$. For every fixed $s$, we may bound the difference between $U_{\text{spf},p}(sT+h,sT)$ and $e^{-ihH(s)}$ by using the error bound of the time-independent product formula in (Childs et al., 2021, Theorem 11), which gives

[TABLE]

for a $p$-dependent constant $\widetilde{c}_{p}$.

According to the eigenvalue perturbation theorem of Ref. Bhatia and Davis (1984); Elsner and He (1993), there exists an ordering $\lambda_{j}(s)$ and $\widetilde{\lambda}_{j}(s)$ of the eigenvalues of $U_{\text{spf},p}(sT+h,sT)$ and $e^{-iH(s)}$, respectively, such that

[TABLE]

Let $\lambda_{0}(s)$ and $\lambda_{1}(s)$ denote the two eigenvalues that determine the gap, then the gap of $U_{\text{spf},p}(sT+h,sT)$ can be bounded by

[TABLE]

and

[TABLE]

In the second line of both estimates, we use $\theta\leq\frac{\pi}{2}\sin\theta$ for all $\theta\in[0,\pi/2]$. The proof is completed by choosing $c_{p}=\pi\widetilde{c}_{p}$. ∎

Now we are ready to prove Corollary 9.

Proof of Corollary 9.

The sketch of the proof is as follows. We will first give a lower bound for the fixed-time spectral gap of the walk operator $W(s)$, then give a lower bound for the multistep gap by verifying the conditions in Lemma 20, and apply the discrete adiabatic theorem to bound the overall error between the actual and ideal evolution, from which we can infer the choices of $T$ and $h$.

Let $T_{d}=T/h$, and we assume $T_{d}$ is an integer. We first bound the fixed-time gap of $W(s)$ defined by Eq. 9. According to the choice of $h$ with proper constant factor, we can make

[TABLE]

where $c_{p}$ is the constant factor in Lemma 27. Then, using Lemma 27, the gap of $W(s)$ is bounded from below by the gap of $\frac{1}{2}h\Delta(s)$, which has a further lower bound $\frac{1}{2}h\Delta_{*}$.

For the multistep gap, we can use Lemma 20, and we need to verify that $T_{d}$ is sufficiently large. Specifically, suppose that our choice of $T$ is

[TABLE]

for a constant $C>0$. Then we have

[TABLE]

Notice that $\Delta_{*}\leq 2\|H(s)\|\leq 2(\|H_{0}\|+\|H_{1}\|)$, so we have

[TABLE]

According to Lemma 24, we have $c_{1}(s)=\mathcal{O}(h(\|H_{0}\|+\|H_{1}\|))=\mathcal{O}(1)$, so we can choose a sufficiently large constant $C$ such that $T_{d}\geq\frac{2\pi}{h\Delta_{*}/2}\sup c_{1}(s)$. Then Lemma 20 ensures that the multistep gap of $W(s)$ is bounded from below by $h\Delta_{*}/4$.

Therefore, the error bound in Theorem 5 becomes

[TABLE]

and the desired choices of $T$ and $h$ directly follow from this error bound and the conditions for the multistep gap. ∎

Appendix F High-order discrete adiabatic theorem

Here we provide more discussion on the high-order discrete adiabatic theorem, which is claimed in Dranov et al. (1998) and we state as 10 in the main text. First, in the statement of 10, we adapt a slightly different statement compared to the original version in Dranov et al. (1998) for the purpose of clearer presentation. We will explain how to get 10 from there. More importantly, we will then point out the missing steps in the analysis of Dranov et al. (1998) and explain why these missing steps imply that the proof in Dranov et al. (1998) is not valid. Meanwhile, we also show a numerical result verifying that, despite the flawed proof, the claim of the high-order discrete adiabatic theorem is very likely to be correct.

F.1 A slight variant

We show how to derive 10 from the original result in Dranov et al. (1998). We need to introduce some operators before proceeding. Let $P(s)$ and $Q(s)$ denote the spectral projections onto $\sigma_{P}(s)$ and $\sigma_{Q}(s)$, respectively. Define

[TABLE]

The following is the original result from Dranov et al. (1998).

Lemma 28** (Higher-order discrete adiabatic theorem, (Dranov et al., 1998, Theorem 3)).**

Suppose that $V(s)-I$ is supported in $(0,s_{\star})$. Then, for any integer $n$ such that $n/T_{d}\notin(0,s_{\star})$ and for any positive integer $k$, we have

[TABLE]

where $C_{k}$ is a constant independent of $n$ and $T_{d}$.

Notice that Lemma 28 introduces a new parameter $s_{\star}$ which could be larger than $1$. Such an extension allows us to bound the discrete adiabatic error when the number of the walk operators is larger than $T$. Now we prove 10 from Lemma 28.

Proof of 10.

The idea is to smoothly extend $W(s)$ to a longer time interval $[0,3]$ to make it satisfy the condition in Lemma 28. Specifically, we consider a walk operator on the time interval $[0,3]$ and force it to be unchanged only on $[1,2]$ to make sure that it satisfies the boundary cancellation condition in the discrete sense as well. On the time interval $[0,3]$, we define

[TABLE]

The operators associated with $\widetilde{W}(s)$ are denoted using the same letter with an upper tilde (for example $\widetilde{P}(s)$ and $\widetilde{Q}(s)$ denote the spectral projections of $\widetilde{W}(s)$). Since $\widetilde{W}^{(k)}(1)=W^{(k)}(0)=0$ and $\widetilde{W}^{(k)}(2)=W^{(k)}(1)=0$ for all $k\geq 1$, the operator $\widetilde{W}(s)$ is a smooth extension of $W(s)$.

Notice that $\widetilde{W}(s)$ and thus $\widetilde{P}(s)$ and $\widetilde{Q}(s)$ remain unchanged on $[0,1]\cup[2,3]$. Then, for any $s\in[0,1-1/T_{d}]\cup[2,3-1/T_{d}]$, we must have

[TABLE]

so $\widetilde{v}(s)=\sqrt{\widetilde{S}(s)\widetilde{S}(s)^{\dagger}}=I$ and $\widetilde{V}(s)=\widetilde{v}(s)^{-1}\widetilde{S}(s)=I$. This means that $\widetilde{V}(s)-I$ is supported in $(0,2)$ and satisfies the assumption of Lemma 28. As a result, we have

[TABLE]

where $C_{k}$ is a constant independent of $T_{d}$.

We now show how we may obtain a bound for the desired quantity from Eq. 161. First, we have

[TABLE]

where we use Eq. 14 in the second equality. Then, since $|\psi\rangle$ is within the eigenspace corresponding to $\sigma_{\widetilde{P}}(0)$, we have

[TABLE]

where $C_{k}$ is the constant in Lemma 28. Finally, notice that, by definition, $\widetilde{P}(2)=P(1)$, and $|\psi\rangle$ is an eigenstate of $\widetilde{W}(s)$ for all $s\in[0,1]$. Suppose the corresponding eigenvalue is $e^{i\theta}$ for a real number $\theta$, then

[TABLE]

Therefore we have

[TABLE]

∎

F.2 Missing steps in the proof of the high-order discrete adiabatic theorem and a numerical validation

We first follow the notations in Dranov et al. (1998) to define several operators. Let $P(s)$ and $Q(s)$ be the spectral projections as before, and $V(s)$ is defined through Eq. 157. Let

[TABLE]

It can be shown that the Volterra equation holds as

[TABLE]

In Dranov et al. (1998), another set of operators $\left\{\Omega_{j}(s)\right\}$ is defined in a fixed-point-iteration fashion – let $\Omega_{0}(s)=1$ for all $s$, and

[TABLE]

By taking the summation over $j$, we can see that $\sum_{j=0}^{\infty}\Omega_{j}(s)$ is expected to converge to $\Omega(s)$.

Lemma 1 in Dranov et al. (1998) attempts to show that $Q_{0}\Omega_{j}(1)P_{0}$ is of order $\mathcal{O}(T_{d}^{-k})$ for arbitrary positive integer $k$ under the boundary cancellation condition, and the high-order discrete adiabatic theorem becomes a consequence of this lemma. This was proven using mathematical induction over $j$. However, in the proof of (Dranov et al., 1998, Lemma 1), the induction assumption (Eq. 40 in Dranov et al. (1998)) was taken to be $Q_{0}\Omega_{j}(n/T_{d})P_{0}=\mathcal{O}(T_{d}^{-k})$ for all $0\leq n\leq T_{d}$, which is much stronger than the original claim which only holds when $n=0$ and $T_{d}$. Furthermore, intuitively, the operator $Q_{0}\Omega_{j}(n/T_{d})P_{0}$ should not scale $\mathcal{O}(T_{d}^{-k})$ inside the time interval since the boundary cancellation condition is expected to accelerate the convergence only near the boundary.

We verify our statements numerically using a $4$-level model. We consider the walk operator to be the first-order exponential integrator with time step size $1$, i.e., $W(s)=e^{-iH(s)}$, where $H(s)=(1-f(s))H_{0}+f(s)H_{1}$. We choose the Hamiltonians $H_{j}=Q_{j}^{\dagger}D_{j}Q_{j}$ for $j=0,1$, where $D_{0}=\text{diag}(0.5,0.8,1.2,1.4)$ and $D_{1}=\text{diag}(0.3,1.0,1.5,1.9)$ are two diagonal matrices, $Q_{0}$ is the orthonormal basis of the matrix

[TABLE]

and $Q_{1}$ is the orthonormal basis of

[TABLE]

The scheduling function is the glue function $f(s)=c^{-1}\int_{0}^{s}\exp\left(-\frac{1}{s^{\prime}(1-s^{\prime})}\right)ds^{\prime}$, where $c$ is the normalization constant such that $f(1)=1$. Notice that such a scheduling function satisfies the boundary cancellation condition.

Figure 5 (a) shows the numerical scaling of $\max_{s}\|Q_{0}\Omega_{1}(s)P_{0}\|$ and $\|Q_{0}\Omega_{1}(1)P_{0}\|$ under different choices of $T_{d}$. Here we test the scaling for $\Omega_{1}$ instead of $\Omega_{0}$ because $\Omega_{0}$ is always the identity matrix so the orthogonal projection is always [math]. Our numerical results suggest that $\max_{s}\|Q_{0}\Omega_{1}(s)P_{0}\|$ is of order $\mathcal{O}(1/T_{d})$ while $\|Q_{0}\Omega_{1}(1)P_{0}\|$ decays super-polynomially, which verifies that the induction assumption used in Dranov et al. (1998) is questionable.

Although the proof presented in Dranov et al. (1998) seems flawed, we believe that the conclusion of the high-order discrete adiabatic theorem still holds true. We verify this using the same example, and our numerical results shown in Fig. 5 (b) confirms that the adiabatic error $\|Q_{0}\Omega(1)P_{0}\|$ at the boundary indeed converges super-polynomially in $1/T_{d}$.

Appendix G Spectral gap of the walk operator in adiabatic Grover search

G.1 Direct computations

Here we determine the spectral gap of the walk operator $W(s)$ for a step of Trotter as applied to a search problem in Section VI. First we diagonalise $H_{0}$ using Eq. 71 as

[TABLE]

That gives the walk operator $W(s)$ as

[TABLE]

The eigenvalues of $W(s)$ satisfy the following quadratic equation

[TABLE]

So the eigenvalues can be represented as $\lambda=e^{-\frac{i}{2}}\mu$, where $\mu$ satisfies

[TABLE]

Solving this quadratic equation yields the eigenvalues

[TABLE]

The spectral gap of $W(s)$ is

[TABLE]

G.2 Proof of Lemma 14

This Lemma lower bounds the gap of the walk operator $W(s)$ for a step of Trotter in terms of the spectral gap of the Hamiltonian. To bound this spectral gap, we may use Lemma 7 with $h=1$ to give

[TABLE]

The commutator in the basis $\left\{|e_{0}\rangle,|e_{1}\rangle\right\}$ can be computed as

[TABLE]

and

[TABLE]

Then Eq. 191 becomes

[TABLE]

Appendix H Proof of Theorem 15

This Theorem bounds the error in using the Trotter formula to solve a search problem. The majority of the proof uses the general form of the discrete adiabatic theorem in Lemma 2 to bound the diabatic errors. To use Lemma 2, we first need to bound $c_{1}$ and $c_{2}$. This can be done by straightforward calculations similar to the proof of Lemma 23 and we can choose

[TABLE]

For general $p$, by the definition of $f(s)$, we have

[TABLE]

and

[TABLE]

where

[TABLE]

Therefore we can choose

[TABLE]

Here

[TABLE]

The error bound in Lemma 2 involves the multistep gap and the maximum/minimum over several subsequent steps. Here we show that, when $T$ is sufficiently large (which can be ensured by the choice of $T$ in Theorem 15), we can ignore the maximum and minimum and consider the single-step gap condition at sacrifice of extra multiplicative constant factors. We first focus on the minimum of the multistep gap of $W$. Since the two relevant eigenvalues determining the gap have different rotational directions on the circle for any time $s$, the multistep gap of $W$ is exactly the minimum single step gap over two subsequent steps, so

[TABLE]

We bound the ratio of gaps at steps $s$ and $s^{\prime}$ as

[TABLE]

Suppose that $T\geq 12d_{N,p}$, then for any $|s-s^{\prime}|\leq 3/T$, we have

[TABLE]

This implies that

[TABLE]

Therefore we can directly replace the minimum of the multistep gap by the single-step gap at fixed time $s$ up to some constant factors. According to Eq. 201 and Eq. 205, the functions $\hat{c}_{1}(s)$ and $\hat{c}_{2}(s)$ can also be replaced, up to constant factors, by $c_{1}(s)$ and $c_{2}(s)$.

We now compute each term in Lemma 2. For simplicity of notation, we omit the dependence on $H$, $f$ and $N$ in the gap and use $\Delta(s;M)$ to denote $\Delta_{H}(s;N,M,f)$, and we ignore all constant factors in the computations. The boundary terms can be bounded as

[TABLE]

For the summations, we first bound them as

[TABLE]

and

[TABLE]

So the overall error is bounded by

[TABLE]

It remains to estimate $d_{N,p}$ and $\sum_{j=0}^{T-1}\frac{1}{T}\frac{\Delta(j/T;1)^{2p-1}}{\Delta(j/T;M)^{2}}$.

Before we proceed, it is useful to preliminarily calculate $d_{N,p}$ for general $1\leq p<2$. By Eq. 202 and the change of variables $f=\frac{1}{2}+\frac{1}{2\sqrt{N-1}}\tan\theta$, we have

[TABLE]

Here the integral in the first line comes from Eq. 198 and the condition $f(1)=1$.

H.1 $p=1$

We first consider the limit case when $p=1$. The error bound Eq. 211 becomes

[TABLE]

For $d_{N,1}$, using Eq. 212 we have

[TABLE]

For the summation, we first use the fact that $\Delta(s;M)$ is symmetric with respect to $1/2$ and monotonically decreases on $[0,1/2]$, and we can obtain

[TABLE]

Now we focus on the dominant summation. Using monotonicity of $\Delta(s;1)$, we have

[TABLE]

We change the denominator of the integrand from $\Delta(j/T;M)^{2}$ to $\Delta(s;M)^{2}$ at the cost of an extra multiplicative factor. Specifically, for any $s\in[(j-1)/T,j/T]$, using

[TABLE]

we have

[TABLE]

Suppose that $T\geq 16\log(N)$, we can conclude that

[TABLE]

Therefore Eq. 215 becomes

[TABLE]

By changing the variable from $s$ to $f$ and using $df=d_{N,1}\Delta(s;1)ds$, we have

[TABLE]

Therefore Eq. 215 becomes

[TABLE]

Plugging Eq. 214 and Eq. 222 back to Eq. 213, we can bound the total error as

[TABLE]

H.2 $1<p<2$

Now we consider $1<p<2$. Continuing with Eq. 212 and using integration by parts, we may obtain

[TABLE]

Both terms on the right hand side of Eq. 224 are bounded independent of $N$, so $\int_{0}^{\arctan\sqrt{N-1}}\cos^{p-2}\theta\,d\theta=\mathcal{O}(1)$ and thus, by Eq. 212,

[TABLE]

For the summation, with the same reasoning of Eq. 215 and Eq. 220, we have

[TABLE]

We focus on the integral. By changing the variable from $s$ to $f$ and using $df=d_{N,p}\Delta(s;1)^{p}ds$, we have

[TABLE]

By the change of variable $f=\frac{1}{2}+\frac{\sqrt{M}}{2\sqrt{N-M}}\tan\theta$, we have

[TABLE]

Notice that, for $\theta\in[0,\pi/2)$, the function $\tan^{2}\theta$ monotonically increases. So we can separate the integral into two parts as

[TABLE]

Using $1<p<2$ and thus

[TABLE]

we have

[TABLE]

Plugging this back into Eq. 226, we have

[TABLE]

Plugging Eq. 225 and Eq. 232 back to Eq. 211, we can bound the total error by

[TABLE]

Appendix I Proof of Theorem 16

In this section we give the complete proof of Theorem 16. We start with summarizing the properties of the function $f(s)$ in the following lemma.

Lemma 29**.**

Let $f(s)$ be the function defined in Eq. 82. Then

$f(s)\in C^{\infty}[0,1]$, 2. 2.

$f(0)=0$, $f(1)=1$, and $f^{(k)}(0)=f^{(k)}(1)=0$ for all $k\geq 1$, 3. 3.

for all $s\in[0,1/2]$ and $k\geq 1$, $f(s)+f(1-s)=1$ and $f^{(k)}(s)=f^{(k)}(1-s)$.

To prove Theorem 16, we will frequently use the following lemma.

Lemma 30**.**

Let $\Delta_{H}(s)$ be the gap defined in Eq. 74. Then for any $s\in[0,1/2]$ and any integer $l\geq 0$, we have

[TABLE]

where $c_{e}$ is the normalization factor defined in Eq. 81.

Proof.

We start with

[TABLE]

Let $t=2s\in[0,1]$, then

[TABLE]

where we define

[TABLE]

It suffices to bound $F(t)$. This can be done with nuanced modifications of the proof of (An and Lin, 2022, Lemma 14), but for completeness we provide the details again. Since $F(0)=F(1)=0$ and $F(t)>0$ for $t\in(0,1)$, there exists a $t_{*}\in(0,1)$ such that $F(t_{*})$ reaches its maximum. At this point, $F^{\prime}(t_{*})=0$, where $F^{\prime}(t)$ can be computed as

[TABLE]

where

[TABLE]

For $t\in[1-\frac{\sqrt{c_{e}}}{(l+1)\log(N/M)},1]$ where $\sqrt{c_{e}}\approx 0.08$, we have

[TABLE]

[TABLE]

and

[TABLE]

Then

[TABLE]

Meanwhile, for the same range of $t$,

[TABLE]

Therefore, for $t\in[1-\frac{\sqrt{c_{e}}}{(l+1)\log(N/M)},1]$, we always have

[TABLE]

Since $F^{\prime}(t_{*})=0$, we must have $G(t_{*})=0$, so Eq. 245 implies that $t_{*}\leq 1-\frac{\sqrt{c_{e}}}{(l+1)\log(N/M)}$. Meanwhile, notice that for $t\in[0,1/2)$, we have $F(t)<F(1-t)$, so $t_{*}\geq 1/2$. Therefore we obtain the range of $t_{*}$ that

[TABLE]

We now bound $F(t)$. Using $G(t_{*})=0$ again, we have

[TABLE]

Therefore

[TABLE]

∎

We are ready to present the main proof of Theorem 16. The idea is separately using Lemma 2 and 10 to obtain two independent error bounds, so the overall error is bounded by the minimum of the two.

Proof of Theorem 16.

By differentiating the product walk operator, we may choose

[TABLE]

We use the similar trick as in Appendix H to get rid of the hat and check notations and the multistep gap in Lemma 2. Specifically, for $|s-s^{\prime}|\leq 3/T$, we have

[TABLE]

Here $\max|f^{\prime}|$ is an absolute constant which is independent of $N$ and $M$, so we may choose $T\geq 6\sqrt{N/M}\max|f^{\prime}|$ and thus

[TABLE]

So we can bound $\check{\Delta}_{2,W}(s;M)$ from below by $\frac{1}{3}\Delta_{H}(s;M)$ and drop the hat notations in $c_{1}$ and $c_{2}$ in the same way as in Appendix H.

Now we bound each term in the upper bound of Lemma 2. The two boundary terms are [math] due to the boundary cancellation. For the summations, we use symmetry of $f(s)$ and Lemma 30 and obtain

[TABLE]

and

[TABLE]

Therefore Lemma 2 implies an error bound

[TABLE]

Meanwhile, since the scheduling function $f(s)$ satisfies the boundary cancellation condition, 10 implies another error bound

[TABLE]

for any integer $k\geq 1$. Here the constant $C_{k}(N,M)$ which is independent of $T$ but may depend on $k,N,M$. Therefore the overall error can be bounded by the minimum of Eq. 254 and Eq. 255. ∎