Time-dependent Hamiltonian simulation with $L^1$-norm scaling

TL;DR

This paper introduces new quantum simulation algorithms for time-dependent Hamiltonians that depend on the $L^1$ norm, offering improved efficiency over previous methods, with potential applications in quantum chemistry.

Contribution

The paper develops novel algorithms exploiting the $L^1$ norm for time-dependent Hamiltonian simulation, including a classical sampling method and a rescaling technique, improving over prior approaches.

Findings

Gate complexity scales with the $L^1$ norm of the Hamiltonian.

Rescaled Dyson-series algorithm is nearly optimal.

Sampling-based approach is more practical for near-term applications.

Abstract

The difficulty of simulating quantum dynamics depends on the norm of the Hamiltonian. When the Hamiltonian varies with time, the simulation complexity should only depend on this quantity instantaneously. We develop quantum simulation algorithms that exploit this intuition. For sparse Hamiltonian simulation, the gate complexity scales with the $L^1$ norm $\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\lVert_{\max}$, whereas the best previous results scale with $t\max_{\tau\in[0,t]}\left\lVert H(\tau)\right\lVert_{\max}$. We also show analogous results for Hamiltonians that are linear combinations of unitaries. Our approaches thus provide an improvement over previous simulation algorithms that can be substantial when the Hamiltonian varies significantly. We introduce two new techniques: a classical sampler of time-dependent Hamiltonians and a rescaling principle for the Schr\"{o}dinger equation. The rescaled Dyson-series algorithm is nearly optimal with respect to all parameters of interest, whereas the sampling-based approach is easier to realize for near-term simulation. These algorithms could potentially be applied to semi-classical simulations of scattering processes in quantum chemistry.

Tables

Table 1. Table 1 : Complexity comparison of previous algorithms (top three) and the algorithms introduced in this paper (bottom two) for simulating time-dependent Hamiltonians. Logarithmic factors are suppressed by O ~ ~ 𝑂 \widetilde{O} notation and the (non-query) gate complexities are compared. The product formula algorithm of [ 50 ] is omitted as its gate complexity scales polynomially with high-order derivatives and is not directly comparable to other algorithms in the table. The complexity of the full Monte Carlo estimation algorithm [ 42 ] is not analyzed explicitly; only its first step is compared. The fractional-query algorithm [ 6 ] does not have an explicit implementation for Hamiltonians in the LCU model, and its implementation in the SM model is streamlined by the Dyson-series approach [ 8 , 36 , 29 ] .

Algorithms	SM	LCU
Monte Carlo estimation (first step) [42]	$\tilde{O}\left({\left(d^2{\parallel H\parallel }_{\max ,\infty }t\right)}^{2}n/ϵ\right)$	$O\left({\left({\parallel \alpha \parallel }_{1,\infty }t\right)}^2{g}_e/ϵ\right)$
Fractional-query [6]	$\tilde{O}\left(d^2{\parallel H\parallel }_{\max ,\infty }tn\right)$	N/A
Dyson series [8, 36, 29]	$\tilde{O}\left(d{\parallel H\parallel }_{\max ,\infty }tn\right)$	$\tilde{O}\left({\parallel \alpha \parallel }_{\infty ,\infty }t{L}^2{g}_c\right)$
Continuous qDRIFT (Section 3.3)	$\tilde{O}\left({\left(d^2{\parallel H\parallel }_{\max ,1}\right)}^{2}n/ϵ\right)$	$O\left({\parallel \alpha \parallel }_{1,1}^2{g}_e/ϵ\right)$
Rescaled Dyson series (Section 4.2)	$\tilde{O}\left(d{\parallel H\parallel }_{\max ,1}n\right)$	$\tilde{O}\left({\parallel \alpha \parallel }_{\infty ,1}{L}^2{g}_c\right)$

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}\exp_{\mathcal{T}}\bigg{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\bigg{)}=-iH(t)\exp_{\mathcal{T}}\bigg{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\bigg{)}.$$

$$\left\lVert\exp_{\mathcal{T}}\bigg{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\bigg{)}-\exp_{\mathcal{T}}\bigg{(}-i\int_{0}^{t}\mathrm{d}\tau\,G(\tau)\bigg{)}\right\rVert_{\infty}\leq\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)-G(\tau)\right\rVert_{\infty}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{E}(t,0)=-iH(t)\mathrm{E}(t,0),\mathrm{E}(0,0)=I.$$

$$\mathrm{E}(t,0)-I=\mathrm{E}(t,0)-\mathrm{E}(0,0)={\int }_0^t\mathrm{d}\tau \frac{\mathrm{d}}{\mathrm{d}\tau }\mathrm{E}(\tau ,0)=-i{\int }_0^t\mathrm{d}\tau H(\tau )\mathrm{E}(\tau ,0).$$

$$\mathrm{E}(t,0)=I-i{\int }_0^t\mathrm{d}\tau H(\tau )\mathrm{E}(\tau ,0).$$

$$\mathrm{E}(t,0)=\mathrm{E}(t,s)\mathrm{E}(s,0).$$

$$\mathrm{E}(0,t):={\mathrm{E}}^{-1}(t,0)={\mathrm{E}}^{†}(t,0).$$

$$\mathcal{E}(t,0)(\rho ):=\mathrm{E}(t,0)\rho {\mathrm{E}}^{†}(t,0)=\mathrm{E}(t,0)\rho \mathrm{E}(0,t).$$

$${\Vert \alpha \Vert }_1:=\sum _{j=1}^L|{\alpha }_j|,{\Vert \alpha \Vert }_2:=\sqrt{\sum _{j=1}^L|{\alpha }_j{|}^2},{\Vert \alpha \Vert }_{\infty }:=\underset{j\in \{1,2,\dots ,L\}}{\max }|{\alpha }_j|.$$

$$\left\lVert A\right\rVert_{1}:=\mathrm{Tr}\big{(}\sqrt{A^{\dagger}A}\big{)},\qquad\left\lVert A\right\rVert_{2}:=\sqrt{\mathrm{Tr}\big{(}A^{\dagger}A\big{)}},\qquad\left\lVert A\right\rVert_{\infty}:=\max_{|\psi\rangle}\left\lVert A|\psi\rangle\right\rVert_{2}.$$

$${\Vert f\Vert }_1:={\int }_0^t\mathrm{d}\tau |f(\tau )|,{\Vert f\Vert }_2:=\sqrt{{\int }_0^t\mathrm{d}\tau |f(\tau ){|}^2},{\Vert f\Vert }_{\infty }:=\underset{\tau \in [0,t]}{\max }|f(\tau )|.$$

$${\Vert \alpha \Vert }_{1,1}:={\int }_0^t\mathrm{d}\tau \sum _{j=1}^L|{\alpha }_j(\tau )|,{\Vert \alpha \Vert }_{1,\infty }:=\underset{\tau \in [0,t]}{\max }\sum _{j=1}^L|{\alpha }_j(\tau )|.$$

$${\Vert A\Vert }_{2,1,\infty }:=\underset{\tau \in [0,t]}{\max }\sum _{l=1}^L{\Vert A_l(\tau )\Vert }_2.$$

$${\Vert A\Vert }_{\max }:=\underset{j,k}{\max }|A_{j,k}|.$$

$${\Vert A\Vert }_{\max }\le {\Vert A\Vert }_{\infty }.$$

$${\Vert A\Vert }_{\max ,1}:={\int }_0^t\mathrm{d}\tau {\Vert A(\tau )\Vert }_{\max },{\Vert A\Vert }_{\max ,\infty }:=\underset{\tau \in [0,t]}{\max }{\Vert A(\tau )\Vert }_{\max }.$$

$${\Vert \mathcal{E}\Vert }_{\diamond }:=\max \{{\Vert (\mathcal{E}\otimes \text{mathds}{1}_{\mathcal{H}})(B)\Vert }_1:{\Vert B\Vert }_1\le 1\},$$

$${\Vert \mathcal{U}-\mathcal{V}\Vert }_{\diamond }\le 2{\Vert U-V\Vert }_{\infty }.$$

$${\mathcal{O}}_{\text{loc}}|j,s⟩$$

$${\mathcal{O}}_{\text{val}}|\tau ,j,k,z⟩$$

$$H(\tau )=\sum _{l=1}^L{\alpha }_l(\tau )H_l,$$

$$H(\tau )=\sum _{l=1}^L{H}_l(\tau ),$$

$$\text{textsc}Select(H):=\sum _{l=1}^L|l⟩⟨l|\otimes H_l,$$

$$U_m{Q}^{{\tau }_m}{U}_{m-1}\cdots U_1{Q}^{{\tau }_1}{U}_0,$$

$$\left\lVert e^{-itG}-\bigg{(}e^{-i\frac{t}{r}\beta_{1}G_{1}}\cdots e^{-i\frac{t}{r}\beta_{L}G_{L}}\bigg{)}^{r}\right\rVert_{\infty}=O\bigg{(}\frac{(\left\lVert\beta\right\rVert_{1}t)^{2}}{r}\bigg{)}.$$

$$\text{textsc}Select(\text{textsc}Exp-G)=\sum _{l=1}^L|l⟩⟨l|\otimes e^{-i\pi G_l}$$

$$\left\lVert e^{-itG}-\bigg{(}e^{-i\frac{t}{r}\beta_{1}G_{1}}\cdots e^{-i\frac{t}{r}\beta_{L}G_{L}}\bigg{)}^{r}\right\rVert_{\infty}=O(\epsilon).$$

$$O\bigg{(}T\frac{\log(T/\epsilon)}{\log\log(T/\epsilon)}\bigg{)}=O\bigg{(}\left\lVert\beta\right\rVert_{1}t\frac{\log(\left\lVert\beta\right\rVert_{1}t/\epsilon)}{\log\log(\left\lVert\beta\right\rVert_{1}t/\epsilon)}\bigg{)}$$

$$O\bigg{(}d^{2}\left\lVert H\right\rVert_{\max,1}\frac{\log(d^{2}\left\lVert H\right\rVert_{\max,\infty}t/\epsilon)}{\log\log(d^{2}\left\lVert H\right\rVert_{\max,\infty}t/\epsilon)}\bigg{)}$$

$$\displaystyle\left\lVert\exp_{\mathcal{T}}\biggl{(}-i\int_{\frac{kt}{r}}^{\frac{(k+1)t}{r}}\mathrm{d}\tau\,H(\tau)\biggr{)}-e^{-i\frac{t}{r}H\big{(}\frac{kt}{r}\big{)}}\right\rVert_{\infty}$$

$$\displaystyle\leq\int_{\frac{kt}{r}}^{\frac{(k+1)t}{r}}\mathrm{d}s\bigg{(}s-\frac{kt}{r}\bigg{)}\left\lVert H^{\prime}\right\rVert_{\infty,\infty}$$

Full text

1. 1 Introduction

2. 2 Preliminaries

3. 2.1 Time-dependent Hamiltonian evolution

4. 2.2 Notation for norms

5. 2.3 Hamiltonian input models

6. 2.4 Simulation algorithms with $L^{1}$-norm scaling

7. 2.5 Query complexity with $L^{1}$-norm scaling

8. 3 Continuous qDRIFT

9. 3.1 A classical sampler of time-dependent Hamiltonians

10. 3.2 Universality

11. 3.3 Complexity of the continuous qDRIFT algorithm

12. 4 Rescaled Dyson-series algorithm

13. 4.1 A rescaling principle for the Schrödinger equation

14. 4.2 Complexity of the rescaled Dyson-series algorithm

15. 5 Applications to chemistry and scattering theory

16. 6 Discussion

17. A Continuous qDRIFT and Hamiltonian averaging

18. A.1 Hamiltonian averaging

19. A.2 Implementing Hamiltonian averaging by continuous qDRIFT

1 Introduction

Simulating the Hamiltonian dynamics of a quantum system is one of the most promising applications of a quantum computer. The apparent classical intractability of simulating quantum dynamics led Feynman [25] and others to propose the idea of quantum computation. Quantum computers can simulate various physical systems, including condensed matter physics [3], quantum field theory [28], and quantum chemistry [2, 14, 48, 37]. The study of quantum simulation has also led to the discovery of new quantum algorithms, such as algorithms for linear systems [27], differential equations [10], semidefinite optimization [11], formula evaluation [23], quantum walk [17], and ground-state and thermal-state preparation [41, 20].

