1 Introduction↩︎

The quantum simulation problem [1], [2] is central to quantum computing, with applications in quantum chemistry, materials science, and condensed matter physics (see, e.g., [3]–[6]). An important special case is the Hamiltonian simulation problem, which concerns the time evolution of a finite-dimensional quantum system governed by the Schrödinger equation \[\label{eqn:HamSim} i \frac{d}{dt} u(t) = H(t) u(t), \quad t \in [0,T],\tag{1}\] where \(u(t) \in \mathbb{C}^N\) and \(H(t) \in \mathbb{C}^{N\times N}\) is Hermitian for each \(t \in [0,T]\). Since 1 describes a finite-dimensional closed quantum system, its propagator is a unitary operator. We therefore refer to 1 as describing unitary dynamics. In this setting, efficient quantum algorithms have been developed for both time-independent [7]–[17] and time-dependent [18]–[23] cases. However, these methods do not directly apply to general linear non-homogeneous differential equations of the form \[\label{non-unitary-ode} \frac{d}{dt} u(t) = -A(t)u(t) + b(t),\tag{2}\] where \(A(t) \in \mathbb{C}^{N \times N}\) may be non-anti-Hermitian and \(b(t) \neq 0\). Such equations arise in broader quantum simulation settings, notably in the simulation of finite-dimensional open quantum systems. We therefore refer to 2 as describing non-unitary dynamics.

Several quantum algorithms have been developed for linear differential equations of the form 2 (see, e.g., [24]–[31]. Most of the existing approaches are based on the Quantum Linear System Algorithm (QLSA) [32]–[35], typically embedding \(A(t)\) into a larger dilated matrix and solving an associated enlarged linear system with additional implementation and complexity overhead. Recently, An et al. [36], [37] and Low and Somma [38] introduced and analyzed the Linear Combination of Hamiltonian Simulation (LCHS) framework, which expresses the solution of 2 as an integral over unitary evolutions. This approach simulates non-unitary dynamics by combining Hamiltonian simulation algorithms developed for unitary dynamics. The resulting algorithms achieve optimal state-preparation cost and optimal scaling in both the evolution time and the target precision.

From a broader perspective, LCHS is a special case of the Linear Combination of Unitaries (LCU) framework of Childs and Wiebe [39], as it represents non-unitary dynamics through linear combinations of unitary simulation primitives. A related LCU-based construction is the Multi-Product Formula (MPF) [14], [16], [17], [23], combining Lie–Trotter–Suzuki product formulas [7] to obtain high-order Hamiltonian simulation algorithms. Specifically, MPFs construct these high-order approximations by linearly combining sequences of low-order Trotter steps with carefully chosen weights. This systematic cancellation of algorithmic errors circumvents the deep exponential nesting required by standard high-order product formulas, achieving near-optimal time and precision dependence while maintaining favorable commutator scaling from product formulas [15]. However, for general non-unitary dynamics, existing approaches, including LCHS and other advanced algorithms, cannot currently achieve both of these desirable properties simultaneously. This limitation naturally motivates the central question of our work: can we design an efficient quantum algorithm for general non-unitary systems with both near-optimal scaling (in time and precision) and commutator scaling?

In this work, we answer this question affirmatively in a quadrature-dependent sense. We develop an LCHS–MPF algorithm for non-unitary dynamics that applies to both time-independent and time-dependent settings, and illustrate its use in applications arising from linear differential equations and open quantum systems. In this algorithm, LCHS reduces non-unitary simulation to a family of Hamiltonian simulations, while MPFs implement these simulations with commutator-sensitive complexity. Our key observation is that the outer LCHS quadrature does not merely control discretization error. Rather, the quadrature nodes also determine the Hamiltonians passed to the inner simulation routine, and therefore affect the commutator structure, step-size restrictions, LCU normalization, weighted error accumulation, and final query complexity. We track this dependence explicitly and obtain a quantum algorithm for non-unitary simulation whose complexity reflects both the quadrature structure of LCHS and the commutator scaling of MPFs.

1.1 Summary of Results↩︎

We briefly summarize our approach and main results in the time-independent homogeneous case. The non-homogeneous and time-dependent cases are discussed in 9. Let \(A=L+iH\), where \(L\succeq0\) and \(L,H\) are Hermitian. The LCHS representation expresses the dissipative evolution \(e^{-AT}\) as an integral over unitary evolutions \[\label{lchs-intro} e^{-AT}\approx \frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}} \hat{f}(k)e^{-iG_k T}\,dk.\tag{3}\] Here \(k\in\mathbb{R}\) is the frequency variable, \(\hat{f}\) is the kernel function and \(G_k:=H+kL\) is the Hamiltonian simulated at frequency \(k\). After applying a quadrature rule \(Q=\{(k_i,w_i):i\in\mathcal{I}_Q\}\) with nodes \(k_i\) and weights \(w_i\), the integral is replaced by a finite sum of the form \[W_Q^{\operatorname{ideal}} = \sum_{i\in\mathcal{I}_Q}v_i e^{-iG_{k_i}T}, \qquad v_i = \frac{w_i\hat{f}(k_i)}{\sqrt{2\pi}}.\] Thus, the quadrature rule determines not only the discretization error but also the Hamiltonians that must be simulated. We replace each \(e^{-iG_{k_i}T}\) with an MPF approximation. The post-quadrature error and complexity bounds are then determined by the quadrature radius \(R_Q := \max_{i\in\mathcal{I}_Q}|k_i|\), which controls the range of Hamiltonians \(H+kL\) appearing in the simulation, and by the LCU normalization factor \(\alpha_Q := \sum_{i\in\mathcal{I}_Q}|v_i|\), which controls the success probability and amplitude-amplification overhead of the outer LCU procedure. Since LCHS is combined with Hamiltonian simulation via MPF, the cost also depends on a commutator scale \(\mu_Q\), which quantifies the MPF simulation difficulty of the quadrature-indexed Hamiltonian family \(\{G_{k_i}\}_{i\in \mathcal{I}_Q}\). Our first result is a post-quadrature complexity estimate that makes this quadrature-dependent commutator scaling explicit. The following result is stated in a simplified informal form. More precise complexity estimates are given in 26 and 27.

The role of \(\alpha_Q\) in 1 is important. In previous LCHS algorithms, the kernel and quadrature rule are typically fixed, so the corresponding LCU normalization is either bounded explicitly or absorbed into constants. In contrast, we keep \(\alpha_Q\) explicit because the kernel and quadrature choices also affect the commutator scale. Consequently, its effect must be evaluated jointly with \(\alpha_Q\) and \(\mu_Q\). We remark that \(\alpha_Q\) and \(\mu_Q\) may also scale as \(\mathcal{O}(\operatorname{polylog}(1/\epsilon))\) for certain quadrature rules. We will state these scalings explicitly in the later discussions.

In 1, \(F(m)\) and \(\mu_Q\) are auxiliary functions that upper bound \(\mu_{m,Q}\), which is defined in 24 , by \(\mu_{m,Q} \leq F(m) \mu_Q\). Here, we only note that \(\mu_{m,Q}\) depends on the quantities \[\Lambda_{m,Q} := \sum_{i\in\mathcal{I}_Q}|v_i|\Phi_m(k_i,N) \qquad \Phi_{m,Q}^{\ast} := \max_{i\in\mathcal{I}_Q}\Phi_m(k_i,N),\] where \(m\) denotes the MPF order and \(\Phi_m(k,N)\) is an abstract MPF error profile for simulating \(G_k\). Hence, 1 is abstract because it is phrased in terms of a general MPF error profile. In [spec-aftab] [spec-mizuta], we specialize this result to two concrete settings, which we summarize below:

1. In [16], Aftab et al. provided the first complete error analysis for the time-independent MPF. Informally, their analysis associates to each pair \((H,L)\) a commutator scale, denoted here by \(\chi_J(H,L)\), which measures how rapidly the relevant nested commutators can grow. For a quadrature rule \(Q\), this scale must be compared with the quadrature radius \(R_Q\). Therefore, we choose \(\rho_Q\) satisfying \[0<\rho_Q<\frac{\chi_J(H,L)}{(J+1)^{1/J}\max\{1,R_Q\}}.\] In this case, we have \[\begin{align} \Phi_{m,\nu}(H,L) :={}& \sum_{\substack{j\in 2\mathbb{Z}_+\\ j\ge 2m}} \sum_{l=1}^{m} \frac{\rho_Q^{\,j+l-(2m+1)}}{l!} \sum_{\substack{j_1,\dots,j_l\in 2\mathbb{Z}_+\\ j_1+\cdots+j_l=j}} \sum_{\substack{\ell_1,\dots,\ell_l\ge0\\ 0\le \ell_\kappa\le j_\kappa+1\\ \ell_1+\cdots+\ell_l=\nu}} \prod_{\kappa=1}^{l} \mathcal{C}_{j_\kappa+1,\ell_\kappa}(H,L) \end{align}\] for \(\nu\ge0\). Here, \(\mathcal{C}_{j_\kappa+1,\ell_\kappa}(H,L)\) is determined by nested commutators of length \(j_\kappa+1\) in \(H\) and \(L\) that contain exactly \(\ell_\kappa\) occurrences of \(L\). If the discrete moments are \(M_\nu^Q := \sum_{i\in\mathcal{I}_Q}|v_i|\,|k_i|^\nu\) then \(\Lambda_{m,Q} := \sum_{\nu\ge0}\Phi_{m,\nu}(H,L)M_\nu^Q\). In this case, \(F(m) \equiv 1\).

2. In [17], Mizuta proved an improved MPF estimate of Aftab et al. for local and extensive Hamiltonians. We use this estimate when it applies to \(\{G_{k_i}\}_{i\in\mathcal{I}_Q}\). For \(p_0\in\mathbb{N}\), define \[\mu_{m,p_0}(G_{k_i}) := \sup_{\substack{j,l\in\mathbb{N}\\ j\ge 2m,\;1\le l\le \lfloor j/2\rfloor}} \left( \sum_{\substack{2\le j_1,\ldots,j_l\le p_0-1\\ j_1+\cdots+j_l=j}} \prod_{\kappa=1}^{l} \alpha_{\operatorname{comm},j_\kappa+1}(G_{k_i}) \right)^{\frac{1}{j+l}},\] where \(\alpha_{\operatorname{comm},j}(G_{k_i})\) denotes the \(j\)-th order nested-commutator. The locality-based estimate considers the quantities \(\Lambda_{m,Q}(p_0) := \sum_{i\in\mathcal{I}_Q} |v_i|\, \mu_{m,p_0}(G_{k_i})^{2m+1}\) and \(\Phi_{m,Q}^{\ast}(p_0) := \max_{i\in\mathcal{I}_Q} \mu_{m,p_0}(G_{k_i})^{2m+1}\). In this case, \(F(m) = \mathcal{O}(\log^2 m)\).

1 compares representative algorithms in the time-independent homogeneous setting. Compared with optimal LCHS, our algorithm replaces the norm-based Hamiltonian simulation cost with a commutator-sensitive MPF cost after quadrature. Thus, the method is useful when the family \(H+kL\) has favorable commutator structure. The complexity also depends on quadrature-induced quantities such as \(R_Q\) and \(\alpha_Q\). Hence, the method is not a black-box improvement over LCHS. Rather, it is a structure-sensitive refinement. The tradeoff is that the dependence on \(T\) and \(1/\epsilon\) is not linear and logarithmic, respectively, as in optimal LCHS. Therefore, our result should be viewed as near-optimal in regimes where the commutator scaling also enters the complexity analysis favorably.

Our second result concerns the choice of quadrature rule. Because the quadrature rule determines the nodes \(k_i\) and weights \(v_i\), it affects both the discretization error and the MPF cost through the quantities \(R_Q\), \(\alpha_Q\), the discrete moments \(M_\nu^Q\), and the resulting weighted MPF profiles. We therefore compare quadrature rules in terms of these quantities, together with the quadrature error \(E_{\operatorname{quad}}(Q)\) and the node count \(|\mathcal{I}_Q|\).

Thus, the sinh–sinh quadrature rule improves the quadrature-cardinality scaling due to its double-exponential decay rate. This reduction in the number of quadrature nodes directly improves the LCU implementation cost of the algorithm.

We further demonstrate the utility of the LCHS-MPF framework through three concrete classes of applications.

1. Fractional diffusion equations with imaginary potentials, including fractional Bloch–Torrey-type dynamics, model anomalous diffusion in complex media and arise in transport and magnetic-resonance-type applications [40]–[44]. After discretization, \(A=L+iH\), where \(L\) comes from the fractional diffusion operator and \(H\) from the potential, and the LCHS–MPF analysis for \(G_k=H+kL\) gives \[\mu_Q=\mathcal{O} (A_V^{(N)} (1+R_QN^{\theta_s/d} ) ).\] Here \(N\) is the number of grid points, \(d\) is the spatial dimension, \(0<s<1\) is the fractional order, \(\theta_s=\max\{0,2s-1\}\), \(R_Q\) is the quadrature radius, and \(A_V^{(N)}\) measures the relevant discrete Fourier regularity of the potential. Hence, the complexity reflects the diffusion–potential commutator structure rather than only a worst-case norm of \(A\).

2. Advection–diffusion equations model heat transfer, mass transport, pollutant dispersion, and related transport phenomena [45]–[48]. In the discretized system, diffusion gives the dissipative part and advection gives the Hamiltonian part, and for \(G_k=H+kL\) the commutator scale is \[\mu_Q=\mathcal{O} (|b|N^{1/d}+aR_QN^{2/d} ).\] Here \(a>0\) is the diffusion coefficient, \(b\) is the advection strength, \(N\) is the number of grid points and \(d\) is the dimension. This separates the effects of diffusion, advection, discretization, and quadrature, giving a sharper estimate than standard norm-based LCHS bounds.

3. No-jump evolution in dissipative transverse-field Ising systems with local spontaneous emission arises in quantum-trajectory descriptions of open many-body systems [49]–[52]. The effective non-Hermitian Hamiltonian yields \(\dot{\psi}=-(L+iH)\psi\), the quadrature Hamiltonians \(G_k=H+kL\) remain \(2\)-local, and in the large-system regime the leading commutator scale is \[\mu_Q=\mathcal{O} \left( \alpha_Q^{3/2}|V|^{1/2}\left(|h|+D|J|+\frac{R_Q\gamma}{2}\right)^{3/2}\right).\] Here \(|V|\) is the number of spins, \(h\) is the transverse-field strength, \(J\) is the Ising coupling, \(D\) is the maximum graph degree and \(\gamma\) is the local emission rate. Thus, the complexity is controlled by Pauli commutator structure rather than only by a generic norm or block-encoding normalization, exposing local many-body structure via commutator scaling.

1.2 Related Works↩︎

Several recent works study when non-unitary quantum simulation can inherit structural advantages from unitary Hamiltonian simulation. Wang et al. [53] recently derived commutator-based Trotter error bounds for Lindbladian simulation. Their work is close in motivation to ours. Both analyses ask whether the favorable commutator scaling of Hamiltonian product formulas persists in dissipative dynamics. The technical settings, however, are distinct. Their bounds apply directly to product-formula approximations of Lindblad generators and use Richardson extrapolation for observable estimation, whereas our approach first uses the LCHS representation to reduce dissipative linear evolution to a quadrature-indexed family of Hamiltonian simulations. Thus, the commutators appearing in our analysis are those of the Hamiltonians generated by the LCHS quadrature, and the main issue is to track how this quadrature dependence enters the MPF error and complexity estimates.

A second related direction concerns quantum algorithms for non-unitary dynamics based on matrix-function and eigenvalue-transformation representations. An et al. [54] extended the LCHS approach from matrix exponentials to broader matrix eigenvalue transformations represented through Laplace transforms. Low and Su [55] introduced quantum eigenvalue processing, an eigenvalue transformation framework for non-normal matrices that addresses related non-Hermitian phenomena from a different algorithmic perspective. These works complement the present paper by broadening the class of matrix functions that can be implemented on quantum computers. By contrast, the focus of the present work is not to enlarge the class of admissible matrix functions, but to analyze how the post-quadrature LCHS representation interacts with a commutator-sensitive Hamiltonian simulation primitive.

