Hamiltonian simulation for hyperbolic partial differential equations by scalable quantum circuits

Solving partial differential equations for extremely large-scale systems within a feasible computation time serves in accelerating engineering developments. Quantum computing algorithms, particularly the Hamiltonian simulations, present a potential and promising approach to achieve this purpose. Actually, there are several oracle-based Hamiltonian simulations with potential quantum speedup, but their detailed implementations and accordingly the detailed computational complexities are all unclear. This paper presents a method that enables us to explicitly implement the quantum circuit for Hamiltonian simulation; the key technique is the explicit gate construction of differential operators contained in the target partial differential equation discretized by the finite difference method. Moreover, we show that the space and time complexities of the constructed circuit are exponentially smaller than those of conventional classical algorithms. We also provide numerical experiments and an experiment on a real device for the wave equation to demonstrate the validity of our proposed method.


I. INTRODUCTION
Partial differential equations (PDEs) serve as essential tools for investigating the dynamic behavior of various physical phenomena, including heat conduction, fluid dynamics, and electromagnetic waves [1].Solving PDEs for extremely large systems within a reasonable computation time is crucial for accelerating engineering developments in industries.Despite remarkable progress in addressing extensive physical systems through the use of supercomputers [2,3], obtaining solutions within a feasible computation time is still intractable.
A potentially promising strategy to substantially reduce the computational expenses for solving PDEs involves the utilization of quantum computing.Quantum computing has attracted considerable interest in recent decades as a prospective avenue for achieving dramatically fast computation compared to classical computing.Although quantum computers currently suffer from limited hardware scalability and less noise resistance, there has been remarkable progress in hardware performance.One of the promising applications of quantum computers is a solver of PDEs.
For steady-state problems, PDEs reduce to a system of linear or non-linear equations and can be solved by a linear system solver.There are mainly two types of quantum algorithms for solving systems of linear equations: one is variational quantum algorithms [4][5][6] and the other is the Harrow-Hassidim-Lloyd (HHL) algorithm [7].Variational quantum algorithms are aimed at the use on near-term quantum devices and have been extensively studied for applications to PDEs [8][9][10][11].On the other hand, HHL-based algorithms focus on the fault-tolerant quantum computers and provide theoretical quantum speedup over classical algorithms under certain conditions [7,12,13].Although HHL-based algorithms require several oracles to be implemented for state preparation, matrix inversion, and extracting solutions [14], there are various studies that can be applied for implementing each part [15][16][17].
For time evolution problems governed by PDEs, there are mainly two types of quantum algorithms as well, i.e., the near-term and long-term algorithms.As for nearterm algorithms, variational quantum simulation [18,19] has been applied to solve PDEs [20,21] while Hamiltonian simulation [22] is a counterpart for long-term ones.
Hamiltonian simulation involves implementing a quantum circuit for the time evolution of a quantum system, exp(−iHτ ), with the time increment τ and the Hamiltonian of the target quantum system H; it is also referred to as quantum simulation.The targeted Hamiltonians are, for example, Ising Hamiltonians [23,24] and molecular Hamiltonians [25].Remarkably, implementation of quantum simulation based on the Ising Hamiltonian on a 127-qubit quantum computer has been recently reported, which is a significant contribution demonstrating the utility of quantum computers [23].On the other hand, there are also reports of applying Hamiltonian simulation by reducing the governing equations of classical systems to the Schrödinger equation [26][27][28][29][30][31][32].Costa et al. [26] proposed quantum simulation of the wave equation where the Hamiltonian is given as the incidence matrix of a graph that represents the discretized target space.Babbush et al. [27] generalized this approach and proposed a quantum algorithm for simulating classical coupled oscil-lators with the rigorous proof of the exponential speedup over any classical algorithms.Jin et al. [28,29] proposed an interesting approach called Schrödingerization where the general ordinary differential equation is transformed to the Schrödinger equation by the warped phase transformation.An et al. [30] proposed a technique of linear combination of Hamiltonian simulation for simulating general nonunitary dynamics.Although these results suggest that quantum algorithms may exhibit exponential speedup even for classical system simulations, these methods rely on an oracle access to the Hamiltonian, which makes their implementation by elementary quantum gates unclear.
This paper proposes a Hamiltonian simulation algorithm for solving a special type of PDEs, i.e., linear hyperbolic PDEs without any source terms, which can be transformed into the Schrödinger equation.Specifically, we derive an explicit quantum circuit representation of time evolution operators for Hamiltonian simulation driven by differential operators.Then we apply the proposed method to the advection and wave equations, which are transformed into the Schrödinger equation and are discretized by the finite difference method (FDM).
The key technique of implementing the time evolution operators given by differential operators is to diagonalize each term of Hamiltonian using the Bell basis.This is similar to the idea of the extended Bell measurement [33] that efficiently estimates the expectation of band matrices derived from the discretization of PDEs.
The contributions of this paper are summarized as follows.
• We provide an algorithm to build quantum circuits of time evolution operators for Hamiltonian simulation by differential operators.This result will contribute to implementing the FDM on a quantum computer.
• We derive the space and time complexities (Theorems 1 and 2); our implementation requires dnqubits and quantum circuits with O(dn 3 T 2 /ε) or O(dn 2.5 T 1.5 /ε 0.5 ) non-local gates, to perform Hamiltonian simulation of hyperbolic PDEs defined on the d-dimensional lattice with 2 n nodes in each dimension up to time T within the additive error ε.
That is, the algorithm enjoys an exponential reduction of resources with respect to the spatial degree of freedom.
• We transform the advection and wave equations into the Schrödinger equation in the real space for solving them by Hamiltonian simulation, based on the proof of self-adjointness of the Hamiltonian derived from these equations.This analysis is followed by thorough numerical simulations and an experiment on a real quantum device.
The rest of this paper is organized as follows.In Sec.II, we briefly introduce the finite difference operators along with their representation in a qubits system.We also discuss the transformation of the advection and wave equations into the Schrödinger equation to make our method applicable to these equations.In Sec.III, we then provide the quantum circuit implementation of the finite difference operators, together with the theoretical error bound between the constructed circuit and the target Hamiltonian evolution.In Sec.IV, we provide several numerical experiments and an experimental result of a real device.Finally, we conclude this study in Sec.V.