Let $H(\tau)$ be a Hamiltonian defined for $0\leq\tau\leq t$. The problem of Hamiltonian simulation is to approximate the evolution \exp_{\mathcal{T}}\big{(}{-}i\int_{0}^{t}\mathrm{d}\tau\ H(\tau)\big{)} using a quantum circuit comprised of elementary quantum gates, where $\exp_{\mathcal{T}}$ denotes the time-ordered matrix exponential. If the Hamiltonian $H(\tau)=H$ does not depend on time, the evolution operator can be represented in closed form as $e^{-itH}$. Then the problem can be greatly simplified and it has been thoroughly studied by previous works on quantum simulation [31, 1, 5, 9, 4, 8, 35, 33, 32, 7, 13, 29, 36, 18, 19, 16].

Simulating a general time-dependent Hamiltonian $H(\tau)$ naturally subsumes the time-independent case, and can be applied to devising quantum control schemes [40, 38], describing quantum chemical reactions [12], and implementing adiabatic quantum algorithms [22]. However, the problem becomes considerably harder and there are fewer quantum algorithms available. Wiebe, Berry, Høyer, and Sanders designed a time-dependent Hamiltonian simulation algorithm based on higher-order product formulas [50]. They assume that $H(\tau)$ is smooth up to a certain order and they give an example in which a desired approximation cannot be achieved due to the non-differentiability of the Hamiltonian. The smoothness assumption is relaxed in subsequent work by Poulin, Qarry, Somma, and Verstraete [42] based on techniques of Hamiltonian averaging and Monte Carlo estimation. The fractional-query algorithm of Berry, Childs, Cleve, Kothari, and Somma can also simulate time-dependent Hamiltonians [6], with an exponentially improved dependence on precision and only logarithmic dependence on the derivative of the Hamiltonian. A related quantum algorithm for time-dependent Hamiltonian simulation was suggested by Berry, Childs, Cleve, Kothari, and Somma based on the truncated Dyson series [8], which is analyzed explicitly in [29, 36].

In this paper, we study time-dependent Hamiltonian simulation based on a simple intuition: the difficulty of simulating a quantum system should depend on the integrated norm of the Hamiltonian. To elaborate, first consider the special case of simulating a time-independent Hamiltonian. The complexity of such a simulation depends on $t\left\lVert H\right\rVert$ [15], where $\left\lVert\cdot\right\rVert$ is a matrix norm that quantifies the size of the Hamiltonian. It is common to express the complexity in terms of the spectral norm $\left\lVert H\right\rVert_{\infty}$ (i.e., the Schatten $\infty$-norm), which quantifies the maximum energy of $H$.

In the general case where the Hamiltonian $H(\tau)$ is time dependent, we expect a quantum simulation algorithm to depend on the Hamiltonian locally in time, and therefore to have complexity that scales with the integrated norm $\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\rVert$. This is the $L^{1}$ norm of $\left\lVert H(\tau)\right\rVert$ when viewed as a function of $\tau$, so we say such an algorithm has $L^{1}$-norm scaling. Surprisingly, existing analyses of quantum simulation algorithms fail to achieve this complexity; rather, their gate complexity scales with the worst-case cost $t\max_{\tau\in[0,t]}\left\lVert H(\tau)\right\rVert$. It is therefore reasonable to question whether our intuition is correct, or if there exist faster time-dependent Hamiltonian simulation algorithms that can exploit this intuition.111For the Dyson-series approach, Low and Wiebe claimed that the worst-case scaling may be avoided by a proper segmentation of the time interval [36]. However, it is unclear how their analysis can be formalized to give an algorithm with complexity that scales with the $L^{1}$ norm. In Section 4, we propose a rescaling principle for the Schrödinger equation and develop a rescaled Dyson-series algorithm with $L^{1}$-norm scaling.

Our work answers this question by providing multiple faster quantum algorithms for time-dependent Hamiltonian simulation. These algorithms have gate complexity that scales with $\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\rVert$, in contrast to the best previous scaling of $t\max_{\tau\in[0,t]}\left\lVert H(\tau)\right\rVert$. As the norm inequality $\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\rVert\leq t\max_{\tau\in[0,t]}\left\lVert H(\tau)\right\rVert$ always holds but is not saturated in general, these algorithms provide a strict speedup over existing algorithms. We further analyze an application to simulating scattering processes in quantum chemistry, showing that our improvement can be favorable in practice.

We introduce notation and terminology and state our assumptions in Section 2. Following standard assumptions about quantum simulation, we consider two different input models of Hamiltonians. The first is the sparse matrix (SM) model common for analyzing Hamiltonian simulation in general, in which the Hamiltonian is assumed to be sparse and access to the locations and values of nonzero matrix elements are provided by oracles. We quantify the complexity of a simulation algorithm by the number of queries and additional gates it uses. The second model, favorable for practical applications such as condensed matter physics and quantum chemistry simulation, assumes that the Hamiltonian can be explicitly decomposed as a linear combination of unitaries (LCU), where the coefficients are efficiently computable on a classical computer and the summands can be exponentiated and controlled on a quantum computer. We ignore the cost of implementing the coefficient oracle and focus mainly on the gate complexity. Quantum simulation algorithms can sometimes work more efficiently in other input models, but we study these two models since they are versatile and provide a fair comparison of the gate complexity.

Reference [6] claims that the fractional-query algorithm can simulate time-dependent Hamiltonians with $L^{\infty}$-norm scaling. However, it is not hard to see that its query complexity in fact scales with the $L^{1}$ norm. While we do not show how to achieve this scaling for the gate complexity, our analysis is simple and suggests that such a result might be possible. We analyze the query complexity of the fractional-query algorithm in Section 2.5.

We develop two new techniques to simulate time-dependent Hamiltonians with $L^{1}$-norm scaling. Our first technique is a classical sampling protocol for time-dependent Hamiltonians. In this protocol, we randomly sample a time $\tau\in[0,t]$ and evolve under the time-independent Hamiltonian $H(\tau)$, where the probability distribution is designed to favor those $\tau$ with large $\left\lVert H(\tau)\right\rVert$. Campbell introduced a discrete sampling scheme for time-independent Hamiltonian simulation [13] and our protocol can be viewed as its continuous analog, which we call “continuous qDRIFT”. We show that continuous qDRIFT is universal, in the sense that any Hamiltonian simulable by [13] can be simulated by continuous qDRIFT with the same complexity. In addition, we shave off a multiplicative factor in the analysis of [13] by explicitly evaluating the rate of change of the evolution with respect to scaling the Hamiltonian. Continuous qDRIFT and its analysis are detailed in Section 3. Our algorithm is also similar in spirit to the approach of Poulin et al. [43] based on Hamiltonian averaging and Monte Carlo estimation, although their algorithm does not have $L^{1}$-norm scaling. We discuss the relationship between these two approaches in Appendix A.

We also present a general principle for rescaling the Schrödinger equation in Section 4. In the rescaled Schrödinger equation, the time-dependent Hamiltonian $H(\tau)$ has the same norm for all $\tau\in[0,t]$, so the norm inequality $\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\rVert\leq t\max_{\tau\in[0,t]}\left\lVert H(\tau)\right\rVert$ holds with equality. Using this principle, we show that the simulation algorithm based on the truncated Dyson series [8, 36, 29] can also be improved to have $L^{1}$-norm scaling.

To illustrate how our results might be applied, we identify a specific problem in quantum chemistry for which our $L^{1}$-norm improvement is advantageous: semi-classical scattering of molecules in a chemical reaction. For such a simulation, $\left\lVert H(\tau)\right\rVert$ changes dramatically throughout the evolution, so its $L^{1}$ norm can be significantly smaller than its $L^{\infty}$ norm. We discuss this application further in Section 5.

Finally, we conclude in Section 6 with a brief discussion of the results and some open questions.

2 Preliminaries

2.1 Time-dependent Hamiltonian evolution

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$. By default, we assume that $H(\tau)$ is continuously differentiable and $H(\tau)\neq 0$ everywhere, and we defer the discussion of pathological cases to Section 6. If the Hamiltonian $H(\tau)=H$ is time independent, the evolution is given in closed form by the matrix exponential $e^{-itH}$. However, there exists no such closed-form expression for a general $H(\tau)$ and we instead represent the evolution by \exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\big{)}, where $\exp_{\mathcal{T}}$ denotes the time-ordered matrix exponential. We have

[TABLE]

If $G(\tau)$ is another time-dependent Hamiltonian, the evolutions generated by $H(\tau)$ and $G(\tau)$ have distance bounded by the following lemma.

Lemma 1** ($L^{1}$-norm distance bound of time-ordered evolutions [45, Appendix B]).**

Let $H(\tau)$ and $G(\tau)$ be time-dependent Hamiltonians defined on the interval $0\leq\tau\leq t$. Then,

[TABLE]

Here, $\left\lVert\cdot\right\rVert_{\infty}$ denotes the spectral norm.

We will abbreviate the evolution operator as \mathrm{E}(t,s):=\exp_{\mathcal{T}}\big{(}{-}i\int_{s}^{t}\mathrm{d}\tau\,H(\tau)\big{)} when there is no ambiguity. In the special case where $H(\tau)=H$ is time independent, the evolution $e^{-itH}$ only depends on the time duration so we denote $\mathrm{E}(t):=\mathrm{E}(t,0)$. Therefore, we have the differential equation

[TABLE]

We may further obtain an integral representation of $\mathrm{E}(t,0)$. To this end, we apply the fundamental theorem of calculus to the Schrödinger equation and obtain

[TABLE]

Equivalently, $\mathrm{E}(t,0)$ satisfies the integral equation

[TABLE]

For any $0\leq s\leq t$, the evolution operator satisfies the multiplicative property

[TABLE]

The operator $\mathrm{E}(0,t)$ with $t\geq 0$ is understood as the inverse evolution operator

[TABLE]

For a thorough mathematical treatment of time-dependent Hamiltonian evolution, we refer the reader to [21]. Finally, the quantum channel corresponding to the unitary evolution $\mathrm{E}(t,0)$ is denoted as $\mathcal{E}(t,0)$ and is defined by

[TABLE]

For time-independent Hamiltonians, we denote $\mathcal{E}(t):=\mathcal{E}(t,0)$.

2.2 Notation for norms

We introduce norm notation for vectors, matrices, operator-valued functions, and linear maps on the space of matrices.

Let $\alpha=\begin{bmatrix}\alpha_{1}&\alpha_{2}&\cdots&\alpha_{L}\end{bmatrix}\in\mathbb{C}^{L}$ be an $L$-dimensional vector. We use $\left\lVert\alpha\right\rVert_{p}$ to represent the vector $\ell_{p}$ norm of $\alpha$. Thus,

[TABLE]

For a matrix $A$, we define $\left\lVert A\right\rVert_{p}$ to be the Schatten $p$-norm of $A$ [47, 51]. We have

[TABLE]

Finally, if $f\colon[0,t]\rightarrow\mathbb{C}$ is a continuous function, we use $\left\lVert f\right\rVert_{p}$ to mean the $L^{p}$ norm of the function $f$. Thus,

[TABLE]

We combine these norms to obtain norms for vector-valued and operator-valued functions. Let $\alpha\colon[0,t]\rightarrow\mathbb{C}^{L}$ be a continuous vector-valued function, with the $j$th coordinate at time $\tau$ denoted $\alpha_{j}(\tau)$. We use $\left\lVert\alpha\right\rVert_{p,q}$ to mean that we take the $\ell_{p}$ norm $\left\lVert\alpha(\tau)\right\rVert_{p}$ for every $\tau$ and compute the $L^{q}$ norm of the resulting scalar function. For example,

[TABLE]