Contour-integral methods provide another point of comparison. Takahira et al. [56] used Cauchy’s integral formula and block encoding to implement matrix functions as linear combinations of shifted matrix inverses. More recently, Jiang and An [57] analyzed the complexity of contour-integral-based quantum eigenvalue transformations, with applications to Hamiltonian simulation, matrix polynomials, and linear differential equations. Wang et al. [58] developed a contour-based matrix-decomposition framework for non-unitary dynamics using a finite decomposition derived from Cauchy’s residue theorem. These approaches are related to LCHS through analytic matrix-function representations, However, their elementary operations are typically shifted resolvents or Hermitian transformations, rather than Hamiltonian evolutions of the form \(e^{-i(H+kL)t}\) whose commutator structure is tracked in the present work.

Other transform-based approaches have also been developed in the literature. Schrödingerization [30], [59]–[61] converts non-unitary linear dynamics into unitary dynamics in one higher dimension, achieving optimal or near-optimal matrix-query dependence after suitable smooth initialization. Jin et al. [62] proposed a transmutation-based method for dissipative diffusion generated by positive semi-definite operators, expressing the diffusion semigroup via the Kannai transform as a Gaussian-weighted superposition of unitary wave propagators. A recent Poisson-summation framework by Wang et al. [63] further emphasizes the role of spectral aliasing and discretization in non-unitary matrix transformations. These works illustrate that non-unitary simulation can often be reduced to superpositions of unitary or resolvent-type primitives.

1.3 Discussion and Open Questions↩︎

Our main result shows that the query complexity is polylogarithmic in time and inverse precision, and at the same time achieves commutator scaling. This matches our previous analysis of commutator scaling for MPF. It remains an open question whether one can obtain optimal scaling in both time and inverse precision while preserving commutator scaling.

Another natural question is whether the quadrature dependence can be further optimized beyond the quadrature rules analyzed here. Although the present analysis accommodates fairly general quadrature rules, the resulting complexity bounds still depend on quantities such as \(\Lambda_{m,Q}\), \(R_Q\), \(\alpha_Q\), and the discrete moments \(M^Q_\nu\). It would be interesting to determine whether there exist quadrature rules that simultaneously improve discretization accuracy, LCU normalization, moment growth, and the post-quadrature MPF cost. A related question is how to optimize the LCHS kernel profile jointly with the quadrature rule. The kernel parameters influence the approximation error, truncation error, quadrature radius, LCU normalization, and the discrete moment bounds entering the MPF error profiles. A more systematic kernel–quadrature optimization may therefore improve the overall complexity of the algorithm.

The LCU implementation introduces an overhead governed by the normalization of the linear combination and by the norm of the desired output state. Although amplitude amplification can boost the success probability, the resulting state-preparation cost still contains normalization-dependent factors. It would be interesting to determine whether alternative block-encoding constructions, quantum signal processing techniques, or dissipative embeddings can reduce this overhead while preserving the commutator-sensitive advantages of the LCHS–MPF framework.

Several extensions of the time-dependent framework remain open. The analysis in 9.2 assumes sufficient regularity of the time-dependent generator on each subinterval. It would be interesting to determine whether similar complexity guarantees can be obtained under weaker assumptions, such as piecewise-smooth, discontinuous, or rapidly varying dissipative generators within the LCHS-compatible setting.

Finally, it would be useful to better understand the practical performance of MPF-based dissipative simulation. The present work focuses primarily on asymptotic complexity bounds, but concrete implementations may exhibit additional structure that is not visible at the worst-case level. Numerical studies could clarify which parameter regimes are most favorable for near-term or fault-tolerant quantum architectures.

1.4 Organization↩︎

The remainder of the paper is organized as follows. 2 reviews the notation and algorithmic preliminaries. 3 analyzes the pre-quadrature approximation, truncation, and inner-simulation errors for the time-independent homogeneous case. 4 develops the post-quadrature LCHS–MPF error and complexity analysis, including the abstract theorem, commutator-scaling specializations, and comparison of quadrature rules. 5 discusses various applications.

Acknowledgments↩︎

Junaid Aftab acknowledges the support by the National Science Foundation under the grant DMS-2231533. Dong An acknowledges funding from Quantum Science and Technology - National Science and Technology Major Project via Project 2024ZD0301900, and the support by The Fundamental Research Funds for the Central Universities, Peking University. Konstantina Trivisa acknowledges support from the National Science Foundation under grants DMS-2008568 and DMS-2231533.

2 Preliminaries↩︎

This section collects the preliminary details used throughout the paper. 2.1 introduces the notation, 2.2 reviews product formulas, 2.3 recalls the Linear Combination of Unitaries (LCU) framework, 2.4 describes Multi-Product Formulas (MPFs), and 2.5 introduces the Linear Combination of Hamiltonian Simulation (LCHS) framework.

2.1 Notation↩︎

We briefly discuss the main notation used throughout this work.

2.1.1 Standard Notation↩︎

Let \(\mathbb{R}\), \(\mathbb{C}\), and \(\mathbb{N}\) denote the sets of real, complex and natural numbers, respectively. We write \(\mathbb{Z}\) for the integers and \(\mathbb{Z}_{+}\) for the positive integers. For \(M \in \mathbb{N}\), let \(\mathbb{Z}^d_M := \{0,\ldots,M-1\}^d\) be the \(d\)-fold Cartesian product of the cyclic group, \(\mathbb{Z}_M\). All logarithms are natural logarithms unless otherwise stated. If \(f, g : \mathbb{N} \to \mathbb{R}^+\) are non-negative functions, we use the following standard notation from complexity theory:

1. \(f(n) = \mathcal{O}(g(n))\) if and only if there exists a constant \(C > 0\) and an integer \(M \in \mathbb{N}\) such that \(f(n) \leq C g(n)\) for all \(n \geq M\). The notation \(\widetilde{\mathcal{O}}\) hides polylogarithmic factors.

2. \(f(n) = \Omega(g(n))\) if and only if there exists a constant \(C > 0\) and an integer \(M \in \mathbb{N}\) such that \(f(n) \geq C g(n)\) for all \(n \geq M\).

The symbol \(\overleftarrow{\prod_{\gamma}}\) denotes a product in which the elements are arranged with increasing indices from right to left. Similarly, the notation \(\overrightarrow{\prod_{\gamma}}\) is defined to denote a product with elements ordered with increasing indices from left to right. Specifically, we have \[\overset{\longleftarrow}{\prod_{i=1,\cdots,k}} A_i = A_k \cdots A_1, \quad \overset{\longrightarrow}{\prod_{i=1,\cdots,k}} A_i = A_1 \cdots A_k.\] \(\| f \|_{L^1[a,b]}\) denotes the \(L^1\)-norm, defined as \(\|f\|_{L^1[a,b]} = \int_a^b |f(s)| \, ds\), on \([a,b]\). We use the convention \(\hat{f}(\xi) = \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} f(x) e^{-i \xi x} dx\) and \(f(x) = \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} \hat{f}(\xi) e^{i \xi x} d\xi\) for the Fourier transform and the inverse Fourier transform, respectively.

2.1.2 Linear Algebra Notation↩︎

All vector spaces considered are finite-dimensional and defined over \(\mathbb{C}\). Vectors are denoted either by lowercase Roman letters with a \(\vec{\cdot}\) symbol on top or by ket notation, depending on context. Linear operators—represented as matrices—are denoted by uppercase Roman letters. The linear algebra notation used throughout this paper is summarized below:

1. \(A^\dagger\) denotes the conjugate transpose of \(A\) and \(I\) denotes the identity matrix of appropriate dimension.

2. \(\| \vec{a} \|_1\) denotes the \(1\)-norm of a vector \(\vec{a}\), defined as \(\| \vec{a} \|_1 = \sum_{i=1}^n |a_i|\).

3. \(\| A \|\) denotes the spectral norm of a matrix \(A\), defined as its largest singular value.

4. \(\| A \|_{L^1[a,b]}\) denotes the \(L^1\)-norm of a time dependent matrix \(A(t)\) on \([a,b]\), defined as \(\int_a^b \|A(s)\| \, ds\).

5. The notation \(A \succeq 0\) denotes that \(A\) is a positive semi-definite matrix.

6. For two matrices \(A\) and \(B\), the commutator is defined by \([A,B] = AB - BA\). More generally, for matrices \(A_1, A_2, \ldots, A_m\), we define the nested commutator \[[A_1, A_2, \ldots, A_m] := [A_1, [A_2, \ldots, [A_{m-1}, A_m] \cdots ]].\] In the special case where \(A_1 = \cdots = A_{m-1} = A\) and \(A_m = B\), this nested commutator is denoted by \(\operatorname{ad}_A^{\,m-1}(B)\).

The formal solution to the matrix-valued linear differential equation \(\frac{dU(t)}{dt} = -i H(t) U(t)\) for \(T_0 \leq t \leq T\) is written as \(U(T_0,T) = \exp_{\mathcal{T}} ( -i \int_{T_0}^T H(s) ds )\), where \(\exp_{\mathcal{T}}\) denotes the time-ordered exponential operator. We also denote the time-ordered exponential operator as \(\mathcal{T}e\). The reader is referred to [64] for more details. In what follows, we set \(T_0 = 0\) and write \(U(T) := U(0,T)\) when convenient.

2.2 Product Formulas↩︎

Consider a time-dependent Hamiltonian \(H(t)\) for \(t \in [0,T]\), decomposed as \(H(t)=\sum_{\gamma=1}^{\Gamma} H_{\gamma}(t)\). For each \(\gamma\), let \(U_{\gamma}(T)\) denote the time-ordered exponential generated by \(H_{\gamma}(t)\) on \([0,T]\). Product-formula methods assume that each \(U_{\gamma}(T)\) can be implemented efficiently. Under this assumption, a generic product formula approximates the exact time-evolution operator \(U(T)\) as \[\label{generic-prod} U_p(T) = \prod_{\xi=1}^{\Xi} \prod_{\gamma=1}^{\Gamma} U_{\pi_\xi(\gamma)}(\beta_{\xi \gamma} T,(\beta_{\xi \gamma} + \alpha_{\xi \gamma})T),\tag{4}\] where \(\pi_\xi\) is a permutation in \(S_{\Gamma}\) and \(\alpha_{\xi,\gamma},\beta_{\xi,\gamma}\in[0,1]\) such that \(\beta_{\xi \gamma} + \alpha_{\xi \gamma} \in [0,1]\). We say that \(U_p\) is a \(p\)th-order product formula if its single-step approximation error satisfies \[\|U_p(T)-U(T)\| = \mathcal{O}(T^{p+1}).\] In practice, the interval \([0,T]\) is partitioned into \(r\) segments, and the evolution on each segment is approximated by a product formula. The parameter \(r\) is referred to as the Trotter number. The resulting approximation is \((U_p(T/r))^r\).

When \(H(t) \equiv H\) is time-independent, the error scaling of product formulas is well understood. In particular, [15] showed that, for a general \(p\)th-order product formula, ensuring the global error bound \(\lVert U_p(T/r)^r - U(T) \rVert \le \epsilon\) requires choosing \(r\) such that \[\label{product-error} r = \mathcal{O} \left( \frac{\alpha_{\operatorname{comm},p+1}(H)^{1/p} T^{1+1/p}}{\epsilon^{1/p}} \right), \quad \alpha_{\operatorname{comm},p+1}(H) = \sum_{j_1,\dots,j_{p+1}=1}^{\Gamma} \|[H_{j_1},\dots,H_{j_{p+1}}]\|,\tag{5}\] where \(\epsilon\) is the target precision. The dependence on \(\alpha_{\operatorname{comm},p+1}(H)\) is commonly referred to as commutator scaling. This property enables product-formula Hamiltonian simulation algorithms to achieve favorable system-size scaling for physically local Hamiltonians (see, e.g., [67]–[69]). By contrast, the error scaling of product formulas for time-dependent Hamiltonians is less completely understood. Early work of Wiebe et al. [18] derived bounds for smoothly varying Hamiltonians based on \(1\)-norm scaling. Subsequent studies [20], [70], [71] analyzed low-order time-dependent product formulas, some of which exhibit commutator scaling. A general treatment was later given by Childs et al. [72], who established commutator scaling for smoothly varying Hamiltonians with finite-range interactions. More recently, Mizuta et al. [23] analyzed generic smoothly varying time-dependent product formulas and proved commutator-scaling bounds. In this work, we use only their results on time-dependent multi-product formulas, which are recalled in 2.4.

2.3 Linear Combination of Unitaries↩︎

The Linear Combination of Unitaries (LCU) algorithm [39] provides a method for implementing linear combinations of unitary operators on a quantum computer. Let \(M \in \mathbb{N}\), let \(a_1,\ldots,a_M \in \mathbb{C}\), and let \(U_1,\ldots,U_M\) be unitary operators. The LCU algorithm aims to implement the operator \(U = \sum_{j=1}^{M} a_j U_j\). The algorithm assumes that two operators can be efficiently implemented. First, it requires access to state preparation (\(\operatorname{PREP}\)) oracles: \[\begin{align} \label{F-LCU} \operatorname{PREP}_{\operatorname{R}} \ket{0} = \frac{1}{\sqrt{\| \vec{a} \|_1 }} \sum_{j=1}^{M} \sqrt{a_{j}} \ket{j}, \quad \operatorname{PREP}_{\operatorname{L}} \ket{0} = \frac{1}{\sqrt{\| \vec{a} \|_1 }} \sum_{j=1}^{M} \overline{\sqrt{a_{j}}} \ket{j} \end{align}\tag{6}\] where \(\sqrt{\cdot}\) denotes the principal branch of the square root. Second, it assumes access to a multi-qubit controlled operation, commonly called a select oracle, of the form \[\label{controlled-LCU} \operatorname{SEL}=\sum_{j=1}^{M} U_j \otimes \ket{j}\bra{j}.\tag{7}\] Assuming access to the required operators, the LCU algorithm implements the unitary operator \(U_{\text{LCU}} = (I \otimes \operatorname{PREP}^\dagger_{\operatorname{L}}) \cdot \operatorname{SEL} \cdot (I \otimes \operatorname{PREP}_{\operatorname{R}})\). If \(\ket{\psi}\) denotes an input state, we have \[U_{\text{LCU}} \ket{\psi} \ket{0}^{\log M} = \frac{1}{\| \vec{a} \|_1} \left(\sum_{j=1}^{M} a_j U_j\right) \ket{\psi} \ket{0}^{\log M} + \ket{\perp},\] where \(\ket{\perp}\) denotes a potentially non-normalized state satisfying \((I \otimes \ket{0}^{\log M}\bra{0}^{\log M}) \ket{\perp} = 0\). Constructing \(U_{\operatorname{LCU}}\) requires a single query to each of the oracles \(\operatorname{PREP}^\dagger_{\operatorname{L}}\), \(\operatorname{PREP}_{\operatorname{R}}\), and \(\operatorname{SEL}\), with the implementation of \(\operatorname{SEL}\) typically dominating the computational cost. The ancilla register, initially prepared in the state \(\ket{0}^{\log M}\), must be measured to obtain \(\sum_{j=1}^{M} a_j U_j\ket{\psi}\). To achieve a constant success probability, \(\mathcal{O}(\|\vec{a}\|_1)\) rounds of amplitude amplification [73] are required.

2.4 Multi-Product Formula↩︎

The Multi-Product Formula (MPF) for time-independent Hamiltonian simulation constructs a linear combination of time-independent product formulas \[\label{multiproduct} U_{\operatorname{MP}}(T) = \sum_{j=1}^{M} a_j \, U^{b_j}_{2}(T/b_j), \quad a_j\in \mathbb{R}, b_j \in \mathbb{N}.\tag{8}\] Low, Kliuchnikov, and Wiebe [14] showed that, for \(M = m\), the coefficients \(a_j\) and exponents \(b_j\) can be chosen such that \(\|\vec{a}\|_1 = \mathcal{O}(\log m)\) and \(\|\vec{b}\|_1 = \mathcal{O}(m^2 \log m)\), ensuring that \(\|U_{\operatorname{MP}}(T) - e^{-i H T} \| = \mathcal{O}(T^{2m+1})\). Consequently, although an MPF is constructed from product formulas of fixed order, it can achieve arbitrarily high orders of convergence. It can be implemented using the LCU algorithm as follows:

1. Consider implementing the MPF over a short time-step, \(\Delta\). Define the state-preparation oracles as in 6 , and let the select oracle \[\operatorname{SEL} = \sum_{j=0}^{m-1} U_2\left(\Delta/b_{j+1}\right)^{b_{j+1}} \otimes \ket{j}\bra{j}\] encode the MPF. Then \((I \otimes \operatorname{PREP}^\dagger_{\operatorname{L}}) \cdot \operatorname{SEL} \cdot (I \otimes \operatorname{PREP}_{\operatorname{R}})\) implements the MPF via LCU, with the ancilla state \(\ket{0}\) encoding the target operation and initial success probability \(\mathcal{O}(1/\|\vec{a}\|_1^2)\).