A. Finite difference operators and its representation in qubits
Let us consider a one-dimensional domain Ω := (0, L) where L is the length of the domain.We discretize the domain Ω by uniformly distributed N points with the interval l := L/(N + 1), where N is a power of two, i.e., N = 2 n for n ∈ N. Let us then consider a scalar field u defined on the domain Ω and discretize the field u using its value at each point (i.e., at each node in the 1-dimensional lattice) as u := [u 0 , u 1 , . . ., u N −1 ].The forward difference operator D + for the spatial derivative acts on u as where u N is determined from the boundary condition (BC).For instance, for periodic BC. ( Note that the prescribed value for the Neumann boundary condition is set to ensure that the spatial derivative (D + u) N −1 on the boundary node is zero.Similarly, the backward difference operator D − is defined as the operator that acts on u as where u −1 is determined from the BC such as The central difference operator D ± and the Laplacian operator D ∆ are the operators acting as where u N and u −1 are defined by Eqs. ( 2) and ( 4), respectively.
In what follows we quantize the above-described difference operators.For this purpose, we need the quantum state corresponding to the discretized field u, on which those quantized operators act; that is, let |u⟩ be an nqubit state on which u is encoded as where |j⟩ := |j n−1 j n−2 . . .j 0 ⟩ with j n−1 , j n−2 , . . ., j 0 ∈ {0, 1} is the computational basis.Here, we assume that u is normalized, i.e., ∥u∥ 2 = 1.Now, it is known that, for the qubit-based system, the finite difference operators can be represented as matrix product operators (MPOs) using the following three 2 × 2 matrices [34]: where σ 01 and σ 10 are the ladder operators.The following two 2 × 2 matrices are also useful: The point of quantization is the MPO representation of the shift operators S − = j=1 |j⟩ ⟨j − 1| as follows: where Here, σ ⊗0 ij is regarded as a scalar 1.For example, for the case n = 2, we have S − = |0⟩ ⟨1| + |1⟩ ⟨2| + |2⟩ ⟨3| = σ 01 ⊗ σ 10 + I ⊗ σ 01 .Using the shift operator, we obtain the forward difference operator with the Dirichlet BC as follows: Actually, the following relationship holds: where u 2 n = 0, which exactly corresponds to the forward difference operator with the Dirichlet BC.To impose the Neumann BC, it suffices to add σ ⊗n 11 /l to D + D as For the periodic BC, we have to add σ ⊗n 10 /l to D + D as Similarly, the other finite difference operators can be represented as Since we have the shift operators S − and S + , we can also construct the finite difference operators of the higherorder approximations as discussed in Appendix A. Quantization of the difference operators for the ddimensional domain Ω = (0, L) d is straightforward, as follows.We discretize each segment (0, L) by uniformly distributed N points with the interval l := L/(N + 1), which results in a d-dimensional lattice with N d nodes.We again assume that N is a power of two, i.e., N = 2 n for n ∈ N. Let |u(t)⟩ be a dn-qubit state on which the discretized field u on the lattice is encoded as where |j α ⟩ := |(j α ) n−1 (j α ) n−2 . . .(j α ) 0 ⟩ with (j α ) n−1 , (j α ) n−2 , . . ., (j α ) 0 ∈ {0, 1} is the computational basis and x jα is the spatial coordinate of the j α -th node along the x α -axis.We again assume that u is normalized, i.e., ∥u∥ 2 = 1.The finite difference operator D µ B defined above can easily be extended to those for |u⟩ in the d-dimensional space, as follows: where µ ∈ {−, +, ±, ∆} and B ∈ {D, N, P}.

B. Transforming partial differential equations to Schrödinger equation
In this study, we particularly focus on the advection equation and the wave equation as hyperbolic partial differential equations.To apply Hamiltonian simulation for these equations, we first need to express them in the form of the Schrödinger equation.

Advection equation
Let the scalar field u be governed by the advection equation with the constant velocity field v as where t is the time, x is the spatial coordinate and ∇ is the spatial differential operator; that is, ∇ = ∂/∂x for the one-dimensional case d = 1.The scalar field u represents, for example, temperature or concentration.
The advection equation can be rewritten as If the operator −iv • ∇ is self-adjoint (meaning that ⟨ũ, (−iv • ∇) u⟩ = ⟨(−iv • ∇) ũ, u⟩ holds for arbitrary scalar fields u and ũ), Eq. ( 24) is exactly the Schrödinger equation; actually, this property holds under an appropriate BC as shown in Appendix B 1. Thus, the advection equation falls into the Schrödinger wave equation with the Hamiltonian that is expressed in terms of differential operators, ∇.The next step is to discretize the field variable u and the Hamiltonian −iv •∇.That is, u is discretized as u := [u 0 , u 1 , . . ., u N −1 ] and each element is encoded into the amplitude of the quantum state |u(t)⟩ given in Eq. (21).Also, the Hamiltonian −iv • ∇ is discretized using the central difference operator (19) as Since (D ± ) † = −D ± , the Hamiltonian with the central difference operator is actually a Hermitian matrix.As a result, we obtain the Schrödinger equation d |u(t)⟩ /dt = −iH |u(t)⟩, and this is what we aim to simulate on a quantum device using the Hamiltonian simulation method.Since the initial condition for the advection equation is given as specifying u(0, x), the state preparation oracle for the Hamiltonian simulation has to prepare a quantum state of |u(0)⟩ = Note that the upwind differencing scheme is preferable for the numerical stability for the advection equation [35], but using D + or D − to discretize ∇ does not retain the Hermitian property of H [36].In fact, discretizing ∇ using D + makes the Hamiltonian H = −iv 1 D + for d = 1, which is no longer Hermitian because H † = iv 1 (D + ) † = −iv 1 D − ̸ = H with the relationship of (D + ) † = D − .This is the reason why we here use the central difference scheme; we will investigate the relationship between finite difference schemes and the resulting Hermitian property to explore the possibility of using D + or D − , in our future research.

Wave equation
Let u be governed by the wave equation with the speed c as The scalar field u corresponds to, for example, the displacement of string and membrane, and the pressure.
The wave equation can be rewritten as by setting, when d = 1, and when d = 3, where d is the number of spatial dimensions and x α is the spatial coordinate.This form of the Hamiltonian is similar to that in Ref. [26].Note that when d = 2, we can also set ψ(t, x) and H by a 3 × 1 vector and a 3 × 3 matrix, respectively, which are obtained by omitting the row and column regarding x 3 -axis of those when d = 3.However, Eq. ( 29) has the advantage in quantizing these quantities because of the size of power of two.We can show that, for the above three cases, H is selfadjoint (meaning that ⟨ ψ, Hψ⟩ = ⟨H ψ, ψ⟩ holds for arbi-trary vector fields ψ and ψ) under appropriate BCs; the proof is given in Appendix B 2. Therefore, Eq. ( 27) is exactly the Schrödinger wave equation; in other words, the wave equation falls into the Schrödinger wave equation with its Hamiltonian containing the differential operators.
The scalar field u(t, x) is discretized and encoded in a quantum state |ψ(t)⟩ on qubits as Note that the initial quantum state |ψ(0)⟩ has to be implemented to satisfy the initial conditions of the wave equation, which are typically the initial conditions of u(0, x) and ∂u(0, x)/∂t.When we have an initial value of u(0, x) as an analytic function, which is a typical case in solving PDEs, we can calculate ∂u(0, x)/∂x for the initial condition of our algorithm.Thus, we assume that we can set the initial condition for the quantum state (i.e., ∂u(0, x)/∂x and ∂u(0, x)/∂t) based on the typical initial conditions of the wave equation (i.e., u(0, x) and ∂u(0, x)/∂t).Also, the Hamiltonian H is discretized by the forward and backward difference operators so that the resulting H can be actually a Hermitian matrix, as follows: As a result, we obtain the Schrödinger equation d |ψ(t)⟩ /dt = −iH |ψ(t)⟩, which can be simulated on a quantum device.Note that the above discretization imposes the Dirichlet BC for ∂u/∂t at x α = 0 and for ∂u/∂x α at x α = L.That is, the mixed BC for u consisting of the Dirichlet BC for x α = 0 and the Neumann BC for x α = L is imposed, which is consistent with the condition of self-adjointness of the Hamiltonian H as discussed in Appendix B 2. We can also use the central difference scheme for periodic BC because it retains the Hermitian property.We would like to conduct our future research to discuss the scheme for implementing arbitrary BCs while keeping the Hermitian property of the Hamiltonian.

III. CIRCUIT IMPLEMENTATION METHOD FOR PDES A. Quantum circuit for time evolution by differential operators
Now, we consider a hyperbolic partial differential equation such that it can be reduced to the Schrödinger equation d |u(t)⟩ /dt = −iH |u(t)⟩, where the Hamiltonian H consists of the difference operators, as exemplified in Eqs. ( 25) and (32) in Section II B. Note that, in addition to those motivating examples, a wide class of partial differential equations also falls into our target based on the technique proposed in Ref. [28,29].Our goal is to provide an efficient method for implementing the time evolution operator exp(−iHτ ) for a time increment τ , on a circuit of qubit-based quantum devices.
Let us consider the following Hamiltonian that con-tains the one-dimensional spatial difference operator: where γ ∈ R is a scale parameter and λ ∈ R is a phase parameter.This Hamiltonian can represent the essential part of the one-dimensional Hamiltonian for the advection equation given in Eq. ( 25).Also, the following procedure for constructing quantum circuits and the scaling of the circuit complexity, based on the Hamiltonian (33), is applicable to the Hamiltonian for the wave equation given in Eq. (32).Since σ 01 = (X + iY )/2 and where Z is the single-qubit Z gate.Here we call (|0⟩|1⟩ ⊗(j−1) ± e −iλ |1⟩|0⟩ ⊗(j−1) )/ √ 2 the Bell basis; the Bell basis is the key for decomposing the Hamiltonian into a sum of polynomial number of terms.Also, U j is the unitary matrix so that where H j is the Hadamard gate acting on the j-th qubit, P j (λ) is the Phase gate acting on the j-th qubit as and CNOT j m is the CNOT gate acting on the m-th qubit controlled by the j-th qubit.Note that we herein use the little endian.Applying the first-order Lie-Trotter-Suzuki decomposition, we can approximate the time evolution operator exp(−iHτ ), as follows: where CRZ 1,...,j−1 ) is the multi-controlled RZ gate acting on the j-th qubit controlled by 1, . . ., (j − 1)th qubits.Note that from the second to the third line, exp(I ⊗ A) = I ⊗ exp(A) and exp(U AU † ) = U exp(A)U † were used.The unitary matrix W j is defined as Then, let V (γτ, λ) denote the approximated time evolution operator as We now have a concrete circuit implementation of this approximating unitary matrix V , as shown in Fig. 1.Note that such explicit form of implementation has not been reported in the previous proposals [28,29].Moreover, thanks to the explicit form of V , we can have a detailed evaluation on the approximation error between exp(−iHτ ) and V in the sense of operator norm.Together with the explicit circuit construction, the following lemma gives the approximation error.
Lemma 1.Consider the Schrödinger equation d |u(t)⟩ /dt = −iH |u(t)⟩ such that the Hamiltonian H is given by Eq. (33).The time evolution operator exp(−iHτ ) with the time increment τ can be approximated by the unitary V in Eq. (39), and its explicit circuit implementation is shown in Fig. 1.Moreover, the approximation error in the sense of the operator norm is upper bounded as Proof.We use the fact that the approximation error of the Lie-Trotter-Suzuki decomposition is upper bounded by the sum of operator norm of commutators of all terms contained in the Hamiltonian [37, Proposition 9].Hence our task is to evaluate the commutators of all terms of the Hamiltonian (33), as follows; the detailed calculation is given in Appendix C. For n > j > j ′ > 1, the terms of the shift operators e iλ s − j + e −iλ s + j and e iλ s − j ′ + e −iλ s + j ′ commute as For j > j ′ = 1, we obtain Thus, the terms of the Hamiltonian, γ(e iλ s − j + e −iλ s + j ), can be grouped into those for j > 1 and j = 1.Since the unitary does not change the operator norm, the Trotter error is upper bounded using the result of Ref. [37,Proposition 9], as follows: Since W j>2 and W j ′ >2 commute, we can easily obtain the second-order formula as The following lemma gives the approximation error between exp(−iHτ ) and V (2) in the sense of operator norm.
Lemma 2. Consider the Schrödinger equation d |u(t)⟩ /dt = −iH |u(t)⟩ such that the Hamiltonian H is given by Eq. (33).The time evolution operator exp(−iHτ ) with the time increment τ can be approximated by the unitary V (2) in Eq. (44).The approximation error in the sense of the operator norm is upper bounded as Proof.The detailed proof is given in Appendix D. Here, we provide a scketch of the proof.Let us group the terms of the Hamiltonian H in Eq. ( 33) into Since all terms in H 2 commute each other, the error of the second-order Suzuki formula is upper bounded [37, Proposition 9] by exp(−iHτ ) − V (2) (γτ, λ) By evaluating the commutators, we obtain We also provide the following lemma about the gate counts of the unitary V .Lemma 3. The approximated time evolution operators V in Fig. 1 and V (2) can be implemented using singlequbit gates and at most 9n 2 − 33n + 34 CNOT gates for n ≥ 3.
Proof.As shown in Fig. 1, the non-local gates included in the operator W j are a multi-controlled Rz gate for j ≥ 3 (a controlled Rz gate for j = 2) and totally 2(j−1) CNOT gates in the unitary U j and U † j for j ≥ 2. It is known that the multi-controlled RZ gate with (j − 1) control qubits can be decomposed into single-qubit gates and at most 16j − 40 CNOT gates [38,Theorem 3].Therefore, the number of CNOT gates required to implement the approximated time evolution operator V is n j=3 Since operators V and V (2) have the difference only in the single-qubit gate by definition in Eq. ( 44), the number of CNOT gates in V (2) is the same as that in V .
We now extend the above discussion to the operators acting on a d-dimensional domain Ω.Lemma 4. Let the Hamiltonian H consist of finite difference operators for a dn-qubit system as where γ ∈ R is a scale parameter, λ α ∈ R is the phase parameter and for µ ∈ {−, +}.The time evolution operator exp(−iHτ ) with the time increment τ can be approximated by the unitary The approximation error is upper bounded in the sense of the operator norm as Proof.Here we sketch the proof; the detail is given in Appendix C. From Eq. ( 34), we obtain where V (γη α τ, λ α ) is given in Eq. ( 39) and represents the time evolution operator for the spatial dimension in the x α direction.The approximation error of the Lie-Trotter-Suzuki decomposition is upper bounded by the operator norm of commutators of Hamiltonian [37, Proposition 9].Together with the fact the terms of the Hamiltonian We can easily extend this lemma to the second-order formula as follows.
To simulate the Hamiltonian dynamics over the total time T , it suffices to divide the total time T into r(:= T /τ ) intervals so that the approximation error occurred in the time interval τ could be small enough to be acceptable.
We remark that the essential part of the Hamiltonian (25) for the advection equation falls into the Hamiltonian (49) by setting γ = 1/(2l), η α = v α and λ α = −π/2.Appendix E 1 gives the procedure for constructing the quantum circuit and its Trotter error, for the full Hamiltonian of the advection equation.Also for the case of wave equation, although the Hamiltonian in Eq. ( 32) does not fit into that in Lemma 1, we can easily obtain the quantum circuit for Hamiltonian simulation as discussed in Appendix E 2.