Note that $\left\lVert\alpha(\tau)\right\rVert_{p}$ is continuous as a function of $\tau$, so $\left\lVert\alpha\right\rVert_{p,q}$ is well defined and is indeed a norm for vector-valued functions. Similarly, we also define $\left\lVert A\right\rVert_{p,q}$ for a continuous operator-valued function by taking the Schatten $p$-norm $\left\lVert A(\tau)\right\rVert_{p}$ for every $\tau$ and computing the $L^{q}$ norm of the resulting scalar function. In rare cases, we will also encounter time-dependent linear combinations of operators of the form $A(\tau)=\sum_{l=1}^{L}A_{l}(\tau)$, and we write $\left\lVert A\right\rVert_{p,q,r}$ to mean that we take the Schatten $p$-norm $\left\lVert A_{l}(\tau)\right\rVert_{p}$ of each summand, and apply the $\ell_{q}$ norm and $L^{r}$ norm to the resulting vector-valued functions. For example,

[TABLE]

We also define $\left\lVert A\right\rVert_{\max}$ as the largest matrix element of $A$ in absolute value,

[TABLE]

The norm $\left\lVert A\right\rVert_{\max}$ is a vector norm of $A$ but does not satisfy the submultiplicative property of a matrix norm. It relates to the spectral norm by the inequality [15, Lemma 1]

[TABLE]

If $A$ is a continuous operator-valued function, we interpret $\left\lVert A\right\rVert_{\max,q}$ in a similar way as above. Therefore,

[TABLE]

Finally, we define a norm for linear maps on the space of matrices. Let $\mathcal{E}\colon A\mapsto\mathcal{E}(A)$ be a linear map on the space of matrices on $\mathcal{H}$. The diamond norm of $\mathcal{E}$ is

[TABLE]

where the maximization is taken over all matrices $B$ on $\mathcal{H}\otimes\mathcal{H}$ satisfying $\left\lVert B\right\rVert_{1}\leq 1$. Below is a useful bound on the diamond-norm distance between two unitary channels.

Lemma 2** (Diamond-norm distance between unitary channels [9, Lemma 7]).**

Let $V$ and $U$ be unitary matrices, with associated quantum channels $\mathcal{V}:\rho\mapsto V\rho V^{\dagger}$ and $\mathcal{U}:\rho\mapsto U\rho U^{\dagger}$. Then,

[TABLE]

The sampling-based algorithm (Section 3) produces a channel close to \mathcal{E}(t,0)(\rho)=\exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau H(\tau)\big{)}\rho\exp_{\mathcal{T}}^{\dagger}\big{(}-i\int_{0}^{t}\mathrm{d}\tau H(\tau)\big{)}, and its error is naturally quantified by the diamond-norm distance. Other simulation algorithms such as the Dyson-series approach (Section 4) produce operators that are close to the unitary \exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau H(\tau)\big{)}, and we quantify their error in terms of the spectral norm. For a fair comparison one may instead describe all simulation algorithms using quantum channels and use the diamond-norm distance as the unified error metric. By Lemma 2, we lose at most a factor of $2$ in this conversion.

2.3 Hamiltonian input models

Quantum simulation algorithms may have different performance depending on the choice of the input model of Hamiltonians. In this section, we describe two input models that are extensively used in previous works: the sparse matrix (SM) model and the linear-combination-of-unitaries (LCU) model. We also discuss other input models that will be used in later sections.

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$. In the SM model, we assume that $H(\tau)$ is $d$-sparse in the sense that the number of nonzero matrix elements within each row and column throughout the entire interval $[0,t]$ is at most $d$. We assume that the locations of the nonzero matrix elements are time independent. Access to the Hamiltonian is given through the oracles

[TABLE]

Here, $\mathrm{col}(j,s)$ returns the column index of the $s$th element in the $j$th row that may be nonzero over the entire time interval $[0,t]$. We quantify the complexity of a quantum simulation algorithm by the number of queries it makes to $\mathcal{O}_{\text{loc}}$ and $\mathcal{O}_{\text{val}}$, together with the number of additional elementary gates it requires. Such a model includes many realistic physical systems and is well-motivated from a theoretical perspective [27].

As the following lemma shows, a $d$-sparse Hamiltonian can be efficiently decomposed as a sum of $1$-sparse terms.

Lemma 3** (Decomposition of sparse Hamiltonians [6, Lemma 4.3 and 4.4]).**

Let $H$ be a time-independent $d$-sparse Hamiltonian accessed through the oracles $\mathcal{O}_{\text{loc}}$ and $\mathcal{O}_{\text{val}}$. Then

there exists a decomposition $H=\sum_{j=1}^{d^{2}}H_{j}$, where each $H_{j}$ is $1$-sparse with $\left\lVert H_{j}\right\rVert_{\max}\leq\left\lVert H\right\rVert_{\max}$, and a query to any $H_{j}$ can be simulated with $O(1)$ queries to $H$; and 2. 2.

for any $\gamma>0$, there exists an approximate decomposition222Reference [6] uses [6, Lemma 4.3] and the triangle inequality to show that \bigl{\|}{H-\gamma\sum_{j=1}^{\eta}G_{j}}\bigr{\|}_{\max}\leq\sqrt{2}\gamma d^{2}. However, this bound can be tightened to $\sqrt{2}\gamma$, since the max-norm distance depends on the largest error from rounding off the $d^{2}$ $1$-sparse matrices. * \bigl{\|}{H-\gamma\sum_{j=1}^{\eta}G_{j}}\bigr{\|}_{\max}\leq\sqrt{2}\gamma, where \eta=O\big{(}d^{2}\left\lVert H\right\rVert_{\max}/\gamma\big{)}, each $G_{j}$ is $1$-sparse with eigenvalues $\pm 1$, and a query to any $G_{j}$ can be simulated with $O(1)$ queries to $H$.*

For the LCU model, we suppose that the Hamiltonian $H(\tau)$ admits a decomposition

[TABLE]

where the coefficients $\alpha_{l}(\tau)\geq 0$ are continuously differentiable and nonzero everywhere, and the matrices $H_{l}$ are both unitary and Hermitian. We assume that the coefficients $\alpha_{l}(\tau)$ can be efficiently computed by a classical oracle $\mathcal{O}_{\text{coeff}}$, and we ignore the classical cost of implementing such an oracle. We further assume that each $|0\rangle\langle 0|\otimes I+|1\rangle\langle 1|\otimes H_{l}$ can be implemented with gate complexity $g_{c}$, and each $|0\rangle\langle 0|\otimes I+|1\rangle\langle 1|\otimes e^{-i\tau H_{l}}$ for an arbitrarily large $\tau$ can be performed with $g_{e}$ gates. Such a setting is common in the simulation of condensed matter physics and quantum chemistry. We quantify the complexity of a simulation algorithm by the number of elementary gates it uses.

Quantum simulation algorithms can sometimes work in other input models. For example, the qDRIFT protocol introduced in Section 3 requires only that the Hamiltonians have the form

[TABLE]

where the Hermitian-valued functions $H_{l}(\tau)$ are continuous, nonzero everywhere, and can be efficiently exponentiated on a quantum computer. We call this the linear combination (LC) model. On the other hand, the Dyson-series algorithm can be described in terms of the Select operation

[TABLE]

irrespective of how this operation is implemented. We consider the SM and LCU models for all the time-dependent simulation algorithms so that we can give a fair comparison of their complexity.

2.4 Simulation algorithms with $L^{1}$-norm scaling

We now explain the meaning of $L^{1}$-norm scaling in the SM and the LCU models. Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$. We say that an algorithm in the SM model simulates $H(\tau)$ with $L^{1}$-norm scaling if, given any continuously differentiable upper bound $\Lambda_{\max}(\tau)\geq\left\lVert H(\tau)\right\rVert_{\max}$, the algorithm has query complexity and gate complexity that scale with $\left\lVert\Lambda_{\max}\right\rVert_{1}=\int_{0}^{t}\mathrm{d}\tau\,\Lambda_{\max}(\tau)$ up to logarithmic factors. The norm bound $\Lambda_{\max}(\tau)$, together with other auxiliary information, must be accessed by the quantum simulation algorithm; we assume such quantities can be computed efficiently.

In the LCU model, we are given a time-dependent Hamiltonian with the decomposition $H(\tau)=\sum_{l=1}^{L}\alpha_{l}(\tau)H_{l}$. We say that an algorithm has $L^{1}$-norm scaling if, for any continuously differentiable vector-valued function $\Lambda(t)$ with $\Lambda_{l}(\tau)\geq\alpha_{l}(\tau)$, the algorithm has query and gate complexity that scale with $\left\lVert\Lambda\right\rVert_{\infty,1}=\int_{0}^{t}\mathrm{d}\tau\max_{l}\Lambda_{l}(\tau)$ up to logarithmic factors.

For better readability, we express the complexity of simulation algorithms in terms of the norm of the original Hamiltonian, such as $\left\lVert H\right\rVert_{\max,1}$ and $\left\lVert\alpha\right\rVert_{\infty,1}$, instead of the upper bounds $\left\lVert\Lambda_{\max}\right\rVert_{1}$ and $\left\lVert\Lambda\right\rVert_{\infty,1}$. We use standard asymptotic notation, with $O$, $\Omega$, and $\Theta$ representing asymptotic upper, lower, and tight bounds, respectively. We also suppress logarithmic factors using the $\widetilde{O}$ notation when the complexity expression becomes too complicated. Table 1 compares the results of this paper with previous results on simulating time-dependent Hamiltonians.

Our goal is to develop simulation algorithms that scale with the $L^{1}$-norm with respect to the time variable $\tau$, for both query complexity and gate complexity. We start by reexamining the fractional-query approach. It was mentioned in [6] that this approach can simulate time-dependent Hamiltonians with $L^{\infty}$-norm scaling, but we find that its query complexity scales with the $L^{1}$ norm. We give this improved analysis in the next section.

2.5 Query complexity with $L^{1}$-norm scaling

We begin by reviewing the result of [6] for simulating time-independent Hamiltonians. We assume that the Hamiltonian is given by a linear combination of unitaries $G=\sum_{l=1}^{L}\beta_{l}G_{l}$ with nonnegative coefficients $\beta_{l}$. Here, $G_{l}$ are both unitary and Hermitian, so they are reflections and their eigenvalues are $\pm 1$.

We say that a quantum operation is a fractional-query algorithm if it is of the form

[TABLE]

where $Q$ is unitary with eigenvalues $\pm 1$, $U_{j}$ are unitary operations, and $\tau_{j}\geq 0$. Here, we regard $Q$ as the oracle and $U_{j}$ as non-query operations, so this algorithm has fractional-query complexity $\sum_{j=1}^{m}\tau_{j}$. A quantum algorithm that makes (discrete) queries to $Q$ is a fractional-query algorithm with $\tau_{j}=1$. Conversely, any fractional-query algorithm can be efficiently simulated in the discrete query model. In particular, an algorithm with fractional-query complexity $T$ can be simulated with error at most $\epsilon$ using O\big{(}T\frac{\log(T/\epsilon)}{\log\log(T/\epsilon)}\big{)} discrete queries [6, Lemma 3.8].

To apply the fractional-query approach, we approximate the evolution under $G$ using the first-order product formula

[TABLE]

Observe that $e^{-i\pi G_{l}}$ are unitary operations with eigenvalues $\pm 1$, so \big{(}e^{-i\frac{t}{r}\beta_{1}G_{1}}\cdots e^{-i\frac{t}{r}\beta_{L}G_{L}}\big{)}^{r} can be viewed as a fractional-query algorithm with query complexity $O(\left\lVert\beta\right\rVert_{1}t)$, provided that we can make fractional queries to multiple oracles $e^{-i\pi G_{1}},\ldots,e^{-i\pi G_{L}}$. This can be realized by a standard fractional-query algorithm accessing the single oracle

[TABLE]

with the same query complexity [6, Theorem 4.1].

To simulate with accuracy $\epsilon$, we set r=O\big{(}(\left\lVert\beta\right\rVert_{1}t)^{2}/\epsilon\big{)} to ensure that

[TABLE]

We now convert this multi-oracle algorithm to a single-oracle algorithm with the same fractional-query complexity $T=O(\left\lVert\beta\right\rVert_{1}t)$ and, with precision $O(\epsilon)$, implement it in the discrete query model. Altogether, this approach makes

[TABLE]

queries to the operation $\textsc{Select}(\textsc{Exp-}G)=\sum_{l=1}^{L}|l\rangle\langle l|\otimes e^{-i\pi G_{l}}$.