2. Consider simulating evolution over total time \(T\). Partition \([0,T]\) into \(r\) equal segments. On each segment, apply the MPF for \(\Delta=T/r\) using LCU. Robust oblivious amplitude amplification is applied at each step to maintain a constant success probability.

The cost of the implementation above is dominated by the oracle \(\operatorname{SEL}\), which uses \(\mathcal{O}(\|\vec{b}\|_1)\) queries to controlled-\(U_2\) operations. Amplitude amplification contributes an additional factor of \(\mathcal{O}(\|\vec{a}\|_1)\). Thus, each time step requires \(\mathcal{O}(\|\vec{a}\|_1\|\vec{b}\|_1)\) queries to \(U_2\), and \(r\) steps require \(\mathcal{O}(r\|\vec{a}\|_1\|\vec{b}\|_1)\) queries in total. Although \(\|\vec{a}\|_1\) and \(\|\vec{b}\|_1\) are specified in 2.4, a rigorous choice of \(r\) was previously unavailable and was expected to exhibit commutator scaling because the construction uses \(U_2\). This gap was recently closed by Aftab et al. [16] and Mizuta [17]. We recall the relevant short-time simulation results below, with time-step size denoted by \(\Delta\).

While 4 applies to general time-independent Hamiltonians, it has the limitation that the infimum is taken over commutators of arbitrarily large depth. If \(\alpha_{\operatorname{comm},j}=\Omega(N^{j^\gamma})\) for some \(\gamma>1\), then this infimum vanishes, and no positive admissible \(\Delta\) exists. Mizuta [17] addressed this issue for \(q\)-local, \(g\)-extensive Hamiltonians by showing that only commutators up to finite depth need to be considered in the short-time MPF bound. If \(H=\sum_{\gamma=1}^{\Gamma} H_\gamma\), let \(\gamma\), let \(X_\gamma\) denote the support of \(H_\gamma\). We say that \(H\) is \(q\)-local if \(|X_\gamma|\le q\) for all \(\gamma\) and \([H_\gamma,H_{\gamma'}]=0\) whenever \(X_\gamma\cap X_{\gamma'}=\varnothing\). Moreover, we say that \(H\) is \(g\)-extensive if \(\sum_{\gamma:\, j\in X_\gamma}\|H_\gamma\| \le g\) for all \(j\).

The Multi-Product Formula for time-dependent Hamiltonian simulation, constructed from the second-order time-dependent product formula, is given by \[\label{time-dep-2MPF} U_{\operatorname{MP}}(T) = \sum_{j=1}^{m} a_j \overset{\longleftarrow}{\prod_{b=1,\cdots,b_j}} U_{2}\left( \frac{(b - 1)T}{b_j}, \frac{bT}{b_j} \right).\tag{9}\] An error analysis was recently given by Mizuta et al. [23], who showed that time-dependent MPFs also exhibit commutator scaling. In particular, they proved the following bound for MPFs constructed from second-order time-dependent Trotter–Suzuki formulas.

2.5 Linear Combination of Hamiltonian Simulation↩︎

The Linear Combination of Hamiltonian Simulation (LCHS) algorithm simulates non-unitary dynamics by expressing them as linear combinations of unitary dynamics. LCHS formulation was first introduced in [36] and improved in [37]. This approach generalizes the Fourier representation of the exponential function on the positive real axis, \[\label{fourier-one} e^{-x}=\int_{\mathbb{R}} \hat{f}(k) e^{-ikx} \, dk, \quad x \geq 0.\tag{10}\] Specifically, for a possibly time-dependent operator \(A(t) = L(t) + i H(t)\), with Hermitian \(L(t) \succeq 0\) and \(H(t)\), they show that \[\mathcal{T}e^{-\int_0^t A(s)\,ds} =\int_{\mathbb{R}} \hat{f}(k) U_k(t)dk, \quad t \in [0,T]\] solves 2 in the case \(b(t) \equiv 0\). The case \(b(t) \neq 0\) is then derived via an application of Duhamel’s principle. The LCHS algorithm then proceeds by truncating the integral to \([-K,K]\) for some \(K > 0\) and implementing it via the LCU algorithm, where each \[\label{Uk} U_k(t) := \mathcal{T} e^{-i \int_0^t (H(s) + k L(s)) \, ds}\tag{11}\] is implemented via a time-dependent Hamiltonian simulation algorithm. The function \(\hat{f}(k)\) is called the LCHS kernel function and should be the inverse Fourier transform of the function \(f(x)\) which takes the value of \(e^{-x}\) for \(x \geq 0\). Previous works have suggested the choices of \(\hat{f}(k)\) to be \(\frac{1}{\pi(1+k^2)}\) [36] or \(\frac{1}{2 \pi e^{-2^\beta} e^{(1+iz)^\beta}}, \beta \in (0,1)\) [37], where the LCHS algorithm with the latter kernel function can solve 2 in near-optimal time and precision complexity. Recently, Low and Somma [38] introduced an approximate LCHS framework that further improves upon the previous formalism. Their key observation is that the requirement that 10 hold exactly can be relaxed so that it need only hold approximately. That is, it suffices to find a function \(f\) such that, for \(\epsilon > 0\), we have \(|e^{-x} - f(x)| \le \epsilon\) for \(x \ge 0\), and such that 10 holds only approximately. This weaker condition admits a broader class of kernels and leads to a flexible and more efficient LCHS algorithm with optimal complexity. Their main result is summarized below.

In what follows, we use the improved approximate formulation of Low and Somma [38].

3 Pre-Quadrature Error Analysis↩︎

In this section, we give a pre-quadrature error analysis of the algorithm in the time-independent homogeneous case. Throughout, let \(A=L+iH\), where \(L\) and \(H\) are Hermitian and \(L\succeq 0\). In the pre-quadrature regime, the algorithm has three main sources of error:

1. (Approximation) The approximation error, \(E_{\mathrm{approx}}(y_0)\), arising by replacing \(e^{-AT}\) with the LCHS representation. That is, \(E_{\mathrm{approx}}(y_0) := \| \frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\hat{f}(k)U_k(T)\,dk - e^{-AT} \|\). By ?? , after taking the limit \(R\to\infty\), it is bounded by \(\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}|\hat{f}(k-iy_0)|\,dk\).

2. (Truncation) The truncation error, \(E_{\mathrm{trunc}}(R)\), arising from restricting \(\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\hat{f}(k)\,U_k(t)\,dk\) to \([-R,R]\). That is, \(E_{\mathrm{trunc}}(R) :=\| \frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\hat{f}(k)U_k(T)\,dk - \frac{1}{\sqrt{2\pi}}\int_{-R}^{R}\hat{f}(k)U_k(T)\,dk \|\). By ?? , since \(U_k(T)\) is unitary, it is bounded by \(\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}\setminus[-R,R]}|\hat{f}(k)|\,dk\).

3. (Inner Simulation) The simulation error \(E_{\mathrm{MPF}}(T)\), arising from replacing each \(U_k(T)\) by its MPF approximation \(\widetilde{U}_k(T)\).

In what follows, we use the four-parameter kernel family introduced by Low and Somma [38]: \[\label{low-somma-kernel} \hat{f}_{a,b}(k;c,d) = \frac{(b+1)^{a-1}}{\sqrt{2\pi}} \frac{e^{d(1-ik)} e^{-\frac{k^2+1}{4c^2}}}{(1-ik)(b+ik)^{a-1}}, \quad a \in \mathbb{N},\; b,c>0,\; d\in\mathbb{R}.\tag{12}\] We call \(\vec{\theta}=(a,b,c,d)\) the kernel profile and write \(\hat{f}(k)\) for \(\hat{f}_{a,b}(k;c,d)\) when clear. 3.1 discusses the approximation error, 3.2 the truncation error, 3.3 the inner simulation error, and 3.4 discusses the coupling between the truncation and inner simulation errors.

3.1 Approximation Error↩︎

We first analyze the approximation error \(E_{\mathrm{approx}}(y_0)\). Low and Somma analyzed \(E_{\mathrm{approx}}(y_0)\) for the subclass \(\hat{f}_{2,1}(k;c,d)\). Here, we analyze \(E_{\mathrm{approx}}(y_0)\) for the full family \(\hat{f}_{a,b}(k;c,d)\). Most proofs are relegated to 6.

The proof of 8 can be found in 6.1. The integral in 8 has an exact representation in terms of the confluent hypergeometric function, but that form is less transparent. We instead use the bound \(((b+y_0)^2+k^2)^{(a-1)/2}\geq (b+y_0)^{a-1}\) to obtain the simpler upper bound \[\label{bessel-integ} E_{\mathrm{approx}}(y_0) \le \frac{(b+1)^{a-1}}{(b+y_0)^{a-1}}\frac{e^{\,d(1-y_0)+\frac{y_0^2-1}{4c^2}}}{2\pi}\int_{\mathbb{R}}\frac{e^{-k^2/(4c^2)}}{\sqrt{(y_0-1)^2+k^2}}\,dk.\tag{13}\] This bound is explicit and retains the essential dependence on \(y_0\) and \(a,b,c,d\). 9, which is proved in 6.2, evaluates this integral form.

Using 9 with \(s=y_0-1\), \(\sigma=c\) and \(K_0(x)\leq \sqrt{\frac{\pi}{2x}}e^{-x}\) for \(x>0\), we obtain \[\label{surrogate} E_{\mathrm{approx}}(y_0) \le \frac{c}{\sqrt{\pi}}\frac{(b+1)^{a-1}}{(y_0-1)(b+y_0)^{a-1}}e^{d(1-y_0)+\frac{y_0^2-1}{4c^2}}.\tag{14}\] 14 motivates minimizing the right-hand side and identifying kernel profiles for which this bound is at most \(\epsilon_{\operatorname{approx}}\). In 10, which is proved in 6.3 , we first study existence and uniqueness of minimizers with respect to \(y_0\).

10 implies that, for a fixed kernel profile \(\vec{\theta}=(a,b,c,d)\), the surrogate cost \(B(y_0)\) admits a unique minimizer \(y_0^\ast\). If \(\epsilon_{\operatorname{approx}}>0\) denotes the approximation error budget, then \[\label{suff-from-surr} \text{there exists } y_0 \in (1,\infty) \text{ such that } \frac{c}{\sqrt{\pi}}B(y_0)\le \epsilon_{\operatorname{approx}} \iff \frac{c}{\sqrt{\pi}}B(y_0^\ast)\le \epsilon_{\operatorname{approx}}.\tag{15}\] 15 gives a sufficient condition for \(E_{\operatorname{approx}}(y_0)\leq \epsilon_{\operatorname{approx}}\). It remains to determine when \(y_0^\ast\) satisfies \(\frac{c}{\sqrt{\pi}}B(y_0^\ast)\leq \epsilon_{\operatorname{approx}}\). 11 gives a condition in terms of the parameter \(d\).

11 shows that taking \(d\) above an explicit lower bound ensures \(E_{\mathrm{approx}}(y_0)\leq \epsilon_{\mathrm{approx}}\). The threshold is determined by the unique solution of \(\Psi(s)=\epsilon_{\mathrm{approx}}\), which can be computed numerically. 12, which is proved in 6.4, gives explicit sufficient conditions on \(c\) and \(d\) and determines a corresponding choice of \(y_0\).

3.2 Truncation Error↩︎

We now compute the truncation error in 13.

The proof of 13 can be founnd in 6.5. The bound in 13 shows that \(E_{\mathrm{trunc}}(R)\to 0\) as \(R\to\infty\). This motivates considering the set of radii for which the truncation error is below a fixed tolerance \(\epsilon'>0\) \[\label{g-def} \mathcal{H}_{\epsilon'}:=\{R>0:E_{\mathrm{trunc}}(R)<\epsilon'\}.\tag{16}\] 14 gives a computable lower bound ensuring \(R\in\mathcal{H}_{\epsilon'}\).

3.3 Inner Simulation Error↩︎

We next compute the inner simulation error. For \(k\in\mathbb{R}\), define \(G_k:=H+kL\). Since \(G_k\) consists of \(H\) and \(L\), each nested commutator is formed by choosing either \(H\) or \(L\) at each position, motivating the following definition.

16 formalizes the fact that the nested commutators of \(G_k\) depend polynomially on \(k\).

We approximate \(U_k(t)=e^{-iG_kt}\) by a \(2m\)-th order MPF built from the second-order product formula, denoted \(U_k^{(2m)}\). We state the error bound abstractly, which is instantiated later using either 4 or 5.

We now quantify how the abstract error bound in 1 propagates through the LCHS representation. We first define the quantities needed for the short-time analysis in 18, prove a one-step MPF error bound for step size \(\Delta\) in 19, and extend it to \([0,T]\) in 20. Since the argument does not use specific properties of \(\hat{f}_{a,b}(k;c,d)\), we write \(\hat{f}(k)\) unless the parameter dependence is needed.

We first establish the error bound for a short time step \(\Delta\) in 19. The proof can be found in 7.1.

We now extend this short-time error bound to the full simulation over \([0,T]\) in 20.

3.4 Coupling of Truncation and Inner Simulation Errors↩︎

The key point of 3.3 is that, after integrating over \(k\in[-R,R]\), the inner simulation error depends on \(\Lambda_m(R,N)\). Hence, \(R\) and the inner simulation error are coupled. Increasing \(R\) reduces truncation error but, when \(\Phi_m(k,N)\) grows with \(|k|\), also increases \(\Lambda_m(R,N)\) and amplifies simulation error. We capture this trade-off through a feasibility problem. Define

\[\label{varf-def} \mathcal{F}_{\epsilon_{\operatorname{comb}}} := \left\{ (R,\delta,r)\in(0,\infty)^2\times\mathbb{N} : E_{\mathrm{trunc}}(R) + E_{\mathrm{MPF}}(T) \le \epsilon_{\operatorname{comb}}, \frac{T}{r} \le \Delta(R,N,\delta) \right\}.\tag{19}\] Here \(\epsilon_{\operatorname{comb}}\) is the total error budget for the truncation and inner simulation errors. If \(E_{\mathrm{trunc}}(R) \geq \epsilon_{\mathrm{comb}}\), then \((R,\delta,r)\notin \mathcal{F}_{\epsilon_{\mathrm{comb}}}\) for every \(r\in\mathbb{N}\) and \(\delta>0\). Hence a positive error budget for the inner simulation is available only when \(R\in\mathcal{H}_{\epsilon_{\mathrm{comb}}}\). For such \(R\), define \(\epsilon_R := \epsilon_{\mathrm{comb}}-E_{\mathrm{trunc}}(R)>0\). After fixing \(R \in \mathcal{H}_{\epsilon_{\mathrm{comb}}}\), the problem therefore reduces to determining pairs \((r,\delta)\in\mathbb{N}\times(0,\infty)\) such that \[E_{\mathrm{MPF}}(T) \le \epsilon_R, \quad \frac{T}{r} \le \Delta(R,N,\delta) \label{eq:feasibility-reduced}\tag{20}\] The choice of \(r\) depends on whether 4 or 5 instantiates 1. We defer this discussion to 4.

4 Post-Quadrature Complexity Analysis↩︎

We now analyze the overall query complexity of our algorithm. We fix an admissible kernel from 12 . 4.1 gives the post-quadrature error bound, 4.2 describes the implementation, 4.3 derives complexity under the abstract MPF estimate in 1, and [spec-aftab] [spec-mizuta] specialize the analysis to 4 and 5, respectively. 4.6 discusses various quadrature rules.

4.1 Post-Quadrature Error Bound↩︎