B. Space and time complexities
Together with the quantization of the field in Eq. ( 21), the explicit circuit construction and Lemmas 4 and 6, we now provide the following theorem giving the space and time complexities of our Hamiltonian simulation, as our main result.
Theorem 1.Let H be the Hamiltonian as defined in Eq. (49), i.e., H = γ The time evolution operator exp(−iHT ) up to time T is implementable on the dn-qubits system using the quantum circuits with O(dn 3 γ 2 T 2 d α=1 η 2 α /ε) non-local gates within the additive error ε.Furthermore, the leading term of the number of non-local gates is 9dn 3 γ 2 T 2 d α=1 η 2 α /(2ε).The quantum circuit for exp(−iHT ) consists of the repetitive applications of the one time step unitary V shown in Fig. 1.
Proof.Lemma 4 states that the additive error of the approximated time evolution operator α /2 with the time increment τ .To suppress the error of the simulation over the total time T within a small value ε, it suffices to divide the total time T into r(:= T /τ ) intervals so that which is rearranged as Since each Trotter step As for the second-order formula, we also provide the following theorem giving the space and time complexities of our Hamiltonian simulation.
Note that since V r = W 1 (−γτ, λ)(V (2) ) r W 1 (γτ, λ), the first-order formula works similarly to the secondorder formula and can exhibit the similar bound to the second-order formula in a practical sense [39].
Remark 1. Classical implementation requires O(2 dn ) memories to store the discretized scalar field u.Since the finite difference operators, such as D µ B for µ ∈ {−, +, ±, ∆} and B ∈ {D, N, P} in Eqs.(1), (3), ( 5) and (6), can be represented by sparse matrices of the size 2 dn × 2 dn with the sparsity denoted by s, the application of the operator D µ B to the vector u requires O(s2 dn ) arithmetic operations.Here, s = 2 for µ ∈ {−, +, ±}, i.e., the first-order derivative and s = 3 for µ = ∆, i.e., the Laplacian in Eq. ( 6).Now, let us consider classical simulation using the forward Euler scheme, which has the additive error bounded by O(τ 2 ) where τ is the time increment.Then, classical simulation up to time T within the additive error ε requires O(T 2 /ε) steps, which results in O(s2 dn (T 2 /ε)) arithmetic operations.When using the second-order formula for the time integration, the complexity will be improved to O(s2 dn (T 1.5 /ε 0.5 )).In addition, the Courant-Friedrichs-Lewy (CFL) condition [40] requires T /r = O(l).As a result, classical simulation up to time T within the additive error ε requires O(s2 dn (T 2 /ε + T /l)) or O(s2 dn (T 1.5 /ε 0.5 + T /l)) arithmetic operations.
Remark 2. Focusing on the spatial dimension d, the number of nodes in each dimension 2 n , the total simulation time T , and the additive error of the simulation, ε, our Hamiltonian simulation requires O(dn 3 T 2 /ε) nonlocal quantum gate operations under the assumption of γ = O(1) and η α = O(1) while classical approaches require O(2 dn T 2 /ε) arithmetic operations.Considering the relationship of n = O(log(L/l)) and assuming γ = O(1/l), which is the case of the advection and wave equations, our Hamiltonian simulation requires O(dT 2 log(L/l) 3 /(l 2 ε)) or O(dT 1.5 log(L/l) 2.5 /(l 1.5 ε 0.5 )) under the assumption of These theorem and remarks suggest that our method has the potential for exponential speedup with respect to the number of nodes on the lattice N (= 2 n ) and the size of the domain L when simulating the classical dynamics governed by differential operators, if the operation of quantum gates can be performed as fast as the arithmetic operations in classical computers.Remark 2 also implies our method will exhibit the polynomial speedup with respect to the interval of the lattice l when d ≥ 2.
Note that the time complexity of our algorithm is polynomial with respect to the additive error ε while the state-of-the-art algorithms [41,42] give the time complexity of poly(log(1/ε)).Such algorithms rely on the oracle access to the block encoding of the Hamiltonian, while our current study focuses on deriving the explicit quantum circuit for Hamiltonian simulation.We would like to conduct our future work to derive the quantum circuit for the block encoding-based Hamiltonian simulation, providing the comparison of the gate counts of the Trotter-based Hamiltonian simulation and the block encoding-based one including the constant factor in our future work.
Finally, we remark that it is impractical to access all components of |u(t)⟩ or |ψ(t)⟩ for the advection and the wave equation, respectively, because it requires O(2 dn ) measurements.Hence, the proposed method should be used in a situation when only some characteristic quantities about the solution are of interest; typically, such quantity is represented by ⟨u(t)|O|u(t)⟩ or ⟨ψ(t)|O|ψ(t)⟩ for an observable O.As we discuss later, one example of such observables is the kinetic energy of the system governed by the wave equation: which leads to ⟨ψ(t)|O|ψ(t)⟩ = x |∂u(t, x)/∂t| 2 .The power spectra of the system is also one of the possible observables [32].We can use various well established methods for such estimation of observables [43,44].We would like to construct meaningful observables in more detail for specific applications and discuss efficient estimation of their expected values in our future work.