As mentioned in [6], the fractional-query approach can also be used to simulate time-dependent Hamiltonians by replacing (24) with a product-formula decomposition of the time-ordered evolution. However, [6] only gives a brief discussion of this issue and the claimed complexity has only $L^{\infty}$ scaling. We now give an improved analysis of this algorithm for the SM model, showing that its query complexity achieves $L^{1}$-norm scaling.

Theorem 4** (Fractional-query algorithm with $L^{1}$-norm scaling (SM)).**

A $d$-sparse time-dependent Hamiltonian $H(\tau)$ acting on $n$ qubits can be simulated for time $\tau\in[0,t]$ with accuracy $\epsilon$ using

[TABLE]

queries to the oracles $\mathcal{O}_{\text{loc}}$, $\mathcal{O}_{\text{val}}$.

Proof.

For readability, we assume that $\left\lVert H\right\rVert_{\max,1}$, $\left\lVert H\right\rVert_{\max,\infty}$, and $\left\lVert H^{\prime}\right\rVert_{\infty,\infty}$ are the norm upper bounds provided to the algorithm. We first decompose \exp_{\mathcal{T}}\bigl{(}{-}i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\bigr{)} into a product of evolutions of time-independent Hamiltonians $H(kt/r)$, each evolving for time $t/r$. By Lemma 1, we have

[TABLE]

which implies

[TABLE]

To approximate with precision $\epsilon$, it suffices to choose

[TABLE]

Note that we use $O$ here because we can choose $r$ to be the minimum integer satisfying (30), giving an upper bound on the number of steps that suffice to achieve error at most $\epsilon$.

We then decompose the evolution under each time-independent sparse Hamiltonian $H(kt/r)$ for time $t/r$ with precision $O(\epsilon/r)$. By Lemma 3, $H(kt/r)$ can be decomposed into a sum of $\eta=O(d^{2}\left\lVert H(kt/r)\right\rVert_{\max}/\gamma)$ terms $G_{j}(kt/r)$ such that

[TABLE]

Furthermore, each $G_{j}$ is $1$-sparse and Hermitian with eigenvalues $\pm 1$ and the value and location of each non-zero matrix element in $G_{j}$ can be accessed using $O(1)$ queries to $H$. We choose \gamma=O\big{(}\epsilon/td\big{)} so that

[TABLE]

In the third line we have used the inequality between the spectral norm and max norm, in the fourth line we have used the bound on the max norm (32), and in the fifth line we have used \gamma=O\big{(}\epsilon/td\big{)}. This implies \eta=O\big{(}d^{3}\left\lVert H(kt/r)\right\rVert_{\max}t/\epsilon\big{)} and the fractional query complexity is

[TABLE]

We apply the first-order product formula to obtain

[TABLE]

Therefore, it is possible to choose $r$ as

[TABLE]

such that the error of the first-order product-formula decomposition is at most

[TABLE]

By choosing $r$ as the maximum of (31) and (35), we ensure that the error in each of the $r$ time steps is $O(\epsilon/r)$, so the total error is $O(\epsilon)$.

Altogether, we find a fractional-query algorithm with total query complexity

[TABLE]

and error

[TABLE]

We now convert this multi-oracle algorithm to a single-oracle algorithm with the same fractional-query complexity and, with precision $O(\epsilon)$, implement it in the discrete query model. The single oracle for the standard fractional-query algorithm is now

[TABLE]

This oracle encodes the time-dependence of $H$ in an ancilla. The operators $U_{j}$ in the fractional-query algorithm then need to increment the time register.

Altogether, we make

[TABLE]

discrete queries.

We now show how the query complexity of this approach achieves $L^{1}$-norm scaling. The intuition is that the total query complexity \widetilde{O}\big{(}\sum_{k=0}^{r-1}d^{2}\left\lVert H(kt/r)\right\rVert_{\max}t/r\big{)} should be close to \widetilde{O}\big{(}d^{2}\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\rVert_{\max}\big{)} when $r$ is sufficiently large. Specifically,

[TABLE]

To achieve an additive error of $\delta$, it suffices to choose

[TABLE]

Since $\delta$ can be made arbitrarily close to [math], we have the total query complexity of

[TABLE]

as claimed. ∎

The above analysis shows that the fractional-query algorithm can simulate a time-dependent Hamiltonian with query complexity that scales with the $L^{1}$-norm. However, this approach does not directly give a useful result for the gate complexity. The difficulty arises from the factor of $g$ in [6, Proof of Theorem 1.1], which corresponds to the complexity of applying a sequence of driving operations $U_{j}$. These operations need to increment $k$ (indexing the time), as well as $l$, which takes $\eta(k)$ values depending on $k$. Applying the sequence of operations $U_{j}$ therefore requires determining new values of $l$ and $k$, which can depend on the sum of $\eta(k)$ over a long sequence of values of $k$. This will introduce significant gate complexity, so a fast algorithm would require a more efficient procedure for implementing the driving operations.

Instead, we develop other quantum algorithms that achieve $L^{1}$-norm scaling for not only the query complexity but also the gate complexity. We employ two main techniques: the continuous qDRIFT sampling protocol (Section 3) and a rescaling principle for the Schrödinger equation (Section 4).

3 Continuous qDRIFT

We show in Section 3.2 that continuous qDRIFT is universal, in the sense that any time-independent Hamiltonian simulable by the algorithm of [13] can be simulated by our protocol. We then discuss the simulation complexity in both the SM and the LCU models in Section 3.3.

The continuous qDRIFT protocol also has similarities with the approach of Poulin et al.  [43] based on Hamiltonian averaging and Monte Carlo sampling, although their approach does not have $L^{1}$-norm scaling. We give a detailed comparison between these two approaches in Appendix A.

3.1 A classical sampler of time-dependent Hamiltonians

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$. For this section only, we relax our requirements on the Hamiltonians: we assume that $H(\tau)$ is nonzero everywhere and is continuous except on a finite number of points. We further suppose that each $H(\tau)$ can be directly exponentiated on a quantum computer. The ideal evolution under $H(\tau)$ for time $t$ is given by \mathrm{E}(t,0)=\exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\big{)} and the corresponding quantum channel is

[TABLE]

The high-level idea of the sampling algorithm is to approximate the ideal channel by a mixed unitary channel

[TABLE]

where $p(\tau)$ is a probability density function defined for $0\leq\tau\leq t$. This channel can be realized by a classical sampling protocol. With a proper choice of $p(\tau)$, this channel approximates the ideal channel and can thus be used for quantum simulation.

We begin with a full definition of $\mathcal{U}(t,0)$. Inspired by [13], we choose $p(\tau)$ to be biased toward those $\tau$ with large $\left\lVert H(\tau)\right\rVert_{\infty}$. A natural choice is

[TABLE]

Note that $\mathcal{U}(t,0)$ is a valid quantum channel (in particular, $p(\tau)$ can never be zero). Furthermore, it can be implemented with unit cost: for any input state $\rho$, we randomly sample a value $\tau$ according to $p(\tau)$ and perform $e^{-i{H(\tau)}/{p(\tau)}}$. Note also that ${H(\tau)}/{p(\tau)}$ in the exponential implicitly depends on $t$. Indeed, $\left\lVert H\right\rVert_{\infty,1}$ includes an integral over time, so $p(\tau)$ decreases with the total evolution time $t$. We call this classical sampling protocol and the channel it implements “continuous qDRIFT”.

This protocol assumes that the spectral norm $\left\lVert H(\tau)\right\rVert_{\infty}$ is known a priori and that we can efficiently sample from the distribution $p(\tau)$. In practice, it is often easier to obtain a spectral-norm upper bound $\Lambda(\tau)\geq\left\lVert H(\tau)\right\rVert_{\infty}$. Such an upper bound can also be used to implement continuous qDRIFT, provided that it has only finitely many discontinuities. Specifically, we define

[TABLE]

so $p_{\Lambda}(\tau)$ is a probability density function. Using $p_{\Lambda}$ to implement continuous qDRIFT, we obtain the channel

[TABLE]

whose analysis is similar to that presented here. For readability, we assume that we can efficiently sample from $p(\tau)=\left\lVert H(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1}$ and we analyze $\mathcal{U}(t,0)$.

We show that continuous qDRIFT approximates the ideal channel with error that depends on the $L^{1}$-norm.

Theorem 5** ($L^{1}$-norm error bound for continuous qDRIFT, short-time version).**

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$; assume it is continuous except on a finite number of points and nonzero everywhere. Define \mathrm{E}(t,0)=\exp_{\mathcal{T}}\bigl{(}-i\int_{0}^{t}\mathrm{d}\tau H(\tau)\bigr{)} and let $\mathcal{E}(t,0)(\cdot)=\mathrm{E}(t,0)(\cdot)\mathrm{E}^{\dagger}(t,0)$ be the corresponding quantum channel. Let $\mathcal{U}(t,0)$ be the continuous qDRIFT channel

[TABLE]

where $p(\tau)=\left\lVert H(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1}$. Then

[TABLE]

Note that this bound is only useful when $t$ is small enough that the right-hand side of (50) is less than $1$ (the norm $\left\lVert H\right\rVert_{\infty,1}$ involves an integral over $t$, so it increases with $t$).

To prove this theorem, we need a formula that computes the rate at which the evolution operator changes when the Hamiltonian is scaled. To illustrate the idea, consider the degenerate case where the Hamiltonian $H$ is time independent. Then the evolution under $H$ for time $t$ is given by $e^{-itH}$. A direct calculation shows that

[TABLE]

so the rate is $-itHe^{-itsH}$ in the time-independent case. This calculation becomes significantly more complicated for a time-dependent Hamiltonian. The following lemma gives an explicit formula for

[TABLE]

We sketch the proof of this formula for completeness, but refer the reader to [21, p. 35] for mathematical justifications that are beyond the scope of this paper.

Lemma 6** (Hamiltonian scaling).**

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$ and assume it has finitely many discontinuities. Denote \mathrm{E}_{s}(t,v)=\exp_{\mathcal{T}}\bigl{(}-i\int_{v}^{t}\mathrm{d}\tau\,sH(\tau)\bigr{)}. Then,

[TABLE]

Proof sketch.

We first consider the special case where $H(\tau)$ is continuous in $\tau$. We invoke the variation-of-parameters formula [30, Theorem 4.9] to construct the claimed integral representation for $\frac{\mathrm{d}}{\mathrm{d}s}\mathrm{E}_{s}(t,v)$. To this end, we need to find a differential equation satisfied by $\frac{\mathrm{d}}{\mathrm{d}t}\frac{\mathrm{d}}{\mathrm{d}s}\mathrm{E}_{s}(t,v)$ and the corresponding initial condition \frac{\mathrm{d}}{\mathrm{d}s}\mathrm{E}_{s}(t,v)\big{\rvert}_{t=v}. We differentiate the Schrödinger equation $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{E}_{s}(t,v)=-isH(t)\mathrm{E}_{s}(t,v)$ with respect to $s$ to get

[TABLE]

Invoking the variation-of-parameters formula, we find an integral representation

[TABLE]

It thus remains to find the initial condition \frac{\mathrm{d}}{\mathrm{d}s}\mathrm{E}_{s}(t,v)\big{\rvert}_{t=v}.

We start from the Schrödinger equation $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{E}_{s}(t,v)=-isH(t)\mathrm{E}_{s}(t,v)$ and apply the fundamental theorem of calculus with initial condition $\mathrm{E}_{s}(v,v)=I$, obtaining the integral representation

[TABLE]

Differentiating this equation with respect to $s$ gives

[TABLE]

which implies

[TABLE]

Combining (55) and (58) establishes the claimed integral representation for $\frac{\mathrm{d}}{\mathrm{d}s}\mathrm{E}_{s}(t,v)$.

Now consider the case where $H(\tau)$ is piecewise continuous with one discontinuity at $t_{1}\in[v,t]$. We use the multiplicative property to break the evolution at $t_{1}$, so that each subevolution is generated by a continuous Hamiltonian. We have

[TABLE]

The general case of finitely many discontinuities follows by induction. ∎