The analysis in 3.4 was carried out before discretizing the LCHS integral. We now show that the same coupling persists after applying a quadrature rule. Let \(Q=\{(k_i,w_i)\}_{i\in\mathcal{I}_Q}\) be a quadrature rule, with \(k_i\in\mathbb{R}\) and \(w_i\in\mathbb{C}\). Define its radius by \(R_Q:=\max_{i\in\mathcal{I}_Q}|k_i|\). The corresponding ideal quadrature operator is \[\label{eq:generic-quadrature-operator} W_Q^{\mathrm{ideal}} := \sum_{i\in\mathcal{I}_Q}v_iU_{k_i}(T), \qquad v_i:= \frac{w_i}{\sqrt{2\pi}}\hat{f}_{a,b}(k_i;c,d)\tag{21}\] and the associated quadrature mass is \(\alpha_Q := \sum_{i\in\mathcal{I}_Q}|v_i|\). Under 1, we also define the discrete weighted MPF profile by \[\label{eq:Lambda-Q} \Lambda_{m,Q} := \sum_{i\in\mathcal{I}_Q}|v_i|\Phi_m(k_i,N).\tag{22}\] This is the post-quadrature analogue of \(\Lambda_m(R,N)\). Define \(\Phi_{m,Q}^\ast:=\max_{i\in\mathcal{I}_Q}\Phi_m(k_i,N)\). Let \(U_{k_i}^{(2m)}(T)\) denote the \(2m\)-th order MPF approximation to \(U_{k_i}(T)\) over \([0,T]\), and define \(W_Q^{\mathrm{MPF}} := \sum_{i\in\mathcal{I}_Q}v_i U^{(2m)}_{k_i}(T)\). Specific quadrature rules are discussed in 4.6. We have the following discrete analog of 20.

21 shows that, after fixing the quadrature rule, the inner simulation error couples to the quadrature error through \(R_Q\), \(\alpha_Q\), and \(\Lambda_{m,Q}\). This is the post-quadrature analogue of 3.4. We end with a total error bound, which is the starting point for the complexity analysis in 4.3.

22 shows that the quadrature rule affects the final error through \(E_{\mathrm{quad}}(Q)\) and through \(R_Q\), \(\alpha_Q\), and \(\Lambda_{m,Q}\). Thus, quadrature couples discretization error to quantum implementation error. This yields the post-quadrature problem. Choose an admissible rule \(Q\) with \(E_{\mathrm{quad}}(Q)<\epsilon_{\mathrm{comb}}\), leaving residual budget \(\epsilon_Q := \epsilon_{\mathrm{comb}}-E_{\mathrm{quad}}(Q) > 0\). Then choose \((r,\delta)\in\mathbb{N}\times(0,\infty)\) such that \[(1+\eta_Q(T/r))^{r-1} \left( C_m\frac{T^{2m+1}}{r^{2m}}\Lambda_{m,Q} + r\alpha_QR_m(\delta ) \right) \le \epsilon_Q, \quad \frac{T}{r}\le \Delta(R_Q,N,\delta ).\] This is the post-quadrature analogue of the discussion in 3.4.

4.2 Quantum Implementation↩︎

The complexity estimates below address two tasks in a fully quantum implementation. The first task is to construct a block-encoding of the post-quadrature approximation \[W_Q^{\operatorname{MPF}} = \sum_{i \in \mathcal{I}_Q} v_i U^{(2m)}_{k_i} (T) = \sum_{i \in \mathcal{I}_Q} v_i \left( \sum_{j=1}^m a_j \left[ U_{2,k_i}\left(\frac{T}{r b_j}\right) \right]^{b_j} \right)^r .\] Here \(v_i\) denotes the quadrature coefficients, while \(a_j,b_j\) denote the coefficients from 2.4. Given \(Q\), define the quadrature preparation oracles \[\operatorname{PREP}^Q_{\operatorname{R}}\ket{0} = \frac{1}{\sqrt{\alpha_Q}} \sum_{i\in\mathcal{I}_Q} \sqrt{v_i}\ket{i}, \quad \operatorname{PREP}^{Q}_{\operatorname{L}}\ket{0} = \frac{1}{\sqrt{\alpha_Q}} \sum_{i\in\mathcal{I}_Q} \overline{\sqrt{v_i}}\ket{i},\] where \(\alpha_Q:=\sum_{i\in\mathcal{I}_Q}|v_i|\). For each \(i\in\mathcal{I}_Q\), let \(\mathcal{U}_Q\) be a unitary block-encoding of \(U_{k_i}^{(2m)}(T)\), normalized to one after amplitude amplification at each of the \(r\) time steps. Hence, we have \[(\bra{0} \otimes I)\mathcal{U}_Q (\ket{0} \otimes I) = U_{k_i}^{(2m)}(T).\] Define the controlled unitary \(\operatorname{SELECT}_Q := \sum_{i\in\mathcal{I}_Q} \ket{i}\bra{i}\otimes \mathcal{U}_Q\). By the outer LCU construction, we obtain a block-encoding of \(W_Q^{\operatorname{MPF}}\) \[(\bra{0} \bra{0} \otimes I) (\operatorname{PREP}_{\operatorname L}^{Q\dagger}\otimes I) \operatorname{SELECT}_Q (\operatorname{PREP}_{\operatorname R}^{Q}\otimes I) (\ket{0} \ket{0} \otimes I) = \frac{W_Q^{\operatorname{MPF}}}{\alpha_Q}.\] Writing \(K_m:=\|\vec{a}\|_1\|\vec{b}\|_1\), one controlled implementation of \(U_{k_i}^{(2m)}(T)\) with \(r\) time steps costs \(\mathcal{O}(rK_m)\) controlled second-order product-formula queries. Hence, one use of \(\operatorname{SELECT}_Q\) has the same cost. One block-encoding of \(W_Q^{\operatorname{MPF}}/\alpha_Q\) uses \(\operatorname{SELECT}_Q\), the quadrature preparation oracles \(\operatorname{PREP}_{\operatorname L}^{Q}\), \(\operatorname{PREP}_{\operatorname R}^{Q}\), and their required inverses in the outer LCU construction.

The second task is to prepare the normalized output state \(W_Q^{\operatorname{MPF}}u(0)/\|W_Q^{\operatorname{MPF}}u(0)\|\). Write \(\ket{\widetilde{u}(0)}:=u(0)/\|u(0)\|\). Applying the block-encoding of \(W_Q^{\operatorname{MPF}}/\alpha_Q\) to \(\ket{\widetilde{u}(0)}\) and post-selecting the ancillas on \(\ket{0}\) gives the un-normalized state \(W_Q^{\operatorname{MPF}}\ket{\widetilde{u}(0)}/\alpha_Q\), whose norm is \[\frac{\|W_Q^{\operatorname{MPF}}\ket{\widetilde{u}(0)}\|}{\alpha_Q} = \frac{\|W_Q^{\operatorname{MPF}}u(0)\|}{\alpha_Q\|u(0)\|}.\] Hence, amplitude amplification requires \(\chi_Q(u(0)) := \frac{\alpha_Q\|u(0)\|}{\|W_Q^{\operatorname{MPF}}u(0)\|}\) uses of the block-encoding, up to absolute constants. Therefore, normalized state preparation costs \(\mathcal{O}(\chi_Q(u(0))rK_m)\) controlled second-order product-formula queries, \(\mathcal{O}(\chi_Q(u(0)))\) uses of the quadrature preparation oracles and their inverses, and \(\mathcal{O}(\chi_Q(u(0)))\) uses of the state-preparation oracle for \(\ket{\widetilde{u}(0)}\) and its inverse.

4.3 Overall Query Complexity↩︎

Split the error tolerance as \(\epsilon_{\operatorname{approx}}+\epsilon_{\operatorname{comb}}\leq\epsilon\), where \(\epsilon_{\operatorname{approx}}\) controls \(E_{\operatorname{approx}}(y_0)\) and \(\epsilon_{\operatorname{comb}}\) controls the combined quadrature and inner-simulation errors. 11 12 gives sufficient conditions on \(c,d\) to ensure the approximation error is bounded blow by \(\epsilon_{\operatorname{approx}}\). Choose a quadrature rule \(Q\) such that \(E_{\operatorname{quad}}(Q)<\epsilon_{\operatorname{comb}}\), and define \(\epsilon_Q:=\epsilon_{\operatorname{comb}}-E_{\operatorname{quad}}(Q)\). Specific quadrature rules are discussed in 4.6. We first prove an abstract fixed-\(2m\)-order overall complexity estimate, based on 1, for arbitrary \(Q\). For \(\delta>0\), define \[\label{abstract-rQ-delta-def} r_Q(\delta) := \left\lceil \max\left\{ \frac{T}{\Delta(R_Q,N,\delta)}, \left(2C_mT^{2m+1}\Phi_{m,Q}^{\ast}\right)^{1/2m}, \left( \frac{2e\,C_mT^{2m+1}\Lambda_{m,Q}}{\epsilon_Q} \right)^{1/2m} \right\} \right\rceil .\tag{23}\]

In what follows, define the \(m\)-dependent post-quadrature effective commutator scale \[\label{abstract-mu-mQ-def} \mu_{m,Q}(\delta) := \max\left\{ \frac{1}{\Delta(R_Q,N,\delta)}, (C_m\Phi_{m,Q}^{\ast})^{1/2m}, \left(C_m\Lambda_{m,Q}\right)^{1/2m} \right\}.\tag{24}\] It can then be checked that 23 implies that \[\label{abstract-rQ-mu-bound} r_Q(\delta) = \mathcal{O}\left( 1+ \mu_{m,Q}(\delta) \max\left\{ T,\, T^{1+\frac{1}{2m}},\, T^{1+\frac{1}{2m}}\epsilon_Q^{-\frac{1}{2m}} \right\} \right),\tag{25}\] for \(\epsilon_Q > 0\). Under the additional assumption \(0<\epsilon_Q<1\), this simplifies to \[\label{abstract-rQ-mu-bound-simple} r_Q(\delta) = \mathcal{O}\left( 1+ \mu_{m,Q}(\delta) \max\left\{ T,\, T^{1+\frac{1}{2m}}\epsilon_Q^{-\frac{1}{2m}} \right\} \right).\tag{26}\] In what follows, assume \(0<\epsilon_Q<1\). In 27, we allow mild MPF-order dependence in the effective commutator scale and optimize over the order to obtain complexity logarithmic in \(T/\epsilon_Q\), up to the growth envelope \(F\).

In the applications considered in this paper, we have \(F(m)=m^u\bigl(\log(e+m)\bigr)^v\) for some \(u,v\geq0\). Then 27 implies the following explicit polylogarithmic form: \[\begin{align} \label{abstract-optimized-block-complexity-polylog} \mathcal{Q}_{\operatorname{block}} = \mathcal{O}\Bigg( &\left( 1+ \mu_Q\max\{1,T\} \left(\log\left(e+\frac{T}{\epsilon_Q}\right)\right)^u \left( \log\left(e+\log\left(e+\frac{T}{\epsilon_Q}\right)\right) \right)^v \right) \notag\\ &\times \left(\log\left(e+\frac{T}{\epsilon_Q}\right)\right)^2 \left( \log\left(e+\log\left(e+\frac{T}{\epsilon_Q}\right)\right) \right)^2 \Bigg). \end{align}\tag{32}\] Indeed, the choice of \(m\) implies \(m = \mathcal{O}\left( \log\left(e+\frac{T}{\epsilon_Q}\right) \right)\). Therefore, we have \[F(m) = \mathcal{O}\left( \left(\log\left(e+\frac{T}{\epsilon_Q}\right)\right)^u \left( \log\left(e+\log\left(e+\frac{T}{\epsilon_Q}\right)\right) \right)^v \right).\] Substituting this into ?? gives 32 . The state complexity follows similarly.

4.4 General Commutator Scaling Specialization↩︎

The analysis in 4.3 allowed arbitrary quadrature under 1. We now specialize to 4, still with arbitrary quadrature. We first discuss the applicability of 4. For each \(G_{k_i}=H+k_iL\), we must compute \(\inf_{j\geq J}\alpha_{\operatorname{comm},j}(G_{k_i})^{-1/j}\). Define \(\mathcal{C}_j(H,L):=\max_{0\leq \ell\leq j}\mathcal{C}_{j,\ell}(H,L)\). We have \[\begin{align} \label{series-alpha-comm-bound} \alpha_{\operatorname{comm},j}(G_{k_i}) \leq \sum_{\ell=0}^{j}|k_i|^\ell\mathcal{C}_{j,\ell}(H,L) \leq \mathcal{C}_j(H,L)(j+1)\max\{1,|k_i|\}^j . \end{align}\tag{33}\] Let \(\chi_J(H,L):=\inf_{j\geq J}\mathcal{C}_j(H,L)^{-1/j}\). For \(J\geq 1\), we have \[\begin{align} \label{eq:comm-bound-root} \inf_{j\geq J}\alpha_{\operatorname{comm},j}(G_{k_i})^{-1/j} &\geq \frac{\inf_{j\geq J} ( (j+1)^{-1/j} \mathcal{C}_j(H,L)^{-1/j})}{\max\{1,|k_i|\}} \\ &\geq \frac{ (\inf_{j\geq J}(j+1)^{-1/j}) (\inf_{j\geq J}\mathcal{C}_j(H,L)^{-1/j})}{\max\{1,|k_i|\}} \\ &= \frac{\chi_J(H,L)}{(J+1)^{1/J}\max\{1,|k_i|\}} \geq \frac{\chi_J(H,L)}{(J+1)^{1/J}\max\{1,R_Q\}} := \kappa_Q. \end{align}\tag{34}\] The last inequality follows from \(|k_i|\leq R_Q\). Thus 4 applies uniformly to all \(G_{k_i}\) whenever \(\chi_J(H,L)>0\) and \(\Delta\leq\rho_Q\) for some \(0<\rho_Q<\kappa_Q\). The parameter \(\rho_Q\) is a uniform admissible step-size and enters the error profile because 4 gives an MPF error series in powers of \(\Delta^{j+\ell}\), while the abstract framework factors out \(\Delta^{2m+1}\). Set \(C_m:=\|\vec{a}\|_1\) and define \[\begin{align} \label{series-Dmnu-def} \Phi_{m,\nu}(H,L) := \sum_{\substack{j\in 2\mathbb{Z}_+\\ j\geq 2m}} \sum_{l=1}^{m} \frac{\rho_Q^{j+l-(2m+1)}}{l!} \sum_{\substack{j_1,\dots,j_l\in 2\mathbb{Z}_+\\ j_1+\cdots+j_l=j}} \sum_{\substack{\ell_1,\dots,\ell_l\geq 0\\ 0\leq \ell_\kappa\leq j_\kappa+1\\ \ell_1+\cdots+\ell_l=\nu}} \prod_{\kappa=1}^{l} \mathcal{C}_{j_\kappa+1,\ell_\kappa}(H,L). \end{align}\tag{35}\] for \(\nu\geq 0\). Define the discrete quadrature moment \(M_\alpha^Q:=\sum_{i\in\mathcal{I}_Q}|v_i|\,|k_i|^\alpha\). The relevant quantities appearing in the complexity estimates are \(\Lambda_{m,Q}:=\sum_{\nu\geq 0}\Phi_{m,\nu}(H,L)M_\nu^Q\) and \(\Phi^*_{m,Q} := \max_{i\in\mathcal{I}_Q} \Phi_{m,\nu}(H,k_iL)\). In 28, which is proved in 7.2, we record the convergence of the profiles before stating the complexity consequence. This is the only point where the strict inequality \(\rho_Q<\kappa_Q\) is used.

It remains to verify that the MPF error has the abstract form required in 1. By 4, for each \(i\in\mathcal{I}_Q\) and \(\Delta\leq\rho_Q\), we have \[\begin{align} \label{series-MPF-step-bound-original} \| U_{\operatorname{MPF},k_i}(\Delta) - e^{-iG_{k_i}\Delta}\| \leq& \|\vec{a}\|_1 \sum_{\substack{j\in 2\mathbb{Z}_+\\ j\geq 2m}} \sum_{l=1}^{m} \frac{\Delta^{j+l}}{l!} \sum_{\substack{j_1,\dots,j_l\in 2\mathbb{Z}_+\\ j_1+\cdots+j_l=j}} \prod_{\kappa=1}^{l} \left( \sum_{\ell=0}^{j_\kappa+1} |k_i|^\ell \mathcal{C}_{j_\kappa+1,\ell}(H,L) \right). \end{align}\tag{36}\] Since \(j\geq 2m\) and \(l\geq 1\), we have \(j+l\geq 2m+1\). Thus, for \(\Delta\leq \rho_Q\), we have \(\Delta^{j+l} \leq \Delta^{2m+1}\rho_Q^{j+l-(2m+1)}\). Let \(S_n(r):=\sum_{\ell=0}^{n}r^\ell\mathcal{C}_{n,\ell}(H,L)\). Hence, 36 gives \[\label{series-MPF-step-bound-abstract} \| U_{\operatorname{MPF},k_i}(\Delta) - e^{-iG_{k_i}\Delta} \| \leq C_m \Delta^{2m+1} \sum_{\substack{j\in 2\mathbb{Z}_+\\ j\geq 2m}} \sum_{l=1}^{m} \frac{\rho_Q^{j+l-(2m+1)}}{l!} \sum_{\substack{j_1,\dots,j_l\in 2\mathbb{Z}_+\\ j_1+\cdots+j_l=j}} \prod_{\kappa=1}^{l} S_{j_\kappa+1}(|k_i|) ,\tag{37}\] Denote the nested sum above as \(\Phi_m(|k_i|)\). Its maximum over \(i\in\mathcal{I}_Q\) is precisely \(\Phi_{m,Q}^\ast\). Summing against the quadrature weights gives \[\begin{align} \label{series-weighted-MPF-bound} \sum_{i\in\mathcal{I}_Q} |v_i| \| U_{\operatorname{MPF},k_i}(\Delta) - e^{-iG_{k_i}\Delta} \| &\leq C_m \Delta^{2m+1} \sum_{i\in\mathcal{I}_Q}|v_i|\Phi_m(|k_i|) = C_m \Delta^{2m+1} \Lambda_{m,Q}. \end{align}\tag{38}\]