IV. NUMERICAL AND EXPERIMENTAL RESULTS
In this section, we provide several results of numerical experiments to demonstrate the validity of our proposed method.We used Qiskit 0.45 [45], the open-source toolkit for quantum computation to implement quantum circuits.

A. Advection equation
We first performed Hamiltonian simulation for solving the advection equation with periodic BC, described in Section II B 1. Here, we compare the simulation result of our proposed method with the following two; the first one is the results calculated by directly applying the matrix exponential operator e −iHτ to |u(0)⟩ by matrix-vector calculation, and the second one is the result calculated by the fully classical finite difference method (FDM) with the central differencing scheme given in Eq. ( 5).
Figure 2 illustrates the simulation results of the advection equation with periodic BC in one dimension.We set n = 7, τ = 0.1, T = 20, l = 1, and v 1 = 1 to formulate the problem and used statevector simulator to run the quantum circuits.Also, we set the time increment parameter to 0.01 for the FDM.As an initial state, we prepare Figure 2 also include this initial state as the amplitude at t = 0 to visualize the dynamics from the initial state.Solid lines in Fig. 2(a) represent the amplitudes of the quantum state |u(t)⟩ generated by the proposed circuits, at time t = 0, 10, 20; note that in this case the probability amplitudes of |u(t)⟩ are always real.Dashed lines represent the amplitudes of the quantum state |u(t)⟩ which is directly generated by applying the matrix exponential operator e −iHτ to |u(0)⟩ via matrix-vector calculation.Dotted lines represents the solutions calculated by the FDM.Since the velocity is chosen to have a positive value v 1 = 1, we observe that the scalar field u(t, x i ) moves toward the positive direction as time passes.Also, the solid lines in Fig. 2(b) illustrate the absolute errors between the amplitudes of the quantum state |u(t)⟩ by the proposed circuits and those calculated by FDM; moreover, the dashed lines illustrate the absolute errors between the amplitudes calculated via applying the matrix exponential operator e −iHt to |u(0)⟩ via matrix-vector calculation and those calculated by FDM.The figures imply that the solutions obtained by the three approaches agree well, which demonstrated that the advection equation discretized by the FDM can be simulated by the Hamitlonian simulation algorithm on quantum circuits under acceptable Trotter errors.The oscillation occurred in some spatial region where u(t, x i ) had the sharp gradients, is due to the central differencing scheme, and it is well-known that it can be resolved by using the upwind differencing scheme [35].However, such differencing scheme does not retain the Hermitian property of the Hamiltonian as we mentioned in Section II B 1, that is, (1) and ( 3) to discretize ∇ because of the relationship of (D + ) † = D − .We will address this issue in the future research.Next, Fig. 3 shows the simulation results of the advection equation with periodic BC in two dimensions.We set n = 6, τ = 0.1, T = 20, l = 1, and to setup the problem and used statevector simulator, to run the quantum circuits.The number of qubits is 2n = 12.We chose the time increment parameter to be 0.01 for the FDM.As an initial state, we prepare Figure 3 also include this initial state as the amplitude at t = 0 to visualize the dynamics from the initial state.The left column in Fig. 3 represents the amplitudes of the quantum state |u(t)⟩ generated by the proposed circuits described in Section II B 1. The center column represents the solutions directly obtained by applying the matrix exponential operator e −iHτ to |u(0)⟩ by matrixvector calculation.The right column represents solutions calculated by the FDM.Each row represents the solutions at t = 0, t = 10 and t = 20, respectively.We observe that the distribution of the scalar field |u⟩ moves toward upper right direction according to the setting of Although several numerical oscillation occurred due to the central differencing scheme, we again find that the solutions obtained by the three approaches agree well, which demonstrate the validity of the proposed method.Since the solution in the center column is obtained by directly applying the matrix exponential operator, it does not include the error in the numerical time integration.The difference between the left and center columns comes from the Trotter error, while the difference between the center and right columns comes from the time integral by the forward Euler scheme for the FDM.These errors can be decreased by using a smaller time increment parameter.That is, the time evolution operator exp(−iHτ ) is directly applied to the (a) (b) FIG. 2. Simulation results of the advection equation with periodic BC in one dimension.(a) Solid lines represent the amplitudes of the quantum state |u(t)⟩ generated by the proposed circuits.Dashed lines also represent the amplitudes of the quantum state |u(t)⟩, but it was directly computed by applying the matrix exponential operator e −iHt to |u(0)⟩ via matrix-vector calculation.Dotted lines represent the solutions calculated by the finite difference method (FDM).These three lines are overlapped very well in almost all the spatial region.(b) Absolute errors of the amplitudes of the quantum state |u(t)⟩ by the proposed circuits (solide lines) and those applying the matrix exponential operator e −iHt to |u(0)⟩ via matrix-vector calculation (dashed lines), to the solutions calculated by FDM, respectively.state |u(t)⟩ in the center column, while the forward Euler scheme [35] is used to proceed time in the right column.