Note that our argument implicitly assumes the existence of the derivatives and that we can interchange the order of $\frac{\mathrm{d}}{\mathrm{d}s}$ and $\frac{\mathrm{d}}{\mathrm{d}t}$. A rigorous justification of these assumptions is beyond the scope of the paper; we refer the reader to [21, p. 35] for details.

Proof of Theorem 5.

Define two parametrized quantum channels

[TABLE]

and observe that

[TABLE]

To bound the diamond-norm error $\left\lVert\mathcal{E}_{1}(t,0)-\mathcal{U}_{1}(t,0)\right\rVert_{\diamond}$, we should take a state $\sigma$ on the joint system of the original register and an ancilla register with the same dimension and upper bound $\left\lVert(\mathcal{E}_{1}(t,0)\otimes\mathds{1})(\sigma)-(\mathcal{U}_{1}(t,0)\otimes\mathds{1})(\sigma)\right\rVert_{1}$. For readability, we instead show how to bound the error $\left\lVert\mathcal{E}_{1}(t,0)(\rho)-\mathcal{U}_{1}(t,0)(\rho)\right\rVert_{1}$, but the derivation works in exactly the same way for the distance $\left\lVert(\mathcal{E}_{1}(t,0)\otimes\mathds{1})(\sigma)-(\mathcal{U}_{1}(t,0)\otimes\mathds{1})(\sigma)\right\rVert_{1}$ and the resulting bound is the same.

Invoking Lemma 6, we have

[TABLE]

Thus, the first derivatives of $\mathcal{E}_{s}(t,0)(\rho)$ and $\mathcal{U}_{s}(t,0)(\rho)$ at $s=0$ agree with each other:

[TABLE]

Applying the fundamental theorem of calculus twice, we obtain

[TABLE]

By properties of the Schatten norms and the definition $p(\tau)=\left\lVert H(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1}$, we find that

[TABLE]

Lemma 6 immediately yields an upper bound on $\left\lVert\frac{\mathrm{d}}{\mathrm{d}v}\mathrm{E}_{v}(t,0)\right\rVert_{\infty}$:

[TABLE]

It thus remains to bound $\left\lVert\frac{\mathrm{d}^{2}}{\mathrm{d}v^{2}}\mathrm{E}_{v}(t,0)\right\rVert_{\infty}$.

Using Lemma 6 twice, we have

[TABLE]

which implies

[TABLE]

We finally obtain the desired bound

[TABLE]

as claimed. ∎

The above error bound works well for a short-time evolution. When $t$ is large, in order to control the error of simulation, we divide the entire evolution into segments $[t_{j},t_{j+1}]$ with $0=t_{0}<t_{1}<\cdots<t_{r}=t$ and apply continuous qDRIFT within each. We employ a variable-time scheme to segment the evolution, so that our $L^{1}$-norm scaling result can be generalized to a long-time evolution. Specifically, we have:

Theorem 7** ($L^{1}$-norm error bound for continuous qDRIFT, long-time version).**

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$. Assume that it is continuous except at a finite number of points and nonzero everywhere. Define \mathrm{E}(t,0)=\exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\big{)} and let $\mathcal{E}(t,0)(\cdot)=\mathrm{E}(t,0)(\cdot)\mathrm{E}^{\dagger}(t,0)$ be the corresponding quantum channel. Let $\mathcal{U}(t,0)$ be the continuous qDRIFT channel

[TABLE]

where $p(\tau)=\left\lVert H(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1}$. Then, for any positive integer $r$, there exists a division $0=t_{0}<t_{1}<\cdots<t_{r}=t$, such that

[TABLE]

To ensure that the simulation error is at most $\epsilon$, it thus suffices to choose

[TABLE]

Proof.

The times $t_{1},\cdots,t_{r-1}$ are selected as follows. We aim to simulate with accuracy

[TABLE]

for each segment. To achieve this, we define $t_{1},\cdots,t_{r-1}$ so that

[TABLE]

The existence of such times is guaranteed by the intermediate value theorem. By telescoping, we find from Theorem 5 that

[TABLE]

which establishes the claimed error bound. ∎

3.2 Universality

We now show that the continuous qDRIFT method introduced above can be applied in the far more general LC model where the Hamiltonian is a sum of time-dependent terms. In this sense it can be regarded as a universal method.

Recall from Section 2.3 that in the general LC model, the Hamiltonian can be expressed as

[TABLE]

where each $H_{l}(\tau)$ is continuous and can be efficiently exponentiated on a quantum computer. This includes many familiar models as special cases:

1. (i)

Campbell considered simulating a time-independent Hamiltonian of the form $H=\sum_{l=1}^{L}\alpha_{l}H_{l}$, $\left\lVert H_{l}\right\rVert_{\infty}\leq 1$ [13], which is subsumed by the LC model with the time dependence dropped; 2. (ii)

if $H(\tau)$ is a time-dependent $d$-sparse Hamiltonian, then Lemma 3 shows that it can be decomposed in the form $H(\tau)=\sum_{j=1}^{d^{2}}H_{j}(\tau)$, which again belongs to the LC model as the exponentiation of $H_{j}(\tau)$ can be performed efficiently; and 3. (iii)

the LC model is naturally more general than LCU as each summand is not necessarily unitary.

It is not hard to design a classical sampler for time-dependent Hamiltonians in the LC model. A natural choice is

[TABLE]

where $p_{l}(\tau)$ is the probability distribution

[TABLE]

To analyze the performance of this sampler, we adapt the analysis in Theorem 5 and Theorem 7, which becomes more complicated as we are now sampling a discrete-continuous probability distribution $p_{l}(\tau)$. Fortunately, a significant amount of effort can be saved with the help of the following universal property.

Theorem 8** (Universality of continuous qDRIFT).**