Since \(R_m(\delta)=0\) in 4, \(r_Q\) and \(\mu_{m,Q}\) are defined as in 23 24 by setting \(\Delta(R_Q,N):=\rho_Q\) and suppressing \(\delta\). The fixed-order and optimized error-complexity bounds then follow from 26 27. We record in 29 an easier sufficient condition for computing commutator scaling in practice.

The upper bound \(\omega\) in 29 is easier to compute in practice. The condition \(\rho_Q\omega<1\) ensures both uniform applicability and convergence of the simplified error profiles. Moreover, \(\mu_{m,Q}\) in ?? is independent of \(m\) in most applications. Indeed, \(\alpha_Q\) is independent of \(m\), \(\omega\) is independent of \(m\) for applications considered in 5, and 30, proved in 7.3, gives \(\sup_{m\geq1} S_m<\infty\).

Since the MPF coefficients satisfy \(\|\vec{a}\|_1=\mathcal{O}(\log m)\), we also have \(\sup_{m\geq1}\|\vec{a}\|_1^{1/(2m)}<\infty\). This shows that \(\sup_{m \geq 1} (\|\vec{a}\|_1 S_m)^{1/(2m)} < \infty\). Thus, the sufficient condition \(\rho_Q\omega<1\) gives not only fixed-\(m\) convergence of the simplified profile, but also uniform-in-\(m\) boundedness required to invoke the optimized-order estimate.

4.5 Local Hamiltonian Specialization↩︎

We next use the locality-based MPF estimate from 5. Unlike 4, it includes the remainder \(R_m(\delta)=K_m\delta\), where \(K_m=\|\vec{a}\|_1\|\vec{b}\|_1\). Since 5 gives \(p_0=\mathcal{O}(\log(3N/\delta))\), set \(\delta_{p_0}:=3Ne^{-p_0}\). We choose \(p_0>\log(3N)\), so \(\delta_{p_0}\in(0,1)\). For \(p_0\in\mathbb{N}\), define \[\label{loc-mu-bar-def} \bar\mu_{m,Q}(p_0) := \max_{i\in\mathcal{I}_Q} \mu_{m,p_0}(G_{k_i}),\tag{39}\] where \(\mu_{m,p_0}\) is defined in ?? . Since \(H\) is \(q_H\)-local and \(g_H\)-extensive, and \(L\) is \(q_L\)-local and \(g_L\)-extensive, set \(q_Q := \max\{q_H,q_L\}\) and \(g_Q := g_H + R_Q g_L\). Then \(G_{k_i}=H+k_iL\) is uniformly \(q_Q\)-local and \(g_Q\)-extensive for all \(i\in\mathcal{I}_Q\). The corresponding step-size threshold is \[\label{loc-step-size} \Delta_{Q,p_0}:=\min\left\{\frac{1}{16e^3p_0q_Qg_Q},\frac{1}{4\bar\mu_{m,Q}(p_0)}\right\},\tag{40}\] Set \(C_m=2e^{1/2}\|\vec{a}\|_1\). The quantities appearing in the complexity estimate are \[\begin{align} \label{loc-Lambda-Phi-def} \Lambda_{m,Q}(p_0):=\sum_{i\in\mathcal{I}_Q} |v_i|\,\mu_{m,p_0}(G_{k_i})^{2m+1}, \quad \Phi^{\ast}_{m,Q}(p_0):= \bar\mu_{m,Q}(p_0)^{2m+1}. \end{align}\tag{41}\] For fixed \(m\ge1\) and \(p_0\in\mathbb{N}\), define \[\label{loc-r-choice} r_{Q,p_0}:=\left\lceil\max\left\{\frac{T}{\Delta_{Q,p_0}}, \left( 2C_mT^{2m+1} \Phi^{\ast}_{m,Q}(p_0) \right)^{1/2m}, \left( \frac{ 2e C_m T^{2m+1} \Lambda_{m,Q}(p_0) }{\epsilon_Q} \right)^{1/2m} \right\} \right\rceil\tag{42}\] In 31, proved in 7.4, we derive a growth condition on \(\mu_{m,p_0}\) needed to control \(r_{Q,p_0}\) before applying the results of 4.3.

By 31, \(\bar\mu_{m,Q}(p_0)\), and hence \(\Delta_{Q,p_0}^{-1}\), grow at most polynomially in \(p_0\). By 41 , we have \(\Phi_{m,Q}^{\ast}(p_0)\leq \bar\mu_{m,Q}(p_0)^{2m+1}\) and \(\Lambda_{m,Q}(p_0)\leq \alpha_Q\bar\mu_{m,Q}(p_0)^{2m+1}\). Thus, \(\Phi_{m,Q}^{\ast}(p_0)\), \(\Lambda_{m,Q}(p_0)\), and \(r_{Q,p_0}\) grow at most polynomially in \(p_0\). Since \(\delta_{p_0}=3Ne^{-p_0}\), we have \(r_{Q,p_0}K_m\delta_{p_0}\to 0\) as \(p_0\to\infty\). Hence, there exists \(p_0>\log(3N)\) such that \[\label{loc-p0-admissible} K_m\delta_{p_0} \le \min\left\{ \frac{1}{2r_{Q,p_0}}, \frac{\epsilon_Q}{2e\,\alpha_Qr_{Q,p_0}} \right\}.\tag{43}\] Under 5, \(R_m(\delta)=K_m\delta\). By 43 , we have \[r_mR_m(\delta_m)=r_mK_m\delta_m\leq \frac{1}{2}, \qquad r_m\alpha_QR_m(\delta_m)=r_m\alpha_QK_m\delta_m\leq \frac{\epsilon_Q}{2e}.\] Thus, the accumulated locality-remainder conditions in 26 hold. By 42 , \[\frac{T}{r_m}\leq \Delta_{Q,p_0}, \qquad \frac{C_mT^{2m+1}\Phi_{m,Q}^{\ast}(p_0)}{r_m^{2m}}\leq \frac{1}{2}, \qquad \frac{C_mT^{2m+1}\Lambda_{m,Q}(p_0)}{r_m^{2m}}\leq \frac{\epsilon_Q}{2e}.\] Therefore, we have \[r_m\eta_Q(T/r_m) = \frac{C_mT^{2m+1}\Phi_{m,Q}^{\ast}(p_0)}{r_m^{2m}} + r_mR_m(\delta_m) \leq 1.\] Combining the accumulated-remainder condition with the third inequality above gives the post-quadrature MPF error budget required in 26. Hence, all hypotheses of 26 hold, and the fixed-order error and complexity estimates follow. Define \[\label{loc-mu-mQ-def} \mu_{m,Q} := \max\left\{ \frac{1}{\Delta_{Q,p_0(m)}}, ( C_m\Phi_{m,Q}^{\ast}(p_0(m)) )^{1/2m}, \left( C_m \Lambda_{m,Q}(p_0(m)) \right)^{1/2m} \right\}.\tag{44}\] For the optimized statement, apply 27 with the admissible parameters \(\delta_m=\delta_{p_0(m)}\), \(r_Q(\delta_m)=r_m\), and \(\mu_{m,Q}(\delta_m)=\mu_{m,Q}\). The growth condition \(\mu_{m,Q}\leq \mu_QF(m)\) verifies the hypothesis in ?? for the chosen \(m\).

4.6 Quadrature Rules↩︎

We now fix a quadrature rule \(Q\). For a family \(\{Q_\eta\}\) of quadrature rules, we choose the quadrature parameters \(\eta\) so that \(E_{\operatorname{quad}}(Q_\eta)<\epsilon_{\operatorname{comb}}\), leaving the residual budget \(\epsilon_{Q_\eta}:=\epsilon_{\operatorname{comb}} - E_{\operatorname{quad}}(Q_\eta)>0\). The final choice of quadrature parameters is then governed by the complexity estimates in 26 27.

4.6.1 Uniform Trapezoidal Rule↩︎

We first specialize to the uniform trapezoidal rule. Following the Low–Somma discretization strategy [38], we compare the infinite LCHS integral with the infinite trapezoidal sum using strip analyticity, and then truncate the lattice to obtain the quadrature rule used in the LCU implementation. We use the notation \[\label{eq:FT-def} F_T(z) := \frac{1}{\sqrt{2\pi}}\hat{f}_{a,b}(z;c,d)U_z(T).\tag{45}\] Then \(W_\infty = \frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}} \hat{f}_{a,b}(k;c,d)U_k(T)\,dk = \int_{\mathbb{R}} F_T(k)\,dk\). We first state 32, which is proved in 7.5, a basic estimate used in the complexity analysis of the uniform trapezoidal quadrature rule.

We use the uniform trapezoidal rule because it converges exponentially for functions analytic in a strip around the real axis [38]. This is formalized below.

We now analyze the uniform trapezoidal rule in our setting. Fix \(R>0\) and \(h>0\), and set \(n_h:= \left\lceil\frac{R}{h}\right\rceil\) and \(R_h:=n_hh\). Then \(R_h\ge R\) and \(R_h<R+h\). Define the finite uniform trapezoidal rule \[\label{eq:uniform-trap-Q-def} Q_{R,h}^{\operatorname{trap}} := \{(k_i,w_i): k_i=qh,\;w_i=h,\;q=-n_h,\ldots,n_h\}.\tag{46}\] For this rule, \(R_{Q_{R,h}^{\operatorname{trap}}}=R_h\). The LCHS quadrature coefficients are \(v_q=\frac{h}{\sqrt{2\pi}}\hat{f}_{a,b}(qh;c,d)\). The corresponding ideal quadrature operator is \[\label{eq:uniform-lattice-rule} W_{R,h}^{\operatorname{trap},\operatorname{ideal}} := \frac{h}{\sqrt{2\pi}}\sum_{q=-n_h}^{n_h}\hat{f}_{a,b}(qh;c,d)U_{qh}(T) = h\sum_{q=-n_h}^{n_h}F_T(qh).\tag{47}\] The associated LCU normalization and discrete moments are \[\label{eq:trap-alpha-def} \alpha_{R,h}^{\operatorname{trap}} := \frac{h}{\sqrt{2\pi}}\sum_{q=-n_h}^{n_h}|\hat{f}_{a,b}(qh;c,d)|, \quad M_{\nu,R,h}^{\operatorname{trap}} := \frac{h}{\sqrt{2\pi}}\sum_{q=-n_h}^{n_h}|\hat{f}_{a,b}(qh;c,d)|\,|qh|^\nu,\tag{48}\] for \(\nu\ge0\). Moreover, the omitted lattice tail is \[\label{eq:lattice-tail} E_{\operatorname{tail},h}(R_h) := \frac{h}{\sqrt{2\pi}}\sum_{|q|>n_h}|\hat{f}_{a,b}(qh;c,d)|.\tag{49}\] 34 bounds the quadrature errors \(E_{\operatorname{tail},h}(R_h)\) and \(E_{\operatorname{trap}}(h,\rho)\), the latter arising from 33.

34 bounds the quadrature error by mesh and tail errors. Then 35 gives sufficient conditions for controlling both. In what follows, we use

\[\begin{align} C_{\operatorname{tail}} := \frac{c}{\sqrt{\pi}}\left(\frac{b+1}{b}\right)^{a-1}e^{d-\frac{1}{4c^2}}, \quad C_{\operatorname{ker}} := \frac{1}{2\pi}\left(\frac{b+1}{b}\right)^{a-1}e^{d-\frac{1}{4c^2}}, \quad \widetilde{C}_{\operatorname{ker}} := \frac{(b+1)^{a-1}}{2\pi}e^{d-\frac{1}{4c^2}}. \end{align}\] 35 bounds the quantities that determine the post-quadrature parameters for the trapezoidal quadrature rule.

12 shows that, when the kernel parameters satisfy \(a,d = \mathcal{O}(1)\), one may take \(c = \Omega(\sqrt{\log(1/\epsilon_{\operatorname{approx}})})\). Assuming that all error-tolerance parameters are of the same order as \(\epsilon > 0\), 35 implies \[|\mathcal{I}_{Q_{R,h}^{\operatorname{trap}}}| = \mathcal{O}\left(\frac{1}{\rho}\log^{2}\left(\frac{1}{\epsilon_q}\right)\right), \quad R_{Q_{R,h}^{\operatorname{trap}}} = \mathcal{O}\left(\log\left(\frac{1}{\epsilon_q}\right)\right).\] Moreover, 35 yields exponential mesh-error decay in \(1/h\), Gaussian tail decay in \(R\), and moment estimates that retain the rational decay of the kernel, improving the crude Gaussian scaling \(\mathcal{O}((\log(1/\epsilon_q))^{(\nu+1)/2})\) to the \(a\)-dependent scaling in ?? .

4.6.2 Sinh–Sinh Trapezoidal Rule↩︎