B. Wave equation
Next, we show Hamiltonian simulation for solving the wave equation with the mixed BC in one dimension and with the periodic BC in two dimensions.We here again compare the simulation results of our proposed method with those obtained by directly using the matrix exponential operator e −iHτ to |ψ(0)⟩ and those by the fully classical finite difference method (FDM) with the central differencing scheme.
Figure 4 illustrates the simulation results of the wave equation with the mixed BC in one dimension, specifically the Dirichlet BC for the left side and the Neumann BC for the right side.We set n = 4, τ = 0.1, T = 20, l = 1 and c = 1 to setup the problem and used statevector simulator to run quantum circuits.The number of qubits for encoding this problem is n + 1 = 5.We use a time increment parameter 0.1 for the FDM.As an initial state, we prepare which corresponds to the initial condition u(0, x j ) = 0 (65) Solid, dashed, and dotted lines represent components of solutions corresponding to ∂u(t, x j )/∂t prepared by the proposed circuits described in Section II B 2, the direct multiplication of the matrix exponential operator by |ψ(0)⟩, and the FDM, respectively.Figure 4 clearly illustrates that the solutions obtained by the three approaches agree well, demonstrating that the wave equation could be simulated as the Schrödinger equation as well and actually be implemented on quantum circuits under acceptable Trotter errors.Figure 5 illustrates the simulation results of the wave equation with the periodic BC in two dimensions.We set n = 6, τ = 0.1, T = 20, l = 1 and c = 1 to setup