Let $H(\tau)=\sum_{l=1}^{L}H_{l}(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$ that is nonzero everywhere. Assume that each $H_{l}(\tau)$ is continuous and nonzero everywhere. Define the probability distribution

[TABLE]

Then there exists a time-dependent Hamiltonian $G(\tau)$ defined for $0\leq\tau\leq t$ with finitely many discontinuities, such that the following correspondence holds:

$\left\lVert G\right\rVert_{\infty,1}=\left\lVert H\right\rVert_{\infty,1,1}$. 2. 2.

$\int_{0}^{t}\mathrm{d}\tau\,G(\tau)=\sum_{l=1}^{L}\int_{0}^{t}\mathrm{d}\tau\,H_{l}(\tau)$. 3. 3.

$\int_{0}^{t}\mathrm{d}\tau\,q(\tau)e^{-i\frac{G(\tau)}{q(\tau)}}\rho e^{i\frac{G(\tau)}{q(\tau)}}=\sum_{l=1}^{L}\int_{0}^{t}\mathrm{d}\tau\,p_{l}(\tau)e^{-i\frac{H_{l}(\tau)}{p_{l}(\tau)}}\rho e^{i\frac{H_{l}(\tau)}{p_{l}(\tau)}}$, where we have the probability distribution $q(\tau):=\left\lVert G(\tau)\right\rVert_{\infty}/\left\lVert G\right\rVert_{\infty,1}$.

Before presenting the proof, we explain how Theorem 8 can be applied to simulation in the LC model. We expect that the mixed-unitary channel $\sum_{l=1}^{L}\int_{0}^{t}\mathrm{d}\tau\,p_{l}(\tau)e^{-i\frac{H_{l}(\tau)}{p_{l}(\tau)}}\rho e^{i\frac{H_{l}(\tau)}{p_{l}(\tau)}}$ approximates the ideal evolution with $L^{1}$-norm scaling as in Theorem 5 and Theorem 7, but direct analysis would be considerably more complicated. However, universality (Statement 3 of Theorem 8) shows that this channel is the same as $\int_{0}^{t}\mathrm{d}\tau\,q(\tau)e^{-i\frac{G(\tau)}{q(\tau)}}\rho e^{i\frac{G(\tau)}{q(\tau)}}$. Thus, the analysis of Section 3.1 can be applied with the help of Theorem 8.

Proof of Theorem 8.

We define $G(\tau)$ to be the piecewise Hamiltonian

[TABLE]

where we use the abbreviation

[TABLE]

for the marginal probability distribution. Statements $1$ and $2$ can both be proved by directly evaluating the integrals

[TABLE]

and

[TABLE]

We use Statement $1$ to deduce that

[TABLE]

Therefore,

[TABLE]

which completes the proof of Statement $3$. ∎

Theorem 5′ ($L^{1}$-norm error bound for continuous qDRIFT (LC), short-time version).**

Let $H(\tau)=\sum_{l=1}^{L}H_{l}(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$ that is nonzero everywhere. Assume that each $H_{l}(\tau)$ is continuous and nonzero everywhere. Define \mathrm{E}(t,0)=\exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\big{)} and let $\mathcal{E}(t,0)(\cdot)=\mathrm{E}(t,0)(\cdot)\mathrm{E}^{\dagger}(t,0)$ be the corresponding quantum channel. Let $\mathcal{U}(t,0)$ be the continuous qDRIFT channel

[TABLE]

where $p_{l}(\tau)$ is the probability distribution $p_{l}(\tau):=\left\lVert H_{l}(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1,1}$. Then,

[TABLE]

In the special case where $H=\sum_{l=1}^{L}H_{l}$ is time independent, our bound reduces to

[TABLE]

This tightens a bound due to Campbell [13, Eq. (B12)] by a multiplicative factor from a tail bound. Note that [13] considered the distance $\left\lVert\mathcal{E}(t,0)-\mathcal{U}(t,0)\right\rVert_{\diamond}/2$, which is different from our definition of the diamond-norm distance $\left\lVert\mathcal{E}(t,0)-\mathcal{U}(t,0)\right\rVert_{\diamond}$.

Proof of Theorem 5*′*.

Consider the channel

[TABLE]

where $q(\tau):=\left\lVert G(\tau)\right\rVert_{\infty}/\left\lVert G\right\rVert_{\infty,1}$ and $G(\tau)$ is defined by (80). By Theorem 8, it suffices to bound $\left\lVert\mathcal{E}(t,0)-\mathcal{G}(t,0)\right\rVert_{\diamond}$.

Define two parametrized quantum channels

[TABLE]

and observe that

[TABLE]

For readability, we only consider the trace norm $\left\lVert\mathcal{E}_{1}(t,0)(\rho)-\mathcal{G}_{1}(t,0)(\rho)\right\rVert_{1}$, whose analysis can be easily adapted to bound $\left\lVert(\mathcal{E}_{1}(t,0)\otimes\mathds{1})(\sigma)-(\mathcal{G}_{1}(t,0)\otimes\mathds{1})(\sigma)\right\rVert_{1}$ and thus the diamond-norm distance $\left\lVert\mathcal{E}_{1}(t,0)-\mathcal{G}_{1}(t,0)\right\rVert_{\diamond}$.

By Lemma 6 and Theorem 8, we find that the first derivatives of $\mathcal{E}_{s}(t,0)(\rho)$ and $\mathcal{G}_{s}(t,0)(\rho)$ at $s=0$ agree with each other:

[TABLE]

Thus, we can apply the fundamental theorem of calculus twice and obtain

[TABLE]

which implies

[TABLE]

∎

Theorem 7′ ($L^{1}$-norm error bound for continuous qDRIFT (LC), long-time version).**

Let $H(\tau)=\sum_{l=1}^{L}H_{l}(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$ that is nonzero everywhere. Assume that each $H_{l}(\tau)$ is continuous and nonzero everywhere. Define \mathrm{E}(t,0)=\exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\big{)} and let $\mathcal{E}(t,0)(\cdot)=\mathrm{E}(t,0)(\cdot)\mathrm{E}^{\dagger}(t,0)$ be the corresponding quantum channel. Let $\mathcal{U}(t,0)$ be the continuous qDRIFT channel

[TABLE]

where $p_{l}(\tau)$ is the probability distribution $p_{l}(\tau):=\left\lVert H_{l}(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1,1}$. Then, for any positive integer $r$, there exists a division $0=t_{0}<t_{1}<\cdots<t_{r}=t$, such that

[TABLE]

To ensure that the simulation error is at most $\epsilon$, it thus suffices to choose

[TABLE]

The proof of Theorem 7*′* follows from Theorem 5*′* using the same reasoning as that used to prove Theorem 7.

3.3 Complexity of the continuous qDRIFT algorithm

As an immediate consequence of universality, we obtain the complexity of the continuous qDRIFT algorithm for simulating time-dependent Hamiltonians in both the SM and the LCU models.

Corollary 9** (Continuous qDRIFT algorithm with $L^{1}$-norm scaling (SM)).**

A $d$-sparse time-dependent Hamiltonian $H(\tau)$ acting on $n$ qubits can be simulated for time $\tau\in[0,t]$ with accuracy $\epsilon$ using

[TABLE]

queries to $\mathcal{O}_{\text{loc}}$, $\mathcal{O}_{\text{val}}$ and an additional

[TABLE]

gates, assuming that the probability distribution $p_{j}(\tau):=\left\lVert H(\tau)\right\rVert_{\max}/d^{2}\left\lVert H\right\rVert_{\max,1},\ j\in\{1,\ldots,d^{2}\}$ can be efficiently sampled.

Proof.

For any $\tau\in[0,t]$, Lemma 3 shows that $H(\tau)$ admits a decomposition $H(\tau)=\sum_{j=1}^{d^{2}}H_{j}(\tau)$, where each $H_{j}(\tau)$ is $1$-sparse and a query to any $H_{j}(\tau)$ can be simulated with $O(1)$ queries to $H(\tau)$. We use the continuous qDRIFT algorithm to simulate $H(\tau)=\sum_{j=1}^{d^{2}}H_{j}(\tau)$. We estimate

[TABLE]

where the second equality follows because $H_{j}(\tau)$ is $1$-sparse, and the inequality follows from Lemma 3. Assuming $\left\lVert H_{j}(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1,1}$ can be sampled efficiently, Theorem 7*′* implies that the algorithm has sample complexity and thus query complexity

[TABLE]

For each elementary exponentiation, we initialize a quantum register in the computational basis state $|\tau,j\rangle$ and use it to control the $1$-sparse term we need to simulate. This can be done with gate complexity \widetilde{O}\big{(}n\big{)}. Since the number of $1$-sparse simulations is the query complexity, we obtain the gate complexity

[TABLE]

as claimed.

Our above argument assumes that $\left\lVert H_{j}(\tau)\right\rVert_{\infty}$ is known a priori and that the distribution $\left\lVert H_{j}(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1,1}$ can be efficiently sampled. However, the argument still works if we replace each $\left\lVert H_{j}(\tau)\right\rVert_{\infty}$ by the upper bound

[TABLE]

which means we sample the distribution $p_{j}(\tau):=\left\lVert H(\tau)\right\rVert_{\max}/d^{2}\left\lVert H\right\rVert_{\max,1},\ j\in\{1,\ldots,d^{2}\}$. The claimed query and gate complexities follow from a similar analysis. ∎

Corollary 9′ (Continuous qDRIFT algorithm with $L^{1}$-norm scaling (LCU)).**

A time-dependent Hamiltonian with the LCU decomposition $H(\tau)=\sum_{l=1}^{L}\alpha_{l}(\tau)H_{l}$, where the controlled exponentiation of each $H_{l}$ can be performed with $g_{e}$ gates, can be simulated for time $\tau\in[0,t]$ with accuracy $\epsilon$ with gate complexity

[TABLE]

assuming that the probability distribution $p_{l}(\tau):=\alpha_{l}(\tau)/\left\lVert\alpha\right\rVert_{1,1}$ can be efficiently sampled.

Proof.

For any $H(\tau)=\sum_{l=1}^{L}\alpha_{l}(\tau)H_{l}$, we estimate

[TABLE]

The claimed complexity then follows from Theorem 7*′*. ∎

4 Rescaled Dyson-series algorithm

In this section, we propose a general principle for rescaling the Schrödinger equation (Section 4.1). We then apply this principle to improve the Dyson-series algorithm (Section 4.2) to achieve $L^{1}$-norm scaling.

4.1 A rescaling principle for the Schrödinger equation

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$. The evolution under $H(\tau)$ for time $t$ is given by the unitary operator \mathrm{E}(t,0)=\exp_{\mathcal{T}}\bigl{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\bigr{)}, which satisfies the Schrödinger equation

[TABLE]

We now propose a rescaling principle that helps to achieve $L^{1}$-norm scaling. The goal is to effectively have a Hamiltonian with constant spectral norm. Recall that for a time-independent Hamiltonian one can multiply the time by a constant and divide the Hamiltonian by the same constant and obtain the same time evolution. We can achieve something similar with a time-dependent Hamiltonian by rescaling the total evolution time to

[TABLE]

and using the rescaled Hamiltonian

[TABLE]

From this definition, it is obvious that the Hamiltonian has constant norm, because

[TABLE]

Moreover, we find that the time-evolution operator satisfies

[TABLE]

Solving this equation shows that we can obtain exactly the same time-evolution operator using the rescaled time and Hamiltonian:

[TABLE]

We also have the norm equality

[TABLE]

so any algorithm that simulates the rescaled Hamiltonian $\widetilde{H}(\varsigma)$ with complexity that scales with the $L^{\infty}$ norm can simulate the original Hamiltonian with $L^{1}$-norm scaling.

While our above discussion considers the spectral norm $\left\lVert\cdot\right\rVert_{\infty}$, other norms may be used depending on the input model of the Hamiltonian. Indeed, in the analysis for the SM model below we use the max-norm instead of the spectral norm.

Note that it may be hard in practice to compute the exact value of $\left\lVert H(\tau)\right\rVert$. However, we can instead use the change of variable

[TABLE]

where $\Lambda(\tau)\geq\left\lVert H(\tau)\right\rVert$ is any upper bound on the norm that can be efficiently computed.

4.2 Complexity of the rescaled Dyson-series algorithm

In this section, we show how the Dyson-series algorithm [8, 36, 29] can be rescaled to have $L^{1}$-norm scaling. We address this first for the SM model of Hamiltonian access before handling the LCU model (see Section 2.3 for definitions of these models).

Unlike continuous qDRIFT, the rescaled Dyson-series algorithm requires additional oracle access to the input Hamiltonian. Specifically, we need oracles that implement the inverse change-of-variable

[TABLE]

and compute the max-norm

[TABLE]

A quantum computer with access to these oracles can simulate time-dependent Hamiltonians with $L^{1}$-norm scaling. Note that because $f(\tau)$ increases monotonically, we can use binary search to compute $f^{-1}(\varsigma)$ up to precision $\delta$ using $O(\log(t/\delta))$ queries to $f$, so we expect it to be straightforward to implement the oracle $\mathcal{O}_{\rm var}$ in practice.

Theorem 10** (Rescaled Dyson-series algorithm with $L^{1}$-norm scaling (SM)).**

For $\tau\in[0,t]$, let $H(\tau)$ be a $d$-sparse Hamiltonian acting on $n$ qubits. Let $f(t):=\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\rVert_{\max}$, and suppose we have an upper bound on the max-norm, denoted \left\lVert H\big{(}f^{-1}(\varsigma)\big{)}\right\rVert_{\max}, that is positive and continuously differentiable. Then $H$ can be simulated for time $t$ with accuracy $\epsilon$ using

[TABLE]

queries to the oracles $\mathcal{O}_{\text{loc}}$, $\mathcal{O}_{\text{val}}$, $\mathcal{O}_{\text{var}}$, $\mathcal{O}_{\text{norm}}$ and an additional

[TABLE]

gates.

Proof.

We simulate the rescaled Hamiltonian \widetilde{H}(\varsigma):=H\big{(}f^{-1}(\varsigma)\big{)}/\left\lVert H\big{(}f^{-1}(\varsigma)\big{)}\right\rVert_{\max} for a total time of $\left\lVert H\right\rVert_{\max,1}:=\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\rVert_{\max}$ using the rescaling function

[TABLE]

Following [36, Theorem 9], we construct a unitary operation that block-encodes

[TABLE]

This construction is similar to [36, Lemma 8], except that the Hamiltonian is rescaled. Specifically, we use oracles $\mathcal{O}_{\text{var}}$ and $\mathcal{O}_{\text{norm}}$ to implement the transformation

[TABLE]

from which we obtain the rescaled Hamiltonian by querying $\mathcal{O}_{\text{val}}$ and re-normalizing the result with $\left\lVert H(f^{-1}(\varsigma))\right\rVert_{\max}$ to compute $\widetilde{H}_{jk}(\varsigma)$:

[TABLE]

We then uncompute the ancilla registers storing $f^{-1}(\varsigma)$ and $\left\lVert H(f^{-1}(\varsigma))\right\rVert_{\max}$. Overall, this implements a rescaled oracle

[TABLE]

The remaining algorithm proceeds as in [36]. As the implementation of each $\widetilde{\mathcal{O}}_{\text{val}}$ requires $O(1)$ queries to the oracles $\mathcal{O}_{\text{loc}}$, $\mathcal{O}_{\text{val}}$, $\mathcal{O}_{\text{var}}$, $\mathcal{O}_{\text{norm}}$, the overall query complexity is obtained by applying [36, Theorem 9] to the rescaled Hamiltonian, giving query complexity

[TABLE]

where

[TABLE]

Using this expression for $T$ gives the query complexity in (116).

We now analyze the gate complexity. If the entries of the Hamiltonian are given to within precision

[TABLE]

then the overall error due to the finite precision is $O(\epsilon)$. Since the maximum value of any matrix entry of $H$ is $\|H\|_{\max,\infty}$, the number of bits required is

[TABLE]

The implementation involves performing arithmetic on these values, which can be performed with complexity333In [29] and [36], the most complicated operations used are additions, which can be performed with complexity $O(n_{p})$. Here we are normalizing the Hamiltonian, so we must also perform multiplication and/or division, for which the straightforward approach has complexity $O(n_{p}^{2})$. While it is possible to perform multiplication and division with lower asymptotic complexity, such algorithms are only advantageous for very large instances, and do not affect the result as presented in Theorem 10, where logarithmic contributions to the gate complexity are suppressed. $O(n_{p}^{2})$. Since this is a logarithmic gate cost for each oracle query, it gives a contribution to the gate complexity of \widetilde{O}\bigl{(}d\left\lVert H\right\rVert_{\max,1}\bigr{)}.

The number of time steps is [36, Corollary 4]

[TABLE]

where $\alpha=d\|\widetilde{H}\|_{\max,\infty}$. The complexity to prepare the time registers is $\log M$ times the query complexity. We may ignore the second term in (127), because it is negligible compared to the complexity of the arithmetic.

We have

[TABLE]

Evaluating this derivative, we get

[TABLE]

so we obtain

[TABLE]

The gate complexity of the preparation of the time registers is $\log M$ times the query complexity, where we have shown that

[TABLE]

where $\|{\widetilde{H}^{\prime}}\|_{\infty,1}$ is polynomial in norms of $H$ and its derivative. Since this a logarithmic cost, the contribution to the complexity from preparation of the time registers is \widetilde{O}\bigl{(}d\left\lVert H\right\rVert_{\max,1}\bigr{)}.

The remaining contribution to the gate complexity comes from acting on the system itself. The cost of this is $O(n)$ for each of the oracle queries, which gives gate complexity \widetilde{O}\big{(}d\left\lVert H\right\rVert_{\max,1}n\big{)} (this is the dominant cost in (117)). ∎

Thus, the rescaled Dyson-series algorithm can simulate time-dependent Hamiltonians in the SM model with $L^{1}$-norm scaling. Next we turn our attention to the LCU model. For an input Hamiltonian $H(\tau)=\sum_{l=1}^{L}\alpha_{l}(\tau)H_{l}$, this approach assumes quantum access to the coefficient oracle

[TABLE]

in contrast to the continuous qDRIFT which only needs classical access. Given a classical circuit that computes the coefficients $\alpha_{l}(\tau)$, we can express it as a sequence of elementary gates and construct a corresponding quantum circuit with the same gate complexity. In our analysis, we ignore the implementation details and count the number of uses of the quantum oracle $\mathcal{O}_{\text{coeff}}$. The definitions of $\mathcal{O}_{\text{var}}$ and $\mathcal{O}_{\text{norm}}$ are similar to the SM case, except that the norm $\left\lVert H(\tau)\right\rVert_{\max}$ is replaced by $\left\lVert\alpha(\tau)\right\rVert_{\infty}$.

Theorem 10′ (Rescaled Dyson-series algorithm with $L^{1}$-norm scaling (LCU)).**

For $\tau\in[0,t]$, let $H$ be a time-dependent Hamiltonian with the decomposition $H(\tau)=\sum_{l=1}^{L}\alpha_{l}(\tau)H_{l}$, where each controlled $H_{l}$ can be performed with $g_{c}$ gates. Let $f(t):=\int_{0}^{t}\mathrm{d}\tau\left\lVert\alpha(\tau)\right\rVert_{\infty}$, and suppose we have an upper bound on the $\ell_{\infty}$ norm of the coefficients, denoted \left\lVert\alpha\big{(}f^{-1}(\varsigma)\big{)}\right\rVert_{\infty}, that is continuously differentiable. Then $H$ can be simulated for time $t$ with accuracy $\epsilon$ using

[TABLE]

queries to the oracles $\mathcal{O}_{\text{coeff}}$, $\mathcal{O}_{\text{var}}$, $\mathcal{O}_{\text{norm}}$ and an additional

[TABLE]

gates.

Proof.

We simulate the rescaled Hamiltonian \widetilde{H}(\varsigma):=H\big{(}f^{-1}(\varsigma)\big{)}/\left\lVert\alpha\big{(}f^{-1}(\varsigma)\big{)}\right\rVert_{\infty} for time $\left\lVert\alpha\right\rVert_{\infty,1}$ using the rescaling function

[TABLE]

The rescaled Hamiltonian takes the form

[TABLE]

where $\widetilde{\alpha}_{l}(\tau):=\alpha_{l}(\tau)/\left\lVert\alpha(\tau)\right\rVert_{\infty}$. Therefore, we have a global upper bound on the absolute value of the coefficients

[TABLE]

The remaining construction is similar to [29, Section V C], except that the Hamiltonian is rescaled. Specifically, we use oracles $\mathcal{O}_{\text{var}}$ and $\mathcal{O}_{\text{norm}}$ to implement the transformation

[TABLE]

from which we obtain the rescaled coefficients by querying $\mathcal{O}_{\text{coeff}}$ and doing arithmetic, giving

[TABLE]

We then uncompute the ancilla registers storing $f^{-1}(\varsigma)$ and $\left\lVert\alpha(f^{-1}(\varsigma))\right\rVert_{\infty}$. Overall, this implements a rescaled oracle

[TABLE]

The remaining algorithm proceeds as in [29]. As the implementation of each $\widetilde{\mathcal{O}}_{\text{coeff}}$ requires $O(1)$ queries to the oracles $\mathcal{O}_{\text{coeff}}$, $\mathcal{O}_{\text{var}}$, $\mathcal{O}_{\text{norm}}$, the overall query complexity is obtained by applying [29, Theorem 2] to the rescaled Hamiltonian. The analysis of the gate complexity proceeds along similar lines to that of Theorem 10. The multiplicative factor of $Lg_{c}$ is the cost of implementing the Select operation

[TABLE]

That complexity may be obtained from [18, Lemma G.7], or the unary iteration procedure of [2, Section III A]. ∎

5 Applications to chemistry and scattering theory

There are numerous cases in physics where one needs to simulate time-dependent quantum systems. Indeed, the pulse sequences that constitute individual quantum gates or adiabatic sweeps are described by time-dependent Hamiltonians. Here, we look at the particular case of simulating semi-classical scattering of molecules within a chemical reaction as an example of time-dependent Hamiltonian dynamics [46, 24].

Chemical scattering problems involve colliding reagents. As the molecules move closer, the electronic configuration changes due to strengthening Coulomb interactions, which is ultimately responsible for either the reagents forming a bond or flying apart depending on the initial conditions and the nature of the reagents. In the non-relativistic case, the Hamiltonian for two colliding atoms $A$ and $B$ at positions $x_{A}$ and $x_{B}$, respectively, and $M$ electrons with positions $x_{m}$ for $m=1,\ldots,M$, can be expressed as

[TABLE]

Here $p_{m}=[[p_{m}]_{x},[p_{m}]_{y},[p_{m}]_{z}]$ and $x_{m}=[[x_{m}]_{x},[x_{m}]_{y},[x_{m}]_{z}]$ are three-dimensional vectors of operators, whereas the corresponding nuclear terms (such as $x_{A}$ and $x_{B}$) are three-dimensional vectors of scalars. We further define $|x_{m}-x_{m}^{\prime}|$ to be the operator

[TABLE]

The wave function can be thought of as having a nuclear as well as an electronic component. First, we assume that the nuclear and the electronic wave functions are decoupled [26, 49]:

[TABLE]

This approximation is justified by the fact that the nuclear mass is substantially greater than the electronic mass. We then follow the time-dependent self-consistent field (TDSCF) approximation, which further treats $x_{A}$ and $x_{B}$ as classical degrees of freedom $x_{A}(t)$ and $x_{B}(t)$ with conjugate momenta ($p_{A}(t),p_{B}(t)$). This simplification is justified by Ehrenfest’s theorem, which states that for a sufficiently narrow quantum wave packet, the equation of motion for the centroid follows the classical trajectory (to leading order in $\hbar$). Under this approximation, the electronic dynamics satisfy

[TABLE]

where we have suppressed the implicit dependence of the electronic wave function on $x_{1},\ldots,x_{M}$. The equation of motion for the two nuclear positions in the time-dependent self-consistent field approximation within the Ehrenfest method is given by the Hamilton-Jacobi equation:

[TABLE]

and similarly for $x_{B}.$ The function $H_{\text{nuc}}(t)$ here is simply the Hamiltonian $H_{\text{nuc}}$ with the classical substitutions $x_{A}\rightarrow x_{A}(t)$, $p_{A}\rightarrow p_{A}(t)$ and similarly for $x_{B}$ and $p_{B}$. Similarly, we define $H_{\text{elec}}(t)$ to be the electronic Hamiltonian under this classical substitution.

The evolution in the Ehrenfest method is governed by a pair of tightly coupled quantum and classical dynamical equations, wherein the full Schrödinger equation only needs to be solved to understand part of the dynamics for the system. Indeed, as the Born-Oppenheimer approximation instantaneously holds under the above approximations, we can further express the electronic dynamics within a second-quantized framework with respect to a basis of molecular orbitals as

[TABLE]

where for some basis of orthonormal molecular orbitals ${\psi_{p}(\vec{x};t)}$ (which are implicitly time dependent if these basis functions are chosen to be functions of the nuclear positions, as would be appropriate for an atomic orbital basis),

[TABLE]

Thus under the above approximations, the dynamics that need to be simulated take the form of a standard second quantized simulation of chemistry, except the Hamiltonian is time dependent. The generalization of this to multiple nuclei is similarly straightforward, with the summation over two nuclear positions replaced by summation over all $L$ positions.

Consider the case where two reagents move towards each other from distant points with large momenta. To get an intuitive understanding of this evolution, it is instructive to examine the case of two molecules colliding using a classical force field. This will give us an expression that is qualitatively accurate for $x_{A}(t)$ and $x_{B}(t)$. To do this, we use a Lennard-Jones potential to model the interaction between two helium nuclei. The potential as a function of separation between the nuclei $r(t)=|x_{A}(t)-x_{B}(t)|$ is assumed to be of the form $V(r)=\epsilon\left(\frac{r_{m}^{12}}{r^{12}}-\frac{2r_{m}^{6}}{r^{6}}\right)$, where $\epsilon\approx 10$ K and $r_{m}\approx 2.6$ Å [44]. Setting the initial radial velocity to be approximately the root mean square (RMS) velocity of helium at $25\,\degree\mathrm{C}$, we solve the classical equations of motion to find the trajectory shown in Figure 1.

From Figure 1, we find that the interaction appears as a brief but intense kick between the two systems. As a result, the norm of the Hamiltonian changes dramatically throughout the evolution and we expect simulation algorithms with $L^{1}$-norm scaling to be advantageous over previous approaches. We leave a detailed study of such an advantage as a subject for future work.

6 Discussion

We have shown that a time-dependent Hamiltonian $H(\tau)$ can be simulated for the time interval $0\leq\tau\leq t$ with gate complexity that scales according to the $L^{1}$ norm $\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\rVert$. We designed new algorithms based on classical sampling and improved the previous Dyson-series approach to achieve this scaling. This is a polynomial speedup in terms of the norm dependence, an advantage that can be favorable in practice. In particular, our result has potential applications to simulating scattering processes in quantum chemistry. Our analysis also matches the intuition that the difficulty of simulating a quantum system should depend on the norm of the Hamiltonian instantaneously. This dual interpretation suggests that the $L^{1}$-norm dependence of our result cannot be significantly improved. However, further speedup might be possible if we know a priori the energy range of the initial state, as is suggested in [43, 34].

The rescaled Dyson-series approach is nearly optimal with respect to all parameters of interest. Indeed, a lower bound of \Omega\big{(}d\left\lVert H\right\rVert_{\max}t+\frac{\log(1/\epsilon)}{\log\log(1/\epsilon)}\big{)} was given in [9, Theorem 2] for simulating time-independent sparse Hamiltonians, which of course also holds for the more general time-dependent case. The query complexity (116) of the rescaled Dyson-series approach matches this dependence on $d\left\lVert H\right\rVert_{\max}t$ and on $\epsilon$, except that it scales as the product of the two terms instead of the sum (so, as in all quantum simulation algorithms prior to the advent of quantum signal processing [33], it does not achieve the optimal tradeoff between $t$ and $\epsilon$). However, this approach requires computing the rescaling function (114) and the Hamiltonian norm (115) in quantum superposition, which may introduce large overhead in practice. In comparison, continuous qDRIFT relies on classical sampling and may be better suited to near-term simulation. Its complexity has no dependence on the parameter $L$ in the LCU decomposition (Corollary 9*′*), which is advantageous for Hamiltonians consisting of many terms.

For most of our analysis, we have assumed that the Hamiltonian $H(\tau)$ is continuously differentiable. This assumption can be relaxed to allow finitely many discontinuities. Indeed, if $H(\tau)$ is discontinuous at the times $0=t_{0}<t_{1}<\cdots<t_{r}=t$ but otherwise continuously differentiable, we may divide the evolution into $r$ segments and apply a time-dependent Hamiltonian simulation algorithm within each time interval $[t_{j},t_{j+1}]$. For the Dyson-series approach, the complexity depends linearly on the $L^{1}$ norm, so concatenation gives a simulation of the entire evolution with $L^{1}$-norm scaling. The assumptions about the Hamiltonian can be even further relaxed: the continuous qDRIFT algorithm works properly provided only that $H(\tau)$ is Lebesgue integrable. Further discussion of this point is beyond the scope of this paper, and we refer the reader to [21] for details.

Our analysis can also be adapted to simulate time-dependent Hamiltonians that have countably many zeros. Indeed, since the equation $H(\tau)=0$ has at most countably many solutions, we can find $c\in\mathbb{R}$ such that $H(\tau)+cI$ is nonzero everywhere. Then, \exp_{\mathcal{T}}\bigl{(}-i\int_{0}^{t}\mathrm{d}\tau\,(H(\tau)+cI)\bigr{)}=e^{-ict}\exp_{\mathcal{T}}\bigl{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\bigr{)}, so the result is only off by a global phase. Note that this assumption can be completely dropped if we use continuous qDRIFT: we define the exceptional set

[TABLE]

and redefine $\mathcal{U}(t,0)$ as

[TABLE]

We note that $\mathcal{U}(t,0)$ is a valid quantum channel and can be implemented with unit cost. Indeed, for any input state $\rho$, we randomly sample a value $\tau$ according to $p(\tau)$ and perform $e^{-i{H(\tau)}/{p(\tau)}}$ if $\tau\in[0,t]\backslash\mathcal{B}_{0}$, and the identity operation otherwise. This implements

[TABLE]

The remaining analysis proceeds as in Section 3.

The qDRIFT protocol that we analyzed here only achieves first-order accuracy. It is natural to ask if sampling a different probability distribution could lead to an algorithm with better performance. The answer seems to be “no” if we only use a univariate distribution. To see this, consider the discrete case where $H=\sum_{l=1}^{L}H_{l}$ is a Hamiltonian consisting of $L$ terms. We sample according to a probability vector $p\in[0,1]^{L}$. Upon getting outcome $l$, we perform the unitary $e^{-itH_{l}/p_{l}}$. Effectively, we implement the quantum channel $\mathcal{U}(t)(\rho):=\sum_{l=1}^{L}p_{l}e^{-it\frac{H_{l}}{p_{l}}}\rho e^{it\frac{H_{l}}{p_{l}}}$, which is a first-order approximation to the ideal evolution $\mathcal{E}(t)(\rho):=e^{-it\sum_{l=1}^{L}H_{l}}\rho e^{it\sum_{l=1}^{L}H_{l}}$. In particular, the difference between $\mathcal{U}(t)(\rho)$ and $\mathcal{E}(t)(\rho)$ admits an integral representation

[TABLE]

To estimate the diamond-norm error $\left\lVert\mathcal{U}(t)-\mathcal{E}(t)\right\rVert_{\diamond}$, we take $\sigma$ to be a state on the joint system of the original register and an ancilla register with the same dimension. We compute

[TABLE]

By Jensen’s inequality,

[TABLE]

with equality if and only if all ${\left\lVert H_{l}\right\rVert_{\infty}}/{p_{l}}$ are equal, implying that the probability distribution $p_{l}:={\left\lVert H_{l}\right\rVert_{\infty}}/{\left\lVert H\right\rVert_{\infty,1}}$ is optimal. A similar optimality result holds for continuous qDRIFT (though the proof is more involved).

However, this does not preclude the existence of a higher-order qDRIFT protocol using more complicated sampling [39]. For example, besides the basic evolutions $e^{-itH_{l}/p_{l}}$, one could evolve under commutators $[H_{j},H_{k}]$ or anticommutators $\{H_{j},H_{k}\}$. We could also use a multivariate distribution and correlate different steps of the qDRIFT protocol. For future work, it would be interesting to find a higher-order protocol, or prove that such a protocol cannot exist.

The fractional-query algorithm described in Section 2.5 provides a natural approach to simulating time-dependent Hamiltonians whose query complexity scales with the $L^{1}$-norm. While we believe such a scaling also holds for the gate complexity, it would be highly nontrivial to give an explicit implementation. In any case, the fractional-query approach is streamlined by the Dyson-series approach and the latter can be rescaled to achieve $L^{1}$-norm scaling.

The rescaling principle that we proposed can potentially be applied to improve other quantum simulation algorithms. For example, we can use the product-formula algorithm [50] to simulate the rescaled Hamiltonian \widetilde{H}(\varsigma):=H\big{(}f^{-1}(\varsigma)\big{)}/\left\lVert H\big{(}f^{-1}(\varsigma)\big{)}\right\rVert_{\infty} for time $s=\left\lVert H\right\rVert_{\infty,1}$. The difficulty here is that the derivative of the rescaled Hamiltonian can be larger than the original one, making the rescaled algorithm perform worse. We leave a thorough study of this issue as a subject for future work.

Finally, it would be interesting to identify further applications of our $L^{1}$-norm scaling result, such as to designing new quantum algorithms and to improving the performance of quantum chemistry simulation. It might also be of interest to demonstrate these approaches experimentally, for applications such as implementing adiabatic algorithms with quantum circuits.

Acknowledgements

We thank Ryan Babbush, Earl Campbell, András Gilyén, Mária Kieferová, Robin Kothari, Tongyang Li, Guang Hao Low, Jarrod R. McClean, Yuval Sanders, and Leonard Wossnig for helpful discussions. DWB is funded by Australian Research Council Discovery Projects DP160102426 and DP190102633. This work was also supported in part by the Army Research Office (MURI award W911NF-16-1-0349), the Canadian Institute for Advanced Research, the National Science Foundation (grant 1813814), and the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum Algorithms Teams and Quantum Testbed Pathfinder programs. In addition, YS is supported by the Google Ph.D. Fellowship program. XW acknowledges support from the Department of Defense. NW is supported in part by a Google Quantum Research Award.

Appendix A Continuous qDRIFT and Hamiltonian averaging

Poulin, Qarry, Somma, and Verstraete developed an algorithm for time-dependent Hamiltonian simulation based on techniques of Hamiltonian averaging and Monte Carlo estimation [42]. In this section, we discuss the relation between their algorithm and our continuous qDRIFT.

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$. Assume that $H(\tau)$ is continuous, nonzero everywhere, and efficiently simulable for each particular $\tau$. Then, Poulin et al.’s approach simulates $H(\tau)$ for time $\tau\in[0,t]$ in two steps: (i) they replace the evolution \exp_{\mathcal{T}}\big{(}{-}i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\big{)} by an ordinary matrix exponential $e^{-itH_{\text{av}}}$ of the average Hamiltonian $H_{\text{av}}:=\frac{1}{t}\int_{0}^{t}\mathrm{d}\tau\,H(\tau)$ with an error that scales like O\big{(}(t\left\lVert H\right\rVert_{\infty,\infty})^{2}\big{)}; (ii) they further implement $\int_{0}^{t}\mathrm{d}\tau\,H(\tau)$ with Monte Carlo estimation by picking $m$ random times and approximating $\int_{0}^{t}\mathrm{d}\tau\ H(\tau)\approx\frac{1}{m}\sum_{k=1}^{m}H(\tau_{k})$ with error O\big{(}t\left\lVert H\right\rVert_{\infty,\infty}/\sqrt{m}\big{)}, the result of which is further approximated by product formulas.

The approach of [42] is essentially a sampling-based algorithm and thus similar in spirit to our continuous qDRIFT, except for a notable difference: their algorithm scales with the $L^{\infty}$ norm instead of the $L^{1}$ norm. Unfortunately, this drawback cannot be remedied merely by a better analysis of the same algorithm. Indeed, they use a uniform distribution to pick random times during the Monte Carlo estimation. This sampling ignores the instantaneous norm $\left\lVert H(\tau)\right\rVert_{\infty}$ of the Hamiltonian and therefore the resulting algorithm cannot scale with the $L^{1}$ norm $\int_{0}^{t}\mathrm{d}\tau\left\lVert H(\tau)\right\rVert_{\infty}$.

Instead, continuous qDRIFT uses a probability distribution that biases toward those times with larger instantaneous norm. In Section 3, we proved that such a sampling gives a direct simulation of time-dependent Hamiltonians with complexity that scales with the $L^{1}$ norm. We now give an indirect implementation: (i’) we show in Appendix A.1 that the error of replacing \exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\big{)} by an ordinary matrix exponential of $H_{\text{av}}$ scales like O\big{(}\left\lVert H\right\rVert_{\infty,1}^{2}\big{)}, improving the analysis of [42]; (ii’) we further prove in Appendix A.2 that the average Hamiltonian can be simulated by continuous qDRIFT with $L^{1}$-norm scaling. Combining these two steps, we see that the Monte Carlo estimation approach of [42] is superseded by continuous qDRIFT.