We next consider a non-uniform quadrature rule obtained by applying the trapezoidal rule after the change of variables \(k=\phi_\eta(x):=\eta\sinh x\), where \(0<\eta\le 2c\) is a free scale parameter. Set \(\lambda_\eta:=\eta/(2c)\). Under this change of variables, the Gaussian factor in \(\hat{f}_{a,b}\) becomes \(e^{-\phi_\eta(x)^2/(4c^2)}=e^{-\lambda_\eta^2\sinh^2 x}\), so the transformed integrand has double-exponential decay along the real axis. The choice \(\eta=2c\) gives the Gaussian-normalized case \(e^{-\lambda_\eta^2\sinh^2 x}=e^{-\sinh^2 x}\). See [74] for more on this quadrature rule. Define \[\label{eq:sinh-transformed-integrand} G_{T,\eta}(z) := F_T(\phi_\eta(z)) \phi_\eta'(z) = \eta\cosh z\, F_T(\eta\sinh z).\tag{55}\] Then, by the change of variables \(k=\phi_\eta(x)\), we have \(W_\infty=\int_{\mathbb{R}}F_T(k)\,dk=\int_{\mathbb{R}}G_{T,\eta}(x)\,dx\). The exponential convergence of the sinh–sinh rule follows by applying the infinite-line trapezoidal estimate to the transformed integrand \(G_{T,\eta}\). The relevant analyticity strip is now an \(x\)-strip rather than the original \(k\)-strip. Let \[\label{eq:sinh-beta-star-def} \beta_{\ast,\eta} := \sup\left\{\beta>0: G_{T,\eta} \text{ is analytic on an open set containing } \{z\in\mathbb{C}:|\operatorname{Im}z|\le\beta\}\right\}.\tag{56}\] Equivalently, \(\beta_{\ast,\eta}\) is limited by the nearest preimage of the poles \(k=-i\) and \(k=ib\) under \(z\mapsto \eta\sinh z\), and by the requirement that the transformed Gaussian decay on horizontal translates. In particular, one may restrict to \(\beta<\pi/4\), since \[\operatorname{Re}(\sinh^2(x+i\beta)) = \sinh^2x\cos(2\beta)-\sin^2\beta\] has positive quadratic growth in \(|x|\) whenever \(|\beta|<\pi/4\). For \(\beta\in(0,\min\{\beta_{\ast,\eta},\pi/4\})\), define \[\label{eq:sinh-denom-distances} \delta_{1,\eta}(\beta) := \inf_{\substack{x\in\mathbb{R}\\ |\omega|\le \beta}}\left|1-i\phi_\eta(x+i\omega)\right|, \qquad \delta_{b,\eta}(\beta) := \inf_{\substack{x\in\mathbb{R}\\ |\omega|\le \beta}}\left|b+i\phi_\eta(x+i\omega)\right|.\tag{57}\] Here \(\delta_{b,\eta}(\beta)\) is only needed when \(a>1\). Since \(\beta<\beta_{\ast,\eta}\), these quantities are strictly positive. Set \[\label{eq:sinh-rho-A-def} \rho_\beta:=\cos(2\beta)>0, \qquad A_{\eta,\beta} := \eta\sin\beta\,\bigl(|d|+T\|L\|\bigr).\tag{58}\] Define \[\label{eq:Msinh-explicit-def} M_{\operatorname{sinh},\eta}(\beta) := \frac{\eta (b+1)^{a-1}}{2\pi\,\delta_{1,\eta}(\beta)\,\delta_{b,\eta}(\beta)^{a-1}}\exp\left(d-\frac{1}{4c^2}+\lambda_\eta^2\sin^2\beta+A_{\eta,\beta}+\frac{A_{\eta,\beta}^2}{2\lambda_\eta^2\rho_\beta}\right)\sqrt{\frac{2\pi}{\lambda_\eta^2\rho_\beta}},\tag{59}\] with the convention \(\delta_{b,\eta}(\beta)^{a-1}=1\) when \(a=1\).

36 gives exponential mesh convergence in \(1/h\). The advantage of the sinh–sinh transformation is that the lattice truncation error decays double-exponentially in the \(x\)-truncation length. Fix \(Y>0\) and \(h>0\). As in 4.6.1, set \(n_h := \left\lceil\frac{Y}{h}\right\rceil\) and \(Y_h:=n_hh\). Define the finite sinh–sinh quadrature rule \[\label{eq:sinh-Q-def} Q_{Y,h}^{\operatorname{sinh},\eta} := \left\{(k_i,w_i): k_i=\eta\sinh(qh),\;w_i=\eta h\cosh(qh),\;q=-n_h,\ldots,n_h\right\}.\tag{60}\] For this rule, \(R_{Q_{Y,h}^{\operatorname{sinh},\eta}}=\eta\sinh(Y_h)\). The corresponding ideal quadrature operator is \[\label{eq:sinh-ideal-operator} W_{Y,h}^{\operatorname{sinh},\eta,\operatorname{ideal}} := h\sum_{q=-n_h}^{n_h}G_{T,\eta}(qh) = \frac{\eta h}{\sqrt{2\pi}}\sum_{q=-n_h}^{n_h}\cosh(qh)\,\hat{f}_{a,b}(\eta\sinh(qh);c,d)U_{\eta\sinh(qh)}(T).\tag{61}\] The LCHS quadrature coefficients are \(v_q=\frac{\eta h\cosh(qh)}{\sqrt{2\pi}}\hat{f}_{a,b}(\eta\sinh(qh);c,d)\). The associated LCU normalization is \(\alpha_{Y,h}^{\operatorname{sinh},\eta}:=\sum_{q=-n_h}^{n_h}|v_q|\), and the discrete moments are \(M_{\nu,Y,h}^{\operatorname{sinh},\eta}:=\sum_{q=-n_h}^{n_h}|v_q|\,|\eta\sinh(qh)|^\nu\), for \(\nu\ge0\). Moreover, the omitted transformed-lattice tail is \(E_{\operatorname{tail},h}^{\operatorname{sinh},\eta}(Y_h):=h\sum_{|q|>n_h}\|G_{T,\eta}(qh)\|\). 37 is the analogue of 34.

37 shows that the sinh–sinh rule has exponential mesh convergence in \(1/h\) and double-exponential truncation decay in the transformed cutoff \(Y\). Finally, 38 bounds the quantities determining the implementation complexity of the sinh quadrature rule. Since the proof is similar to that of 35, it is given in 7.6.

If one takes \(a,d=\mathcal{O}(1)\), \(c=\Theta(\sqrt{\log(1/\epsilon_{\operatorname{approx}})})\), and assumes that all error parameters are of order \(\epsilon>0\), then 38 shows that, for a fixed free-scale choice \(\eta=\Theta(\min(1,b))\), the transformed strip width satisfies \(\beta^{-1}=\mathcal{O}(1)\). Therefore \[|\mathcal{I}_{Q_{Y,h}^{\operatorname{sinh},\eta}}| = \mathcal{O}\left( \log\left(\frac{1}{\epsilon_q}\right) \log\log\left(\frac{1}{\epsilon_q}\right) \right),\] which is asymptotically smaller than \(\log^{3/2}(1/\epsilon_q)\).

5 Applications↩︎

In this section, we present applications of our algorithm and analyze its implementation complexity. 5.1 discusses the fractional diffusion equation with an imaginary potential, 5.2 discusses the advection–diffusion equation, and 5.3 the no-jump dynamics regime in a dissipative Ising model.

5.1 Fractional Diffusion Equation↩︎

Let \(0<s<1\). Consider the following partial differential equation: \[\label{eq:fractional-diffusion-pde} \begin{cases} \partial_t u(t,x) = -(-\Delta)^s u(t,x) + (a(t,x)-iv(t,x))u(t,x) + b(t,x), & (t,x)\in (0,T]\times\mathbb{T}^d,\\[0.5em] u(0,x)=u_0(x), & x\in\mathbb{T}^d. \end{cases}\tag{62}\]

Here \((-\Delta)^s\) is a fractional power of the positive semidefinite Laplacian on \(\mathbb{T}^d\), \(a,v:[0,T]\times\mathbb{T}^d\to\mathbb{R}\) are real-valued potentials, and \(b:[0,T]\times\mathbb{T}^d\to\mathbb{C}\) is a source term. If \(v\equiv0\), then 62 reduces to a linear fractional reaction–diffusion equation. The operator \((-\Delta)^s\) models non-local diffusion, \(a(t,x)\) gives a linear reaction term, and the imaginary potential \(-iv(t,x)\) produces local phase modulation, analogous to a Larmor-frequency-offset term in Bloch–Torrey-type models. Thus, 62 is a forced fractional Bloch–Torrey-type equation with an additional real potential. Fractional Bloch–Torrey models have been used for anomalous diffusion in complex media, especially in magnetic resonance applications [40]–[44]. In what follows, we take \(a,b\equiv0\) and assume that \(v\) is time-independent.

Discretize 62 on an \(M^d\) periodic grid with mesh size \(h=1/M\) and \(N=M^d\). Let \(F_M\) be the one-dimensional discrete Fourier transform, \(F_N:=F_M^{\otimes d}\), and let \(-\Delta_{h,1}\) be the one-dimensional periodic discrete Laplacian, with eigenvalues \(4M^2\sin^2(\pi m/M)\), \(m\in\mathbb{Z}_M\). The \(d\)-dimensional periodic discrete Laplacian is given by the Kronecker sum \[\label{eq:periodic-laplacian-kronecker-sum} -\Delta_{h,d} = \sum_{j=1}^d I^{\otimes(j-1)} \otimes (-\Delta_{h,1}) \otimes I^{\otimes(d-j)}, \quad -\Delta_{h,1} = \frac{1}{h^2} \begin{pmatrix} 2 & -1 & & & & -1\\ -1 & 2 & -1 & & & \\ & -1 & 2 & \ddots & & \\ & & \ddots& \ddots & -1 & \\ & & & -1 & 2 & -1\\ -1 & & & & -1 & 2 \end{pmatrix}.\tag{63}\] Consequently, \(-\Delta_{h,d}\) is diagonalized by \(F_N=F_M^{\otimes d}\), with eigenvalues \(\lambda_m = 4M^2 \sum_{\alpha=1}^d \sin^2\left(\frac{\pi m_\alpha}{M}\right)\) for \(m=(m_1,\dots,m_d)\in \mathbb{Z}_M^d\). Therefore, the discretized fractional Laplacian is given by \[\label{eq:discrete-fractional-laplacian} L_s := (-\Delta_{h,d})^s = F_N^\dagger \operatorname{diag}\left[ \left( 4M^2 \sum_{\alpha=1}^d \sin^2\left(\frac{\pi m_\alpha}{M}\right) \right)^s \right]_{m\in\{0,\dots,M-1\}^d} F_N .\tag{64}\] Discretize the potential as \(H:=\operatorname{diag}(v(x_i/M))_{i=1}^N\), with \(x_i\in\mathbb{Z}_M^d\). If \(\vec{u}(t)\) approximates \(u(t,x)\) on the grid and \(\vec{u}_0\) is the initial condition, then the semi-discrete equation is \[\label{eq:fractional-diffusion-discrete} \begin{cases} \dfrac{d\vec{u}}{dt} = -(L_s+iH)\vec{u}(t), \\ \vec{u}(0)=\vec{u}_0 . \end{cases}\tag{65}\]

The estimate in 5 does not apply to 65 , since \(L_s\) is non-local for \(0<s<1\). On \(\mathbb{R}^d\), the fractional Laplacian has the singular integral representation \[\label{eq:frac-lap-singular-integral} (-\Delta)^s u(x) = \frac{2^{2s}\Gamma (s+\frac{d}{2})}{\pi^{d/2}|\Gamma (-s)|} \operatorname{p.v.} \int_{\mathbb{R}^d} \frac{u(x)-u(y)}{|x-y|^{d+2s}}\,dy \qquad 0<s<1.\tag{66}\] On \(\mathbb{T}^d\), the kernel is obtained by periodization. Thus, \((-\Delta)^s u(x)\) depends on all spatial values \(u(y)\), with algebraically decaying interactions. This non-locality is inherited by \(L_s\), so the standard qubit encoding of \(L_s\) is non-local. On the other hand, 4 does apply. Using \(\|[A,B]\|\leq 2\|A\|\|B\|\), we obtain \[\label{eq:fractional-crude-Cjl-bound} C_{j,\ell}(H,L_s) \leq \binom{j}{\ell} 2^{j-1} \|H\|^{j-\ell} \|L_s\|^\ell .\tag{67}\] With \(C_j(H,L_s):=\max_{0\leq \ell\leq j} C_{j,\ell}(H,L_s)\), we have \(C_j(H,L_s) \leq 2^{j-1}\bigl(\|H\|+\|L_s\|\bigr)^j\). It follows that \[\chi_J(H,L_s) := \inf_{j\geq J} C_j(H,L_s)^{-1/j} \geq \frac{1}{2(\|H\|+\|L_s\|)} > 0 .\]

Since \(L_s\succeq0\) and \(H=H^\dagger\), our algorithm applies. We compute the commutator scaling used in the complexity analysis. Since \(L_s\) is diagonalized by \(F_N\), we have \((F_NL_sF_N^\dagger)_{p,q}=\lambda_p^s\delta_{p,q}\), where \(\lambda_p\) is an eigenvalue of \(-\Delta_{h,d}\). The Fourier-basis matrix of \(H\) is \[\label{eq:fractional-FHF} (F_NHF_N^\dagger)_{p,q} = \frac{1}{N}\sum_{x\in\mathbb{Z}_M^d}v(x/M)e^{-2\pi i(p-q)\cdot x/M} =: \widehat V_{p-q}.\tag{68}\] 40, which is proved in 8.1, estimates differences of fractional-Laplacian eigenvalues. Define \(\beta_s:=\min\{2s,1\}\), \(\theta_s:=2s-\beta_s=\max\{0,2s-1\}\), and let \(|\cdot|_{\operatorname{per}}\) denote distance on \(\mathbb{Z}_M^d\).

Using 40, 41 estimates the scaling of \(\operatorname{ad}_{L_s}^j(H)\) for \(j\geq1\). This serves as a prototype for bounding general commutator words in 42.

We now extend 41 to arbitrary nested commutator words in 42, which is proved in 8.2. In 42, we use the weaker consequence of 40 \[\label{eq:fractional-symbol-difference-weak} |\lambda_p^s-\lambda_q^s| \leq C_{s,d}N^{\theta_s/d}\bigl(1+|p-q|_{\operatorname{per}}\bigr)^{\beta_s}.\tag{72}\] The constant in 42 is not sharp since the proof uses the crude bound \(\|\operatorname{ad}_H(A)\|_{a,\mathcal{F}}\leq2\mathcal{V}_a^{(N)}\|A\|_{a,\mathcal{F}}\) in 131 . For special words, direct Fourier-basis formulas give sharper constants. For example, 41 uses \((F_N\operatorname{ad}_{L_s}^{\ell}(H)F_N^\dagger)_{p,q}=(\lambda_p^s-\lambda_q^s)^\ell\widehat V_{p-q}\) together with fractional-Laplacian eigenvalue differences. More general nested commutators finite differences of \(p\mapsto\lambda_p^s\), which may yield sharper constants. We do not pursue these refinements here.

We now analyze the implications of 42 for the complexity of our algorithm. Since \(\mathcal{C}_{j,\ell}(H,L_s)\) collects all commutator words with exactly \(\ell\) occurrences of \(L_s\), 42 implies \[\label{eq:fractional-Cjell-application-bound} \mathcal{C}_{j,\ell}(H,L_s) \leq 2^j 2^{j-\ell} C_{s,d}^{\ell} (\mathcal{V}_{j\beta_s}^{(N)})^{j-\ell} N^{\ell\theta_s/d}\tag{73}\] for \(1\leq \ell\leq j-1\). Define \[\label{eq:fractional-Aj-def} A_{V,j}^{(N)} := \max_{1\leq \ell\leq j-1} \left( 2^j 2^{j-\ell} C_{s,d}^{\ell} (\mathcal{V}_{j\beta_s}^{(N)})^{j-\ell} \right)^{1/j}.\tag{74}\] 73 gives \(\mathcal{C}_{j,\ell}(H,L_s) \leq (A_{V,j}^{(N)})^j N^{\ell\theta_s/d}\). Using the definition of \(\alpha_{Q,j}\) from 29, we obtain \[\begin{align} \alpha_{Q,j}\le \sum_{\ell=0}^{j} R_Q^\ell \mathcal{C}_{j,\ell}(H,L_s)\le (A_{V,j}^{(N)})^j\sum_{\ell=1}^{j-1}(R_QN^{\theta_s/d})^\ell\le \bigl(2A_{V,j}^{(N)}(1+R_QN^{\theta_s/d})\bigr)^j, \end{align}\] where the last inequality follows from \(\sum_{\ell=1}^{j-1}x^\ell\leq (2(1+x))^j\) for \(x\geq0\). Hence, if \(A_V^{(N)} := \sup_{j\ge 2} A_{V,j}^{(N)} < \infty\), we have \(\alpha_{Q,j} \leq ( 2A_V^{(N)} (1+R_QN^{\theta_s/d} ) )^j\) Therefore, the root-growth scale \(\omega\) from 29 satisfies \[\label{eq:fractional-omega-bound} \omega = \mathcal{O} ( A_V^{(N)} (1+R_QN^{\theta_s/d}) ).\tag{75}\] Choose \(\rho_Q=c/\omega\) for some fixed \(0<c<1\). Then \(\rho_Q\omega=c<1\), so 29 applies, and \(\rho_Q^{-1}=\mathcal{O}(\omega)\). Combining 75 with ?? gives \[\label{eq:fractional-muQ-scale} \mu_Q = \mathcal{O} ( A_V^{(N)} (1+R_QN^{\theta_s/d})),\tag{76}\] where the \(\mathcal{O}\)-notation absorbs the fixed finite factor \(S\) from 29. More explicitly, the optimized MPF overhead contributes the polylogarithmic factor \[\begin{align} \label{eq:fractional-polylog-factor} \operatorname{polylog}_Q(T,\epsilon_Q,N) &:= \mathcal{O} \left( \log\left(\frac{ A_V^{(N)} (1+R_QN^{\theta_s/d}) T}{\epsilon_Q}\right)^2 \left( \log\log\left(\frac{ A_V^{(N)} (1+R_QN^{\theta_s/d}) T}{\epsilon_Q}\right) \right)^2\right), \\ &= \mathcal{O} \left( \log\left(\frac{ A_V^{(N)} (1+R_QN^{\theta_s/d}) T}{\epsilon_Q}\right)^4\right), \end{align}\tag{77}\] where we have assumed that \(A_V^{(N)} (1+R_QN^{\theta_s/d}) T/\epsilon_Q \geq e^e\). Substituting 76 into ?? with \(F(m)\equiv1\), we obtain \[\mathcal{Q}_{\operatorname{block}} = \mathcal{O} ( A_V^{(N)} (1+R_QN^{\theta_s/d} )T \operatorname{polylog}_Q(T,\epsilon_Q,N) ).\] Similarly, provided \(W_Q^{\operatorname{MPF}}u_0\neq0\), substituting 76 into ?? , we obtain \[\label{eq:fractional-state-complexity-final} \mathcal{Q}_{\operatorname{state}}(u_0) = \mathcal{O} ( \chi_Q(u_0) A_V^{(N)} (1+R_QN^{\theta_s/d} )T \operatorname{polylog}_Q(T,\epsilon_Q,N) ).\tag{78}\]

5.2 Advection-Diffusion Equation↩︎

Let \(a>0\) and \(b\in\mathbb{R}\). Consider the following partial differential equation: \[\label{eq:generalized-heat-pde} \begin{cases} \partial_t u(t,x) = a\Delta u(t,x) + b \nabla\cdot(\vec{1} u(t,x)) - c(t,x)u(t,x) + f(t,x), & (t,x)\in (0,T]\times (0,1)^d,\\ u(0,x)=u_0(x), & x\in [0,1]^d,\\ u(t,x)=0, & (t,x)\in[0,T]\times\partial [0,1]^d. \end{cases}\tag{79}\] 79 is a \(d\)-dimensional advection–diffusion equation with homogeneous Dirichlet boundary conditions, where \(a\Delta u\) is diffusion, \(b\nabla\cdot(\vec{1}u)\) is drift, \(c(t,x)\geq0\) is absorption, and \(f(t,x)\) is the source. Such equations model heat and mass transfer, reactive systems, drug delivery, and related transport phenomena [45]–[48]. We take \(c,f\equiv0\) for simplicity. Discretize using \(M+1\) equidistant grid points per direction, mesh size \(h=1/M\), and \(N=(M+1)^d\). The \(d\)-dimensional Dirichlet finite-difference Laplacian is the Kronecker sum \[-\Delta_{h,d} := \sum_{j=1}^d I^{\otimes(j-1)} \otimes (-\Delta_{h,1}) \otimes I^{\otimes(d-j)}, \quad -\Delta_{h,1} := \frac{1}{h^2} \begin{pmatrix} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & -1 & 2 & \ddots & \\ & & \ddots & \ddots & -1 \\ & & & -1 & 2 \end{pmatrix}.\] By abuse of notation, we use the same notation for the Dirichlet finite-difference Laplacian and the periodic finite-difference Laplacian from 5.1. Similarly, the \(d\)-dimensional Dirichlet finite-difference drift matrix is the Kronecker sum \[D_{h,d} := \sum_{j=1}^d I^{\otimes(j-1)} \otimes D_{h,1} \otimes I^{\otimes(d-j)}, \quad D_{h,1} = \frac{1}{2h} \begin{pmatrix} 0 & 1 & & & \\ -1 & 0 & 1 & & \\ & -1 & 0 & \ddots & \\ & & \ddots & \ddots & 1 \\ & & & -1 & 0 \end{pmatrix}.\] Here \(\Delta_{h,d}\) discretizes \(\Delta=\sum_{i=1}^d\partial_{x_i}^2\), while \(D_{h,d}\) discretizes \(\nabla\cdot(\vec{1}u)=\sum_{i=1}^d\partial_{x_i}u\). Note that Since \(\Delta_{h,d}=\Delta_{h,d}^\dagger\preceq0\) \(D_{h,d}^\dagger=-D_{h,d}\). Define \(L := -a\Delta_{h,d}\) and \(H := ibD_{h,d}\). If \(\vec{u}(t)\) approximates \(u(t,x)\) on the grid and \(\vec{u}_0\) is the initial condition, then the semi-discrete equation is in dissipative-plus-Hamiltonian form is \[\label{eq:generalized-heat-lchs-form} \begin{cases} \dfrac{d\vec{u}}{dt} = -\left(L+iH\right) \vec{u}(t),\\ \vec{u}(0)=\vec{u}_0 . \end{cases}\tag{80}\]

We now derive commutator bounds. For \(d=1\), define \(A_1:=-h^2\Delta_{h,1}\) and \(B_1:=2hD_{h,1}\).

44 shows that the non-commutativity of \(A_1\) and \(B_1\) is localized at the boundary. Consequently, repeated commutators remain boundary-supported. This is made precise in 45, which is proved in 8.3.

We lift the \(d=1\) estimate to the \(d\)-dimensional tensor-product discretization in 46.

Since \(N=(M+1)^d\) and \(h^{-1}=M\leq N^{1/d}\) 46 implies \[\label{eq:dirichlet-Cjell-N-scaling} \mathcal{C}_{j,\ell}(H,L) \leq d 2^{2j+\ell-2} a^\ell |b|^{j-\ell} N^{(j+\ell)/d}.\tag{83}\] We now translate 46 into a complexity bound. Let \(\alpha_{Q,j} \leq \sum_{\ell=1}^{j-1} R_Q^\ell \mathcal{C}_{j,\ell}(H,L)\) be as in 29. Using ?? , for every \(j \geq 2\) we have

\[\begin{align} \label{eq:advection-diffusion-alphaQj-bound} \alpha_{Q,j} &\leq d2^{2j-2} \sum_{\ell=1}^{j-1} (2aR_Qh^{-2} )^\ell (|b|h^{-1} )^{j-\ell} \leq d2^{2j-2} ( |b|h^{-1} + 2aR_Qh^{-2} )^j . \end{align}\tag{84}\] Since \(h^{-1}=M\leq N^{1/d}\), the root-growth scale \(\omega\) from 29 satisfies \[\label{eq:advection-diffusion-omega-N-scale} \omega = \mathcal{O} ( |b|N^{1/d} + aR_QN^{2/d}).\tag{85}\] Choose \(\rho_Q=c/\omega\) for fixed \(0<c<1\). Then \(\rho_Q\omega=c<1\), so 29 applies and \(\rho_Q^{-1}=\mathcal{O}(\omega)\). Combining 85 with ?? gives \[\label{eq:advection-diffusion-muQ-scale} \mu_Q= \mathcal{O} ( |b|N^{1/d} + aR_QN^{2/d}),\tag{86}\] where the \(\mathcal{O}\)-notation absorbs the fixed finite factor \(S\) from 29. For the optimized-order estimate, the MPF overhead contributes the polylogarithmic factor \[\begin{align} \label{eq:advection-diffusion-polylog-factor} \operatorname{polylog}_Q(T,\epsilon_Q,N) &:= \mathcal{O}\left( \log\left( \frac{ \left(|b|N^{1/d}+aR_QN^{2/d}\right)T }{\epsilon_Q} \right)^2 \left( \log\log\left( \frac{ \left(|b|N^{1/d}+aR_QN^{2/d}\right)T }{\epsilon_Q} \right) \right)^2 \right) \nonumber\\ &= \mathcal{O}\left( \log\left( \frac{ \left(|b|N^{1/d}+aR_QN^{2/d}\right)T }{\epsilon_Q} \right)^4 \right), \end{align}\tag{87}\] where in the second line we assume \((|b|N^{1/d}+aR_QN^{2/d})T/\epsilon_Q \geq e^e\). Substituting 86 into ?? , we obtain \[\label{eq:advection-diffusion-block-complexity-final} \mathcal{Q}_{\operatorname{block}} = \mathcal{O} ( ( |b|N^{1/d} + aR_QN^{2/d} ) T \operatorname{polylog}_Q(T,\epsilon_Q,N) ).\tag{88}\] Similarly, provided \(W_Q^{\operatorname{MPF}}\vec{u}_0\neq0\), substituting 86 into ?? gives \[\label{eq:advection-diffusion-state-complexity-final} \mathcal{Q}_{\operatorname{state}}(\vec{u}_0) = \mathcal{O} ( \chi_Q(\vec{u}_0) ( |b|N^{1/d} + aR_QN^{2/d} ) T \operatorname{polylog}_Q(T,\epsilon_Q,N) ).\tag{89}\]

5.3 No-Jump Dynamics for a Dissipative Ising Model↩︎

We next consider the no-jump regime of a dissipative transverse-field Ising model with local spontaneous emission, formulated as a Lindblad dynamics. The no-jump regime is obtained from the standard quantum-trajectory formalism [49]–[52]. Let \(\mathcal{G}=(V,E(V))\) be a finite graph. The closed-system transverse-field Ising Hamiltonian is \[\label{eq:tfim-hamiltonian} H = \sum_{j\in V} h_j X_j + \sum_{\langle j,\ell\rangle\in E(V)} J_{j\ell} Z_j Z_\ell,\tag{90}\] where \(h_j \in \mathbb{R}\) are transverse-field strengths and \(J_{j\ell} \in \mathbb{R}\) are Ising couplings. Local spontaneous emission is modeled by Lindblad jump operators \(L_j=\sqrt{2\gamma_j}\,\sigma_j^-\), where \(\gamma_j\geq0\) and \(\sigma_j^-=\ket{0}\!\bra{1}\). Thus \(L_j\) represents local decay at site \(j\) with rate \(\gamma_j\). The associated Lindblad equation for \(\rho\) is \[\label{eq:tfim-lindblad} \frac{d\rho}{dt} = -i[H,\rho] + \sum_{j\in V} ( L_j\rho L_j^\dagger - \frac{1}{2}\{L_j^\dagger L_j,\rho\} ).\tag{91}\]

91 gives the unconditional evolution averaged over all emission histories. In the quantum-trajectory picture, \(L_j\rho L_j^\dagger\) represents a detected jump at site \(j\). Conditioning on records with no detected emissions gives the no-jump evolution, generated by the effective Hamiltonian \[\label{eq:tfim-effective-Hamiltonian} H_{\operatorname{eff}} = H-\frac{i}{2}\sum_{j\in V}L_j^\dagger L_j = \sum_{j\in V} h_j X_j + \sum_{\langle j,\ell\rangle\in E(V)} J_{j\ell} Z_j Z_\ell - i\sum_{j\in V}\gamma_j n_j,\tag{92}\] where \(n_j=\sigma_j^+\sigma_j^-= |1\rangle \langle 1|=(I-Z_j)/2\). The final term reflects the no-jump conditioning, where configurations more likely to emit are suppressed when no photon is detected. The corresponding dynamics are given by \[\label{eq:tfim-no-jump-schrodinger-LCHS-derivation} i\frac{d}{dt}\ket{\psi(t)} = H_{\mathrm{eff}}\ket{\psi(t)} .\tag{93}\] Let \(L=\sum_{j\in V}\gamma_j n_j = \sum_{j\in V}\gamma_j \frac{I-Z_j}{2}\). Then \(L\succeq0\) and \(H_{\operatorname{eff}} = H-iL\). Multiplying 93 by \(-i\), we obtain \[\label{eq:tfim-no-jump-LCHS-form} \frac{d}{dt}\ket{\psi(t)} = -i(H-iL)\ket{\psi(t)} = -(L+iH)\ket{\psi(t)}.\tag{94}\]

In what follows, take \(h_j\equiv h\), \(J_{j\ell}\equiv J\), and \(\gamma_j\equiv\gamma\), and assume \(\mathcal{G}\) has maximum degree \(D\), independent of \(|V|\). We first verify locality and extensivity. The Hamiltonian \(H\) is \(2\)-local:. Indeed, \(hX_j\) is supported on on one site and \(JZ_jZ_\ell\) is supported on two sites. For each site \(j\), at most one onsite term and \(D\) interaction terms contain \(j\). Since \(\|X_j\|=\|Z_jZ_\ell\|=1\), \[\sum_{\alpha:\,j\in\operatorname{supp}(H_\alpha)} \|H_\alpha\| \leq |h|+D|J|. \label{eq:tfim-H-extensive-bound}\tag{95}\] Hence, \(H\) is \(g_H\)-extensive with \(g_H:=|h|+D|J|\). Similarly, \(L\) is \(1\)-local because each term \(\gamma n_j\) is supported only on site \(j\). Hence, \(q_L=1\). Since \(n_j\) is a projection, \(\|n_j\|=1\), and therefore \[\sum_{\alpha:\,j\in\operatorname{supp}(L_\alpha)} \|L_\alpha\| = \|\gamma n_j\| = \gamma . \label{eq:tfim-L-extensive-bound}\tag{96}\] Thus \(L\) is \(g_L\)-extensive with \(g_L:=\gamma\). Hence, 4.5 implies that \(\{G_{k_i}\}_{i\in\mathcal{I}_Q}\) is \(q_Q\)-local and \(g_Q\)-extensive, with \(q_Q=\max\{q_H,q_L\}=2\) and \(g_Q\leq |h|+D|J|+R_Q\gamma\).

We now derive commutator bounds. Our approach uses the Pauli incompatibility graph, whose vertices are local Pauli summands and whose edges connect anti-commuting Pauli strings. Such commutativity-sensitive counting is in graph-based ordering methods in Trotter-error analyses [75]. Since \(n_j=(I-Z_j)/2\), the identity part commutes with all terms and does not affect nested commutators. Thus, for commutator estimates, \(G_k=H+kL\) may be replaced by \[\overline{G}_k = h\sum_{j\in V}X_j + J\sum_{\langle j,\ell\rangle\in E(V)}Z_jZ_\ell - \frac{k\gamma}{2}\sum_{j\in V}Z_j . \label{eq:tfim-Gbar-k-def}\tag{97}\] This replacement leaves every nested commutator of length at least two unchanged. In what follows, we will use the following standard computation for commutators of Pauli strings \(P\) and \(Q\) \[\|[xP,yQ]\| = \begin{cases} 0, & PQ=QP,\\ 2|xy|, & PQ=-QP, \end{cases} \label{eq:pauli-commutator-norm}\tag{98}\] for \(x,y \in \mathbb{R}\).

47 gives estimates for each commutator length \(p\). To apply the estimate in 4.5, we need a finite-\(p_0\) growth scale controlling all \(2\leq p\leq p_0\). Define \[\label{eq:tfim-ThetaQp0-def} \Theta_{Q,p_0} := \max_{2\leq p\leq p_0} \left( 4^{p-1}(p-1)! \,2|h| \left(D|J|+\frac{R_Q\gamma}{2}\right) \left( |h|+D|J|+\frac{R_Q\gamma}{2} \right)^{p-2} \right)^{1/p}.\tag{101}\] 47 gives \(\alpha_{\operatorname{comm},p}(G_{k_i})\leq |V|\Theta_{Q,p_0}^p\) for \(2\leq p\leq p_0\). 48 converts this finite-order commutator estimate into a bound on \(\overline{\mu}_{m,Q}(p_0)\) for the locality-based complexity analysis.

We now instantiate the locality-based complexity estimate. By 48 and the definitions of \(\Phi_{m,Q}^{\ast}(p_0)\) and \(\Lambda_{m,Q}(p_0)\) in 41 , we have \[\label{eq:tfim-Phi-Lambda-bound} \Phi_{m,Q}^{\ast}(p_0)\leq (|V|^{1/3}(1+p_0)^{1/3}\Theta_{Q,p_0} )^{2m+1}, \quad \Lambda_{m,Q}(p_0)\leq \alpha_Q (|V|^{1/3}(1+p_0)^{1/3}\Theta_{Q,p_0} )^{2m+1}.\tag{102}\] Together with the locality and extensivity bounds, \(q_Q=2\) and \(g_Q\leq |h|+D|J|+R_Q\gamma\), so 4.5 applies. Choose \(p_0\), \(r_m\), and \(\delta_m\) as in 4.5. By 101 , we have \[\label{eq:tfim-Theta-polynomial-bound} \Theta_{Q,p_0}\leq Cp_0\left(|h|+D|J|+\frac{R_Q\gamma}{2}\right)\tag{103}\] for some \(C>0\). Hence, we have \[\label{eq:tfim-mubar-polynomial-bound} \overline{\mu}_{m,Q}(p_0)\leq C|V|^{1/3}\left(|h|+D|J|+\frac{R_Q\gamma}{2}\right)(1+p_0)^{4/3}.\tag{104}\] Since \(\delta_{p_0}=3|V|e^{-p_0}\) and \(K_m=\mathcal{O}(m^2(\log m)^2)\), the admissibility condition in 43 holds with \(p_0(m)=\mathcal{O}(\log(e+m))\). Thus, the commutator-profile contribution grows as \((\log(e+m))^{4/3}\). For the optimized-order condition, the relevant quantity is \(\mu_{m,Q}\), not only \(\overline{\mu}_{m,Q}(p_0)\).The rooted \(\Phi_{m,Q}^{\ast}\)- and \(\Lambda_{m,Q}\) terms in 44 contribute powers of \(\overline{\mu}_{m,Q}(p_0)^{1+1/(2m)}\). Hence a safe growth function is \(F(m):=(\log(e+m))^2\). More explicitly, one may take \[\label{eq:tfim-muQ-choice} \mu_Q:=C_Q\max\left\{1,\;q_Qg_Q,\;|V|^{1/3}\left(|h|+D|J|+\frac{R_Q\gamma}{2}\right),\;\alpha_Q^{1/2}\left[|V|^{1/3}\left(|h|+D|J|+\frac{R_Q\gamma}{2}\right)\right]^{3/2}\right\},\tag{105}\] where \(C_Q>0\) is independent of \(m\), \(T\), and \(\epsilon_Q\), and absorbs the implicit constants. With this choice, the optimized-order growth condition analogous to ?? holds as \(\mu_{m,Q}\leq \mu_Q(\log(e+m))^2\). Choose \(m=\max\left\{m_0,\left\lceil\log\left(e+\frac{T}{\epsilon_Q}\right)\right\rceil\right\}\). Applying the optimized complexity estimate with \(F(m)=(\log(e+m))^2\) gives the block encoding complexity to be \[\label{eq:tfim-optimized-block-complexity} \mathcal{Q}_{\operatorname{block}}=\mathcal{O}\left(\left(1+\mu_Q\max\{1,T\}(\log(e+m))^2\right)\left(\log\left(e+\frac{T}{\epsilon_Q}\right)\right)^2\left(\log\log\left(e^e+\frac{T}{\epsilon_Q}\right)\right)^2\right).\tag{106}\] Similarly, the normalized output state complexity is \[\label{eq:tfim-optimized-state-complexity} \mathcal{O}\left(\chi_Q(u(0))\left(1+\mu_Q\max\{1,T\}(\log(e+m))^2\right)\left(\log\left(e+\frac{T}{\epsilon_Q}\right)\right)^2\left(\log\log\left(e^e+\frac{T}{\epsilon_Q}\right)\right)^2\right).\tag{107}\]

6 Omitted Proofs in 3↩︎

6.1 Proof of 8↩︎

6.2 Proof of 9↩︎

6.3 Proof of 10↩︎

6.4 Proof of 12↩︎

6.5 Proof of 13↩︎

7 Omitted Proofs in 4↩︎

7.1 Proof of 19↩︎

7.2 Proof of 28↩︎

7.3 Proof of 30↩︎

7.4 Proof of 31↩︎

7.5 Proof of 32↩︎

7.6 Proof of 38↩︎

8 Omitted Proofs in 5↩︎

8.1 Proof of 40↩︎

8.2 Proof of 42↩︎

8.3 Proof of 45↩︎

9 Extension to Non-Homogeneous and Time-Dependent Cases↩︎

In 3 and 4, we analyzed the time-independent homogeneous case in detail. In this section, we explain how the same framework extends to more general cases, namely the non-homogeneous case in 9.1 and the time-dependent case in 9.2.

9.1 Non-Homogeneous Case↩︎

We first consider the time-independent non-homogeneous problem \[\label{eq:non-homogeneous-time-independent} \frac{d}{dt}u(t)=-Au(t)+b(t).\tag{134}\] As before, write \(A=L+iH\), with \(L,H\) Hermitian and \(L\succeq 0\). By Duhamel’s principle, the solution of 134 is \[\label{eq:duhamel-time-independent} u(T)=e^{-AT}u(0)+\int_0^T e^{-A(T-s)}b(s)\,ds.\tag{135}\] 135 shows that the non-homogeneous term is a superposition of homogeneous evolutions over times \(T-s\). Applying the LCHS representation to each homogeneous propagator in the Duhamel term gives the two-variable representation \[\label{eq:duhamel-LCHS-double-integral} \int_0^T e^{-A(T-s)}b(s)\,ds \approx \frac{1}{\sqrt{2\pi}}\int_0^T\int_{\mathbb{R}}\hat{f}_{a,b}(k;c,d)U_k(T-s)b(s)\,dk\,ds.\tag{136}\] Thus, the non-homogeneous contribution involves two discretizations: the LCHS quadrature \(Q\) in the Fourier variable \(k\), and the time quadrature \(P\) in the Duhamel variable \(s\). This viewpoint is used below. For \(\tau\in[0,T]\), define the ideal and MPF-implemented post-quadrature operators by \[\label{eq:WQ-tau-ideal-mpf} W_Q^{\operatorname{ideal}}(\tau):=\sum_{i\in\mathcal{I}_Q}v_iU_{k_i}(\tau), \qquad W_Q^{\operatorname{MPF}}(\tau):=\sum_{i\in\mathcal{I}_Q}v_iU^{(2m)}_{k_i}(\tau),\tag{137}\] The quantities \(\alpha_Q\), \(\Lambda_{m,Q}\), and \(\Phi_{m,Q}^{\ast}\) are defined as before. Let \(P=\{(s_p,\omega_p):p\in\mathcal{I}_P\}\) be a quadrature rule on \([0,T]\) for the integral in 136 . We write \[\tau_p:=T-s_p, \qquad B_P:=\sum_{p\in\mathcal{I}_P}|\omega_p|\,\|b(s_p)\|.\]

We assume \(B_P<\infty\). The assumption \(b\in L^1([0,T])\) suffices for 135 , but does not define point values \(b(s_p)\). The quadrature implementation is therefore stated under the additional source-access assumption that the selected values \(b(s_p)\) are well-defined and accessible via a state-preparation oracles. Discretizing 136 by \(Q\) in \(k\) and \(P\) in \(s\), and then replacing each \(U_{k_i}(\tau_p)\) by its MPF approximation, gives \[\label{eq:duhamel-LCHS-MPF-double-sum} \frac{1}{\sqrt{2\pi}}\int_0^T\int_{\mathbb{R}}\hat{f}_{a,b}(k;c,d)U_k(T-s)b(s)\,dk\,ds \approx \sum_{p\in\mathcal{I}_P}\sum_{i\in\mathcal{I}_Q}\omega_p v_iU^{(2m)}_{k_i}(\tau_p)b(s_p).\tag{138}\] Equivalently, the non-homogeneous LCHS–MPF approximation is \[\label{eq:non-homogeneous-LCHS-MPF-approx} u_{Q,P}^{\operatorname{MPF}}(T):=W_Q^{\operatorname{MPF}}(T)u(0)+\sum_{p\in\mathcal{I}_P}\omega_pW_Q^{\operatorname{MPF}}(\tau_p)b(s_p).\tag{139}\] We state the error estimate in terms of two quadrature errors. First, define the uniform-in-time LCHS quadrature error by \[\label{eq:uniform-quadrature-error} E_{\operatorname{quad}}^{[0,T]}(Q):=\sup_{0\leq \tau\leq T}\left\|\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\hat{f}_{a,b}(k;c,d)U_k(\tau)\,dk-W_Q^{\operatorname{ideal}}(\tau)\right\|.\tag{140}\] The quadrature estimates in 4.6 also bound \(E_{\operatorname{quad}}^{[0,T]}(Q)\). Indeed, the strip estimates used there are monotone in the final time. Replacing \(T\) by any \(\tau\in[0,T]\) only decreases the factor involving \(T\|L\|\). Hence the same quadrature rule \(Q\), chosen using the endpoint time \(T\), satisfies \[E_{\operatorname{quad}}^{[0,T]}(Q)\le E_{\operatorname{mesh}}(Q;T)+E_{\operatorname{tail}}(Q;T),\] with the right-hand side given by the corresponding estimates in 35 or 38. Second, define the time-quadrature error by \[\label{eq:time-quadrature-error} E_{\operatorname{time}}(P):=\left\|\int_0^T e^{-A(T-s)}b(s)\,ds-\sum_{p\in\mathcal{I}_P}\omega_p e^{-A(T-s_p)}b(s_p)\right\|.\tag{141}\] The quantity \(E_{\operatorname{time}}(P)\) is kept abstract below. Once \(P\) and the regularity of the representative \(b\) are specified, it can be bounded by the corresponding quadrature estimate.

We now combine the three error sources in the non-homogeneous construction. The homogeneous propagator error is needed uniformly for all \(\tau\in[0,T]\), since the source term contains propagators \(e^{-A(T-s)}\). 50 separates this uniform LCHS–MPF error from the additional time-quadrature error used to discretize the source term.

50 gives an abstract error decomposition. We next choose the time-quadrature, approximation, and post-quadrature error budgets so that their total contribution is at most the target accuracy \(\epsilon\). This gives the fixed-order non-homogeneous analogue of 26.

We now describe the implementation cost. Assume state-preparation oracles for \(\ket{\widetilde{u}(0)}:=u(0)/\|u(0)\|\), when \(u(0)\neq0\), and for \(\ket{\widetilde{b}(s_p)}:=b(s_p)/\|b(s_p)\|\), for each \(p\in\mathcal{I}_P\) with \(b(s_p)\neq0\). Terms with zero vector coefficients are omitted from the outer LCU, and the phases of \(\omega_p\) are included in the outer preparation. Let \(\beta_P:=\|u(0)\|+B_P\). Define the source-preparation oracle \[\operatorname{PREP}^{P}\ket 0 = \frac{1}{\sqrt{\beta_P}}\left(\sqrt{\|u(0)\|}\ket{\star}+\sum_{p\in\mathcal{I}_P}\sqrt{|\omega_p|\,\|b(s_p)\|}\ket p\right),\] with the phases of \(\omega_p\) included in the corresponding left-preparation state. The vector in 139 is then implemented by an outer LCU over \(\{\star\}\cup\mathcal{I}_P\), with normalization \(\beta_P\). Since each term applies a block-encoding of \(W_Q^{\operatorname{MPF}}(\tau)\) with normalization \(\alpha_Q\), the total normalization is \(\alpha_Q\beta_P\). Thus, the amplitude amplification factor is \(\chi_{Q,P}(u(0),b):=\frac{\alpha_Q\beta_P}{\|u_{Q,P}^{\operatorname{MPF}}(T)\|}\). 52 gives the complexity estimates.

We see that no new commutator quantities arise. The source \(b(t)\) affects only source-state preparation, the time-quadrature error \(E_{\operatorname{time}}(P)\), and the normalization \(B_P\). The post-quadrature quantities \(R_Q\), \(\alpha_Q\), \(\Lambda_{m,Q}\), \(\Phi_{m,Q}^{\ast}\), \(r_Q\), \(\mu_Q\) are inherited from the homogeneous analysis.

9.2 Time-Dependent Case↩︎

We now extend the analysis to the time-dependent case for the special case local and extensive Hamiltonians. For simplicity, we consider the homogeneous problem. \[\label{eq:td-homogeneous-problem} \frac{d}{dt}u(t) = -A(t)u(t), \qquad A(t) = L(t) + iH(t),\tag{143}\] where, for every \(t \in [0,T]\), the matrices \(L(t)\) and \(H(t)\) are Hermitian and \(L(t) \succeq 0\). The extension to the non-homogeneous case follows as in 9.1. Let \(U_A(s,t) := \mathcal{T} \exp (-\int_s^t A(\tau) d\tau )\) denote the non-unitary propagator from \(s\) to \(t\). As before, for \(k \in \mathbb{R}\), define \(G_k(t) := H(t) + kL(t)\) and let \[U_k(s,t) := \mathcal{T} \exp \left(-i\int_s^t G_k(\tau)d\tau \right).\] The approximation and truncation estimates from [approx-error-time-indept] 13 continue to apply after replacing \(e^{-AT}\) by \(U_A(0,T)\) and \(U_k(T)\) by \(U_k(0,T)\). In particular, for the kernel profile \(\vec{\theta}=(a,b,c,d)\), the same admissible choices of \(d\) and \(y_0\) from 11 12 continue to hold. Motivated by 6, define \[C_m:=4^{2m+1}, \quad R_m(\delta):=\|\vec{a}\|_1\|\vec{b}\|_1\delta, \quad \Phi_m(k,N,\delta) := \mu_{k,p_0}^{2m+1}, \quad \mu_{R,p_0}:=\sup_{|k|\le R}\mu_{k,p_0},\] where \(\mu_{k,p_0}\) denotes the quantity \(\mu_{p_0}\) from 6 applied to the Hamiltonian \(G_k(t)\), \(p_0\) is defined as in 6, and \(\|\vec{a}\|_1\) and \(\|\vec{b}\|_1\) are the \(1\)-norms of the MPF coefficient vectors discussed above. Throughout this subsection, we assume that the family \(\{G_k(t): |k|\le R\}\) admits a decomposition into \(q_R\)-local, \(g_R\)-extensive terms, and that the derivative bound in 6 holds with a constant \(f_R\), uniformly in \(k\) and in the initial time of the sub-interval. Define \[\label{eq:td-admissible-step-size} \Delta(R,N,\delta):=\min\left\{\frac{1}{8e^3p_0(2q_Rg_R+2\Gamma f_R)}, \frac{1}{8\mu_{R,p_0}} \right\}.\tag{144}\] We first record the pointwise short-time estimate obtained from 6.

The passage from short-time to long-time simulation is not identical to the time-independent case. This is because the exact propagator on \([0,T]\) is an ordered product of generally different interval propagators \(U_k(t_j,t_{j+1})\), rather than a power of a single fixed short-time propagator. Let \(t_j=j\Delta\), where \(\Delta=T/r\). Then \[\label{eq:td-exact-ordered-product} U_k(0,T)=\overset{\longleftarrow}{\prod_{0 \leq j \leq r-1}} U_k(t_j,t_{j+1}).\tag{145}\] The factors in 145 are generally distinct, so one cannot use the identity \(U_k(0,T)=U_k(0,\Delta)^r\). We instead use the following telescoping estimate for non-stationary products.

We now pass to the post-quadrature setting. The construction remains the same as before, but we recall the relevant details for completeness. Let \(Q=\{(k_i,w_i)\}_{i\in\mathcal{I}_Q}\) be a quadrature rule, and define \(v_i\), \(\alpha_Q\), and \(R_Q\) as in 21 . Similarly, define \(\Lambda_{m,Q}\) as in 22 , and set \(\Phi_{m,Q}^{\ast} := \max_{i\in\mathcal{I}_Q}\Phi_m(k_i,N,\delta)\). The ideal and MPF-implemented quadrature operators, \(W_Q^{\operatorname{ideal}}\) and \(W_Q^{\operatorname{MPF}}\), respectively, are defined as above. As in 21, we have \[\begin{align} \label{eq:td-postquad-inner-error} \| W_Q^{\operatorname{ideal}}-W_Q^{\operatorname{MPF}} \| &\le \left(1+\eta_Q(T/r)\right)^{r-1} \left( \|\vec{a}\|_1C_m \frac{T^{2m+1}}{r^{2m}} \Lambda_{m,Q} + r\alpha_QR_m(\delta) \right), \end{align}\tag{146}\] where \(\eta_Q(\Delta) := \max_{i \in \mathcal{I}_Q}\left(\|\vec{a}\|_1 C_m\Delta^{2m+1}\Phi_m(k_i,N,\delta) + R_m(\delta)\right)\). Define the quadrature error \(E_{\operatorname{quad}}(Q)\) as before. Combining the approximation, quadrature, and inner simulation errors gives \[\begin{align} \label{eq:td-total-error-bound} \left\| U_A(0,T)-W_Q^{\operatorname{MPF}} \right\| &\le E_{\operatorname{approx}}(y_0) + E_{\operatorname{quad}}(Q) \\ &\quad+ \left(1+\eta_Q(T/r)\right)^{r-1} \left( \|\vec{a}\|_1C_m \frac{T^{2m+1}}{r^{2m}} \Lambda_{m,Q} + r\alpha_QR_m(\delta) \right). \end{align}\tag{147}\] The implementation is the same as in 4.2, except that the controlled operation for a fixed quadrature node \(k_i\) implements the ordered product \(U_{k_i}^{(2m)}(0,T) = \overset{\longleftarrow}{\prod_{0\leq j \leq r-1}} \widetilde{U}_{k_i,j}^{(2m)}\) Thus one controlled implementation of \(U_{k_i}^{(2m)}(0,T)\) still uses \(\mathcal{O}(rK_m)\) controlled second-order product-formula queries, provided each time-dependent second-order product-formula segment is available with the same oracle cost model as in 6. The complexity of the algorithm follows from the preceding analysis in 4.3 and 4.5. We do not repeat those details here.

References↩︎

---

1. junaida@umd.edu↩︎

2. Consequently, the symbol \(p\) will be used freely for other purposes.↩︎

3. In the time-dependent case, we say that \(H\) is \(g\)-extensive if \(\sup_{\tau\in[0,\Delta]}\sum_{\gamma:\,j\in X_\gamma}\|H_\gamma(\tau)\|\le g\) for all \(j\).↩︎