Proposed
Matrix exponential FDM  the problem and used statevector simulator to run quantum circuits.The number of qubits to implement this problem is 2n + 1 = 13.Here, we used the central difference operator D ± to define the Hamiltonian.We use a time increment parameter 0.1 for the FDM.As an initial state, we prepare

Proposed
Matrix exponential FDM which reflectes the initial condition of u(t, x j1 , x j2 ) = 0 (68) Figure 5 also include this initial state as the amplitude at t = 0 to visualize the dynamics from the initial state.The left and center columns in Fig. 5 represent components of solutions corresponding to ∂u(t, x j1 , x j2 )/∂t obtained by the proposed quantum circuit and the di-rect calculation of the matrix exponential operator, respectively.The right column represents solutions calculated by the FDM.Each row represents the solutions at t = 0, t = 10 and t = 20, respectively.The results of left and center columns clearly coincide, implying that the proposed quantum circuit could accurately simulate the Schrödinger equation derived from the wave equation.However, these result slightly differ from that by the FDM; the results by the proposed method look rougher than those by the FDM.This comes from the difference of finite difference schemes.For the fully classical FDM, we discretized the Laplacian ∇ 2 of the wave equation by D ∆ while for the Schrödinger equation derived from the wave equation, we discretized the gradient operator ∇ by D ± , which corresponds to discretizing the Laplacian ∇ 2 by (D ± ) 2 .Thus, the discretized equations do not exactly match and accordingly have errors that can be reduced by decreasing the spatial interval l.Although there exists such numerical error, overall, the results well grasp the property of wave propagation.Finally, we provide an experimental result conducted on a real quantum device.We here again study the wave equation with the mixed BC in one dimension setting n = 2, τ = 0.2, T = 2, l = 1 and c = 1.The number of qubits for encoding this problem is n + 1 = 3.Although the number of qubits used here is quite small due to the limitation of noise resilience of current hardware, we show the experimental result to demonstrate that our quantum algorithm is implementable on a current device owing to the explicit construction of quantum circuits by elementary gates of the hardware.We compare the simulation results obtained from a real device, statevector simulator, and the direct multiplication of the matrix exponential operator e −iHτ by |ψ(0)⟩.Although we illustrated all probability amplitudes of quantum states prepared by each approach in the preceding results for validation, accessing all amplitudes is unrealistic for practical use because of the need of the quantum state tomography as we mentioned in Section III B. Here, we evaluated the expectation value of the following specific observable: which is the kinetic energy of the system.As an initial state, we simply used which corresponds to the initial condition u(0, x j ) = 0 (72) ∂u(0, x j ) ∂t = 1 for j = 1 0 otherwise.(73) Figure 6 illustrates the experimental results.We used 4000 shots to estimate the expectation ⟨ψ(t)|O|ψ(t)⟩ at each time.We used the dynamical decoupling (DD) techniques [46][47][48] with super-Hahn echo [46,49] to suppress the decoherence, and also used the readout error mitigation technique, specifically Twirled Readout Error eXtinction (TREX) [50].The number of non-local gates in the quantum circuit was 120 at the maximum simulation time t = 2 for this problem.Note that error bars representing the 95% confidence interval of the expectation value under the shot noise are small enough to be visible since the number of shots, 4000, is large enough to estimate the observable of the one-local Pauli operator.That is, the difference between solutions obtained by the real device (red circles) and the simulator (blue + symbols) is attributed to the hardware noise and the error mitigation.Nonetheless, Fig. 6 clearly illustrates that the solutions obtained by the real device, ibm kawasaki, captures the result of simulator, which demonstrated that proposed method actually be implementable on a real device.