A.1 Hamiltonian averaging

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$ and assume that it is continuous and nonzero everywhere. Define

[TABLE]

where $H_{\text{av}}:=\frac{1}{t}\int_{0}^{t}\mathrm{d}\tau\,H(\tau)$ is the average Hamiltonian. Our goal is to bound the distance between $\mathrm{E}(s,0)$ and $\mathrm{E}_{\text{av}}(s)$ at $s=t$. Using the initial condition $\mathrm{E}(0,0)=\mathrm{E}_{\text{av}}(0)=I$, we have

[TABLE]

By the Schrödinger equation

[TABLE]

we obtain

[TABLE]

which implies, by telescoping, that

[TABLE]

By the fundamental theorem of calculus, the first term of the integrand can be bounded as

[TABLE]

To handle the second term, we use Lemma 1. Observe that the generator of $\mathrm{E}(s,0)$ is $H(u),\ 0\leq u\leq s$, whereas the generator of $\mathrm{E}(\tau,0)$ is $H(u),\ 0\leq u\leq\tau$. So they only differ on the interval \big{[}\min\{s,\tau\},\max\{s,\tau\}\big{]}. Consequently,

[TABLE]

Altogether, we have

[TABLE]

Theorem 11** (Hamiltonian simulation by averaging (spectral-norm distance)).**

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$ and assume that it is continuous and nonzero everywhere. Define \mathrm{E}(t,0):=\exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau\,H(\tau)\big{)} and $\mathrm{E}_{\text{av}}(t):=e^{-itH_{\text{av}}}$, where $H_{\text{av}}:=\frac{1}{t}\int_{0}^{t}\mathrm{d}\tau\,H(\tau)$ is the average Hamiltonian. Then,

[TABLE]

The above bound on the spectral-norm error can be converted to a bound on the diamond-norm error using Lemma 2.

Theorem 11′ (Hamiltonian simulation by averaging (diamond-norm distance)).**

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$ and assume that it is continuous and nonzero everywhere. Define unitary operators \mathrm{E}(t,0):=\exp_{\mathcal{T}}\big{(}-i\int_{0}^{t}\mathrm{d}\tau\ H(\tau)\big{)}, $\mathrm{E}_{\text{av}}(t):=e^{-itH_{\text{av}}}$ and let $\mathcal{E}(t,0)(\cdot):=\mathrm{E}(t,0)(\cdot)\mathrm{E}(t,0)^{\dagger}$, $\mathcal{E}_{\text{av}}(t)(\cdot)=\mathrm{E}_{\text{av}}(t)(\cdot)\mathrm{E}_{\text{av}}(t)^{\dagger}$ be the corresponding channels. Then,

[TABLE]

A.2 Implementing Hamiltonian averaging by continuous qDRIFT

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$ and assume that it is continuous and nonzero everywhere. We have showed that the ideal evolution can be approximated by an evolution under the average Hamiltonian with error that scales with the $L^{1}$ norm. We now show that such a Hamiltonian averaging can be implemented by continuous qDRIFT, again with $L^{1}$-norm scaling. This improves over the algorithm of [42] which scales with the $L^{\infty}$ norm.

Theorem 12** (Hamiltonian averaging by continuous qDRIFT).**

Let $H(\tau)$ be a time-dependent Hamiltonian defined for $0\leq\tau\leq t$ and assume that it is continuous and nonzero everywhere. Define $\mathrm{E}_{\text{av}}(t):=e^{-itH_{\text{av}}}$ and let $\mathcal{E}_{\text{av}}(t)(\cdot)=\mathrm{E}_{\text{av}}(t)(\cdot)\mathrm{E}_{\text{av}}(t)^{\dagger}$ be the corresponding channels. Let $\mathcal{U}(t,0)$ be the continuous qDRIFT channel

[TABLE]

where $p(\tau)=\left\lVert H(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1}$. Then,

[TABLE]

Note that by applying the triangle inequality to Theorem 5 and Theorem 11*′*, we obtain

[TABLE]

Theorem 12 improves the constant prefactor from $8$ to $4$.

Proof of Theorem 12.

We parametrize the two channels $\mathcal{E}_{\text{av}}(t)$, $\mathcal{U}(t)$ and define

[TABLE]

Since $\mathcal{E}_{\text{av},0}(t)(\rho)=\rho$, $\mathcal{E}_{\text{av},1}(t)(\rho)=\mathcal{E}_{\text{av}}(t)(\rho)$, $\mathcal{U}_{0}(t,0)(\rho)=\rho$, and $\mathcal{U}_{1}(t,0)(\rho)=\mathcal{U}(\rho)$, the first derivative of $\mathcal{E}_{\text{av},u}(t)(\rho)$ and $\mathcal{U}_{u}(t,0)(\rho)$ agrees with each other at $u=0$

[TABLE]

Applying the fundamental theorem of calculus twice, we obtain

[TABLE]

We take $\sigma$ to be a state on the joint system of the original register and an ancilla register with the same dimension. Using properties of the Schatten norms, we have

[TABLE]

Using the definition $p(\tau)=\left\lVert H(\tau)\right\rVert_{\infty}/\left\lVert H\right\rVert_{\infty,1}$, the second term of the integrand can be further simplified as

[TABLE]

giving

[TABLE]

Optimizing over $\sigma$ proves the claimed bound. ∎