V. CONCLUSIONS
This paper proposed a scalable quantum circuit implementation method of Hamiltonian simulation for partial differential equations.The key technique for efficiently implementing the time evolution operators on a circuit is to diagonalize each term of Hamiltonian using the Bell basis.We provided concrete quantum circuit representation for Hamiltonian simulation driven by differential operators along with the space and time complexity estimation.For demonstrating the validity of the proposed method, we focused on two partial differential equations, namely the advection and wave equations, which were converted into the Schrödinger equation to be simulated by Hamiltonian simulation.In numerical experiments, we confirmed that the solutions obtained by the proposed method agreed well with the solution calculated by the fully classical finite difference method.This means that the advection and wave equations could be simulated as the Schrödinger equation and actually be implemented on quantum circuits under acceptable Trotter errors.
In the present study, we focused on the advection and wave equations, because they are conservative under appropriate BCs and thus can be exactly converted to the Schrödinger equation.However, because most practical systems are non-conservative, it is required to extend the present approach to systems which cannot directly be described as the Schrödinger equation.Possible approaches are, for example, to describe the target systems as an imaginary time evolution of the Schrödinger equation [51,52] or as the open quantum system governed by the Lindblad equation [53], both of which realize the nonunitary time evolution.We would like to conduct our future research toward such direction.Also, we would like to investigate the relationship among BCs, finite difference schemes, and the resulting Hermitian property in our future research.
where H 1 (Ω) is the Sobolev space defined on the domain Ω.The first term of the last line in Eq. (B1) also vanishes under the Dirichlet BC, i.e., For such U, however, only trivial solution u = 0 ∈ U satisfies the advection equation.

Operator of the wave equation
Let ψ(t, x) and ψ(t, x) be vector fields defined in Sec.II B 2. Applying the integration by parts, we obtain the following relationship about the inner product: where * is the complex conjugate and n is the normal vector pointing outward at the boundary.If the first term in the last line is zero, the operator H in Sec.II B 2 for the wave equation is self-adjoint.This term vanishes when the domain Ω is a unit cell and the periodic BC is imposed on the boundaries, i.e., u ∈ U where where H 2 (Ω) is the Sobolev space defined on the domain Ω.The condition for the self-adjointness of the operator H is also satisfied under the mixed BCs of the Dirichlet and Neumann BCs for u, i.e., when u ∈ U where For j > j ′ = 1, the commutator is rearranged as Thus, the terms of the Hamiltonian, γ(e iλ s − j + e −iλ s + j ), can be grouped into those for j > 1 and j = 1.Since the unitary does not change the operator norm, the Trotter error is upper bounded using the result of Ref. [37,Proposition 9], as follows: Next, we recall Lemma 4.
The terms for the periodic BC commute with the other terms for j > 1 in the Hamiltonian, as follows: and its operator norm is one, which is verified by similar calculation to that of the proof of Lemma 1.Therefore, the approximation error of the Trotter decomposition is upper bounded as We can also obtain the second-order formula as whose approximation error is upper bounded by exp(−iHτ ) − V (2) v 1 τ 2l , − π 2 ≤ v 3 1 τ 3 (2n − 1) 48l 3 , (E9) based on the the similar calculation to that in the proof of Lemma 2. For a d-dimensional case, we can easily extend the discussion here based on the Lemma 4.

Implementation and Trotter error for the wave equation
Here, we consider the following Hamiltonian H for the wave equation in one dimension.The Hamiltonian H is discretized by the forward and backward difference operators so that the resulting H can be actually Hermitian matrix, as follows: which is the extension of the unitary U j (λ) in Eq. (35).Furthermore, we obtain the upper bound of the Trotter error, as follows: which can be obtained by the similar calculation to that in the proof of Lemma 1.We can also obtain the second-order formula as whose error is upper bounded by exp(−iHτ ) − Ṽ (2) cτ l , 0 ≤ c 3 τ 3 (2n − 1) 6l 3 , (E15) based on the the similar calculation to that in the proof of Lemma 2.

FIG. 1 .
FIG.1.Quantum circuit to implement the time evolution operator V , where qj represents the j-th qubit.

FIG. 3 .
FIG. 3. Simulation results of the advection equation with periodic BC in two dimensions.The color plot represents the values of the scalar field |u(t)⟩ calculated by the three approaches.

FIG. 4 .
FIG.4.Simulation results of the wave equation with the mixed BC in one dimension.Solid and dashed lines represent the amplitudes of the quantum state |ψ(t)⟩ prepared by the proposed circuits and that obtained by direct applying the matrix exponential operator e −iHτ to |ψ(0)⟩, respectively.Dotted lines represent the solutions calculated by the finite difference method (FDM).

FIG. 5 .FIG. 6 .
FIG. 5. Simulation results of the wave equation with the periodic BC in two dimensions.Each plot in the column represents the scalar field ∂u/∂t calculated by each approach.

H
= c σ 01 ⊗ D + D − σ 10 ⊗ D − ⊗ s − j + σ 10 ⊗ s + j ) − (σ 01 + σ 10 ) ⊗ I ⊗n   .(E10)By similar calculation to Eq. (37), we obtain exp(−iHτ ) ≈ RX n+1 − 2cτ l n+1 and X n+1 are the RX and X gates acting on the (n + 1)-th qubit, respectively, where X and Y are Pauli matrices, we can naively represent each term of the Hamiltonian (33) by Pauli strings.However, such representation of the Hamiltonian yields the exponentially large number of terms because s − n and s + n are global.Here, we use the Bell basis instead of the Pauli matrices, which can efficiently diagonalize each term of the Hamiltonian (33), as follows: