Quantum algorithm for the Vlasov simulation of the large-scale structure formation with massive neutrinos

Investigating the cosmological implication of the fact that neutrino has finite mass is of importance for fundamental physics. In particular, massive neutrino affects the formation of the large-scale structure (LSS) of the universe, and, conversely, observations of the LSS can give constraints on the neutrino mass. Numerical simulations of the LSS formation including massive neutrino along with conventional cold dark matter is thus an important task. For this, calculating the neutrino distribution in the phase space by solving the Vlasov equation is a suitable approach, but it requires solving the PDE in the $(6+1)$-dimensional space and is thus computationally demanding: Configuring $n_\mathrm{gr}$ grid points in each coordinate and $n_t$ time grid points leads to $O(n_\mathrm{gr}^6)$ memory space and $O(n_tn_\mathrm{gr}^6)$ queries to the coefficients in the discretized PDE. We propose a quantum algorithm for this task. Linearizing the Vlasov equation by neglecting the relatively weak self-gravity of the neutrino, we perform the Hamiltonian simulation to produce quantum states that encode the phase space distribution of neutrino. We also propose a way to extract the power spectrum of the neutrino density perturbations as classical data from the quantum state by quantum amplitude estimation with accuracy $\epsilon$ and query complexity of order $\widetilde{O}((n_\mathrm{gr} + n_t)/\epsilon)$. Our method also reduces the space complexity to $O(\mathrm{polylog}(n_\mathrm{gr}/\epsilon))$ in terms of the qubit number, while using quantum random access memories with $O(n_\mathrm{gr}^3)$ entries. As far as we know, this is the first quantum algorithm for the LSS simulation that outputs the quantity of practical interest with guaranteed accuracy.

Investigating the cosmological implication of the fact that neutrino has finite mass is of importance for fundamental physics.In particular, massive neutrino affects the formation of the large-scale structure (LSS) of the universe, and, conversely, observations of the LSS can give constraints on the neutrino mass.Numerical simulations of the LSS formation including massive neutrino along with conventional cold dark matter is thus an important task.For this, calculating the neutrino distribution in the phase space by solving the Vlasov equation is a suitable approach, but it requires solving the PDE in the (6 + 1)-dimensional space and is thus computationally demanding: Configuring ngr grid points in each coordinate and nt time grid points leads to O(n 6 gr ) memory space and O(ntn 6  gr ) queries to the coefficients in the discretized PDE.We propose a quantum algorithm for this task.Linearizing the Vlasov equation by neglecting the relatively weak self-gravity of the neutrino, we perform the Hamiltonian simulation to produce quantum states that encode the phase space distribution of neutrino.We also propose a way to extract the power spectrum of the neutrino density perturbations as classical data from the quantum state by quantum amplitude estimation with accuracy ϵ and query complexity of order O((ngr + nt)/ϵ).Our method also reduces the space complexity to O(polylog(ngr/ϵ)) in terms of the qubit number, while using quantum random access memories with O(n 3 gr ) entries.As far as we know, this is the first quantum algorithm for the LSS simulation that outputs the quantity of practical interest with guaranteed accuracy.

I. INTRODUCTION A. Background
Quantum computing is an emerging technology and has the potential to speedup numerical tasks intractable by classical computers, today's ordinary computers including supercomputers.Witnessing the recent rapid advance of quantum computing, people are now trying to find use cases of quantum computers with the quantum advantage in various fields.In this paper, we focus on the simulation of the large-scale structure (LSS) of the universe with massive neutrino, an important task in cosmology.
In the standard cosmological model, all the rich structures of the present-day Universe formed through gravitational instability of tiny density fluctuations seeded in the early universe [1].The structure at the largest scale probed by cosmological observations is called the LSS.The evolution and the resultant LSS have been shaped by the nature of the mysterious constituents of the universe.Interestingly, the major components of the universe are largely unidentified.In terms of the energy * miyamoto.kouichi.qiqb@osaka-u.ac.jp fraction, about 69% is contributed by the so-called dark energy, perhaps some unknown form of energy different from matter, and about 26% is dark matter (DM), which is often thought to be unknown elementary particles that are not predicted to exist in the Standard Model (SM) of particle physics [2].Observations of the cosmic LSS can shed light on the nature of such dark components and, eventually, provide an important clue on physics beyond the SM.
An important issue toward a better understanding of particle physics through the LSS is that neutrinos have finite masses.Neutrinos are massless particles in the SM, but the detection of neutrino flavor oscillation [3] has now established that they have nonzero mass.The current constraint from the neutrino oscillation experiments is given as the lower bound on the neutrino mass.Although the estimated mass is much smaller than other SM particles such as electron with 0.51MeV, the nonzero neutrino mass is definite evidence that there is physics beyond the SM.Intriguingly, astronomical observations of the cosmic LSS provide independent constraints on the neutrino mass 1 .Neutrinos produced in the early universe exist even today as relics, and constitute a part of matter component, along with cold dark matter (CDM), the conventional picture of DM as particles with much larger mass, e.g., TeV order for supersymmetric particles.In the LSS formation, neutrino behaves differently from CDM and ordinary matter (baryons), which we hereafter call CDM collectively since they behave similarly under the action of mutual gravity.Neutrinos have non-relativistic but extremely large velocities and thus stream almost freely in the gravitational potential of CDM.Hence neutrinos used to be thought as hot dark matter (HDM) distinguishing from CDM.The distribution and gravitational dynamics of neutrinos affect the formation of LSS differently from the conventional picture with CDM only.This gives hope that, by comparing the theoretical prediction on the LSS with massive neutrino to the results of cosmological observations, one can derive constraints on the neutrino mass.We also note that some particle physics models beyond the SM predict other light elementary particles such as axion [8][9][10][11][12], and they can behave as HDM and be constrained from LSS observations too [13].
This background gives us a strong motivation for the numerical simulation of the LSS with massive neutrino.However, it is a computationally demanding task.Nbody simulations are often chosen as a powerful method for the study of LSS, where we employ a large number of particles and follow their gravitational dynamics.Massive neutrinos have been also incorporated into this method, by either adopting approximate corrections or by directly represented by "light" particles [14][15][16][17][18][19][20][21][22][23].However, there remains a concern that such N -body simulations with massive neutrino may lead to imprecise results because, unlike the conventional matter components, CDM and baryons, neutrinos have typically a much larger velocity dispersion, and the so-called shot noise can be significant unless an extremely large number of simulation particles are employed.
An alternative approach is the Vlasov simulation, that is, solving directly the collisionless Boltzmann equation, also known as the Vlasov equation2 : where f (t, x, v) is the neutrino's distribution function in the 6D phase space consisting of the 3D position x and the 3D velocity v at time t, and F(t, x) is the gravitational force per unit mass on a neutrino at position x at time t.However, solving this is also a challenging task.Since there is no exact analytic solution of Eq. ( 1) in realistic settings, we need to resort to some numerical method.A straightforward way to solve a partial differential equation (PDE) is discretization: Setting grid points in the phase space, we can convert the PDE to a linear system of ordinary differential equations (ODEs), and then apply some algorithm for solving ODE.This approach, though, has limited feasibility because of the so-called curse of dimensionality.If we take n gr grid points in one direction, then the total grid number in the 6D phase space becomes n 6 gr and so does the dimension of the ODE.Then, in the time integration of the ODE, the memory space used is of order O(n 6 gr ), and if we take n t time steps in the integration, then O(n 6 gr n t ) queries to the coefficients in the ODE are made in total.This sixthorder scaling makes the Vlasov simulation heavier than the N -body simulation, where we deal with equations of motion in the 3D space.Although there are some studies that try to solve Eq. ( 1) by supercomputers [24][25][26], the room to increase the grid number to improve accuracy is limited.

B. Our contribution
Motivated by the above things, in this paper, we propose the quantum algorithm for the Vlasov simulation of neutrino run on a fault-tolerant quantum computer (FTQC).
(1) into a linear PDE as follows.Since neutrino accounts for a much smaller fraction than CDM, we can neglect the gravity from neutrino to itself and CDM unless the neutrino density is extremely non-uniform.Then, we obtain the gravitational field F CDM (t, x) by CDM using e.g., the N -body simulation in advance, and approximate F(t, x) with F CDM (t, x).The resulting linear PDE can be transformed into the ODE by discretization.Importantly, the ODE has an antisymmetric coefficient matrix A and thus is considered as a Schrödinger equation with the Hermitian iA being the Hamiltonian.We then apply the Hamiltonian simulation, a methodology to generate a quantum state evolved under a given Hamiltonian, and yield the quantum state |f (T )⟩ encoding the value of f on the grid points at the terminal time T in the amplitudes.
As seen later, this takes the O(n gr + n t )3 complexity in terms of the number of queries to the oracle to access the entries in A, which indicates a large speedup from O(n 6 gr n t ).It should be noted that extracting some quantities of interest from the quantum state encoding f in the amplitudes is another issue [63].Then, we also present how to obtain a typical target quantity in the LSS simulation, the power spectrum P ν (k) of the neutrino density perturbation, which indicates the magnitude of the perturbation at the specified scale k, from the quantum state.This is done by some additional unitary operations on |f (T )⟩ and quantum amplitude estimation (QAE) [64].In total, the proposed method outputs an estimate on the power spectrum of accuracy ϵ with O ((n gr + n t )/ϵ) query complexity.
We also note that the proposed quantum algorithm can provide the advantage on space complexity, too.In the above task for calculating P ν (k), our algorithm uses O(log 5/2 (n gr /ϵ)) qubits in total, which is exponentially smaller than the O(n 6 gr ) memory space in the aforementioned classical method.However, our method uses some quantum random access memories (QRAMs) [65] with O(n 3 gr ) entries to store the precomputed CDM gravitational field F CDM (t, x) on the grid points.
In addition, the proposed method can cope with the following complication in the practical problem setting.Reflecting the stochastic nature of the initial value of the perturbation, P ν (k) is defined as an ensemble average.In the classical LSS simulation, it is estimated through multiple runs of the N -body or Vlasov simulation with different initial conditions.In our quantum method, we do not need multiple runs: We can generate the quantum state that encodes the results from different initial conditions in superposition, and estimate P ν (k) with the single quantum state.

C. Comparison to previous studies
Quantum algorithms for solving the Vlasov equation have been considered in previous studies.Most of them give discussions in the context of plasma physics, and there has been no study focusing on the gravitational LSS simulation with massive neutrino as far as we know.Also in the technical aspect, the existing studies are in directions different from ours.
Refs. [66][67][68] presented FTQC algorithms to solve the Vlasov-Poisson equation, in which the force induced by the particles themselves is considered, in the context of plasma physics.They considered the linearized Vlasov-Poisson equation, which describes the perturbative solution on the zeroth-order analytical solution.Compared to this, our approach for linearization, which approximates the force field by only that from CDM based on the small neutrino mass fraction, is different in the following points.First, the CDM gravity is externally given by, e.g., N -body simulation, in which the non-linear dynamics is incorporated, and it is reflected to the neutrino dynamics solved in the current approach.Second, the condition for our approximation to be valid is that neutrino density non-uniformity is not so large that the neutrino self-gravity is negligible compared to the CDM gravity.This is different from the condition for the perturbative approach to be valid, the perturbation being smaller than the zeroth-order solution, which is the neutrino background density in this case.The fractional mass (and energy) density of neutrino is less than one percent of that of CDM, as given in Eq. ( 3) later.Also, recent fully nonlinear simulations show the maximum local over-density of neutrino is of the order unity [24,69].These values depend slightly on the neutrino mass but remain extremely small compared to those of CDM whose maximum over-densities reach 500-1000 in nonlinear "halos."Hence the local gravitational potential is dominated by CDM in most of the regime of interest.
As a difference from the other aspect, although Refs.[66][67][68] also used the Hamiltonian simulation, they worked in the Fourier space instead of working in the position space x like us.This does not fit our setting that the position-wise CDM gravity F CDM (t, x) are given.
There are also studies on quantum algorithms to solve the Vlasov equation in the nonlinear form.Ref. [70] considers approaches via some linearization methods such as Carleman linearization [71].However, like the existing quantum solvers for nonlinear differential equations [50,[56][57][58][59][60][61][62], the method in [70] has some application conditions such as weakness of the nonlinearity.If such conditions are satisfied, then quantum nonlinear differential equation solvers based on Carleman linearization provide a solution with space and query complexities of the same order as linear ones except for some logarithmic overheads [50].Ref. [72] summarizes the prospect of quantum algorithms to solve plasma dynamics in both linearized and non-linear settings, and also mentions the variational quantum algorithms (VQAs) [73].They might be able to run on noisy intermediate-scale quantum (NISQ) devices in the near term, but they are genuinely heuristic.
When it comes to the simulation of self-gravitating systems, in which LSS simulation is included, [74] proposes a VQA to solve the nonlinear governing equation.Ref. [75] also presents a VQA, with the fuzzy DM [76], a specific scenario for DM, in mind.Although their numerical experiment shows a promising result in the proof-of-concept problem, it is unclear whether their methods scale to the larger problem.[77] proposes an FTQC algorithm for the Vlasov equation based on the reservoir method in the context of fluid dynamics, and [78] proposes a similar algorithm for self-gravitating systems.By their method, the Vlasov equation is simulated by appropriately arranging quantum circuits performing increment and decrement operations.While classically controlling the arrangement of the quantum circuits, it allows the delegation of 6-dimensional operations to quantum computation.As a result, their method can reduce the complexity scaling on n gr to O(n 3 gr ), and thus our method achieves larger speedup.
We also comment that the method to extract the power spectrum from the quantum state has not been proposed to our knowledge.

D. Organization
The rest of this paper is organized as follows.In Sec.II, we will explain some basics of the LSS simulation and quantum algorithms used in this paper.Sec.III is the main part, where we will explain each part of our method one by one: discretizing the Vlasov equation, obtaining the solution-encoding quantum state by Hamiltonian simulation, extracting the power spectrum from the quantum state by QAE, coping with the ensemble average, and so on.To illustrate our algorithm, we present a demonstrative numerical experiment on the Hamiltonian simulation-based time evolution of f (t, x, v) in Sec.IV.Sec.V summarizes this paper.

A. Notation
We use I to denote an identity operator.To avoid ambiguity, we may write it as I n with n ∈ N when its size is n × n.We define R + by the set of all positive real numbers.For every n ∈ N, we define [n] 0 := {0, 1, ..., n − 1}.If a matrix A has at most s ∈ N nonzero elements in each row and each column, then we say that A is s-sparse and the sparsity of A is s.
For a vector x ∈ C n , ∥x∥ denotes its Euclidean norm.For an (unnormalized) quantum state |ψ⟩ on a multiqubit system, ∥ |ψ⟩ ∥ denotes the Euclidean norm of its state vector.For a matrix A ∈ C m×n , ∥A∥ denotes its spectral norm, and ∥A∥ max denotes its max norm, the maximum of the absolute values of its entries.
For ϵ ∈ R + , we say that x ′ ∈ R is an ϵ-approximation of x ∈ R, if |x ′ − x| ≤ ϵ holds.We use log and lg for the natural and binary logarithms, respectively.
We label the position of an entry in a vector and the row and column in a matrix with an integer starting from 0. For example, we write a vector v ∈ C n entrywise as v = (v 0 , v 1 , ..., v n−1 ) and call v i the i-th entry, and a matrix A ∈ C m×n as A = (a ij ) i∈[m]0,j∈[n]0 and call a ij the (i, j)-th entry.
We denote by 1 C the indicator function, which takes 1 if the condition C is satisfied and 0 otherwise.Here, we briefly explain some basics of the LSS simulation and review recent developments of simulations that include massive neutrinos.
1. N -body simulation N -body simulations are often employed to simulate the LSS formation.In this approach, the mass distribution is represented by a collection of a large but tractable number of superparticles.By populating the simulation volume with N p superparticles, we numerically solve the time evolution equation for them: where x i (t) is the space coordinate of the i-th super particle and F i is the gravitational force per unit mass on it caused by the other ones.Then, the important statistical quantities such as the power spectrum of the density perturbation are estimated from the distribution of the superparticles.When simulating the LSS that consists of multiple species of particles, e.g., CDM and neutrino, one may naively attempt to use different sets of superparticles for each component.We need to distribute the particles sufficiently densely in the 3D space in order to estimate the quantities of interest accurately.A common setting on N p is N p = n 3 p with large n p , say O(10 2 ), which means N p can be of order of million.Thus, the N -body simulation is a rather heavy calculation, which may need a supercomputer in classical computing.Nevertheless, it is tractable compared to the Vlasov simulation, which deals with the 6D space-velocity phase space.In particular, it is commonly considered that if the velocity dispersion of the particles is small as is the case for CDM, then the N -body simulation yields accurate result.On the other hand, for particles with large velocity dispersion such as massive neutrino, it is difficult to use a sufficiently large number of superparticles so that the distribution in the velocity space is represented, and often the N -body simulation may lead to imprecise results.

Vlasov simulation for massive neutrinos
In light of the above issue, a more desirable approach for LSS simulation with massive neutrinos is solving the Vlasov equation (1) directly.In particular, for neutrino, we can simplify Eq. (1) using the fact that it accounts for only a small part of the matter: The ratio of the neutrino energy fraction Ω ν to Ω m that of the whole of the matter is [24,79] where h is the present value of the Hubble parameter in units of 100km/s and we use Ω m h 2 ≃ 0.14 [2].Based on this fact, we hereafter neglect the self-gravity of neutrino.Then, Eq. ( 1) approximately boils down to the linear equation (5) is the gravitational force per unit mass on a neutrino particle at position x and time t exerted by CDM.This type of approximation can be found also in previous studies [80].We also neglect the gravitational back reaction from neutrino to CDM, and then obtain F CDM from the dynamics of CDM only by, e.g., the N -body simulation.In other words, we consider the neutrino dynamics given the gravity by CDM as an external force field.This leads to the difference between the equation we try to solve and those in the previous studies on the Vlasov-Poisson equation in plasma physics and LSS simulation [66-68, 70, 72, 74].Although both consider the Vlasov equation ( 1) for some particles, in the latter, the force field F is generated by the particles themselves and given through the Poisson equation, and thus the system of the Vlasov and Poisson equations is solved.On the other hand, in our setting, F is explicitly given as a known function F CDM , and only Eq. ( 4) is solved.
Solving Eq. ( 4) numerically is still a difficult task.As explained in Introduction, taking n gr grid points in each of the 6 directions in the phase space leads to the total number of the grid points being n 6 gr .The total grid number rapidly increases with n gr , and so do the query and space complexity in this approach.Although there are some studies in this direction4 that use supercomputers and take hundreds to thousands of grids in each position direction and tens of grids in each velocity direction [24][25][26], increasing the grid point number largely is hard to desire in classical computing.

Neutrino power spectrum
Once we get the neutrino distribution function f (T, x, v) at the terminal time T of the simulation, we can derive the quantities of interest from it as follows.
A quantity we are typically interested in is the neutrino power spectrum P ν (k).It is defined as Here, δν (t, k) := δ ν (t, x)e −ik•x d 3 x is the Fourier component of the neutrino density perturbation at time t, where is the neutrino density and ρν (t) is its average with respect to position.⟨•⟩ denotes the ensemble average with respect to the randomness of the inflationary perturbation, which is assumed to be homogeneous.δ (3) (•) denotes the 3D Dirac's delta function.Because of the assumption that the universe is isotropic, P ν depends on k only through k := ∥k∥ the norm of the wavenumber vector.
Although neutrino and CDM behave as nonrelativistic matters in the present universe, the fluctuations of the neutrino density and the CDM density leave different imprints in cosmological observations.It is implied that the galaxy number density traces the mass density which is mostly contributed by CDM [19,20,[81][82][83][84][85][86].The effect of neutrino is imprinted in LSS over a wide range of length scales [1].Combining the results from, for instance, galaxy surveys and gravitational lensing observations, we can obtain the information on neutrino distribution and on the shape of P ν .Comparing this to the result of the Vlasov simulation, we eventually obtain constraints on the neutrino mass.

C. Quantum algorithms
We now introduce some basics of quantum computing and quantum algorithms used as the building blocks in this paper.For more general basics on quantum computing, we refer to Ref. [87].

Arithmetic circuits
In this paper, we consider computation on the system with multiple quantum registers.We use the fixed-point binary representation for real numbers and, for each x ∈ R, we denote by |x⟩ the computational basis state on a quantum register where the bit string corresponds to the binary representation of x.We assume that every register has a sufficient number of qubits and thus neglect errors by finite-precision representation.For z = x + iy ∈ C, we define |z⟩ := |x⟩ |y⟩, and the real and imaginary parts are held on separate registers.
We can perform arithmetic operations on numbers represented on qubits.For example, we can implement quantum circuits for four basic arithmetic operations: addition For concrete implementations, see [88] and the references therein.In the finite-precision binary representation, these operations are immediately extended to those for real numbers, and then complex numbers.Hereafter, we collectively call these circuits arithmetic circuits.

Block-encoding
Block-encoding, which was first introduced in [36], is a technique to encode a non-unitary matrix as an upperleft block of a unitary matrix.Overcoming the restriction of the unitarity of quantum circuits, this block-encoding technique makes various types of matrix computation efficiently implementable, and thus is widely used in various quantum algorithms as a fundamental building block.We will give the rigorous definition of the block-encoding in Appendix A 1. In some situations, there are some ways to implement the block-encoding unitary.For example, if we have sparse-access to a matrix A, that is, if A is sparse and we have access to oracles that provide positions and values of nonzero entries in A, then we can construct the block-encoding of A, whose implementation cost is shown as Theorem 3 in Appendix A 1. The Vlasov simulation applies to this case, since, as we will see later, A in this case corresponds to the sparse matrix resulting from discretizing the Vlasov equation and we can calculate its entries for given F CDM (t, x).

Hamiltonian simulation
With H a Hamiltonian of a multi-qubit system and an initial state |ψ 0 ⟩, the Hamiltonian simulation is the task to generate the quantum state |ψ t ⟩ after the evolution under H for time t.In other words, we aim to implement the time evolution operator exp(−iHt) as a quantum circuit.This task has been investigated as a core quantum algorithm for a long time [27,28,31], and recently it has been shown that, given a block-encoding of H, we can implement a quantum circuit for the simulation of H with the optimal query complexity [35,[89][90][91].Roughly speaking, we can construct a block-encoding of exp(−iHt) with O(∥H∥t) queries to a block-encoding of H.The rigorous statement is given as Theorem 4 in Appendix A 1.

Quantum amplitude estimation
Even if we obtain some quantum state by a quantum algorithm such as Hamiltonian simulation, extracting information of interest from the state is another issue [63].There is no general prescription unfortunately, and we need to devise a way to extract the quantity we want on a case-by-case basis.In some cases, the necessary quantity is encoded in the quantum state |Φ⟩ as the amplitude of a specific basis state, and then we can use the quantum algorithm called QAE [64] to estimate the amplitude.Roughly speaking, we can estimate the amplitude with accuracy ϵ, using the quantum circuit to generate |Φ⟩ O(1/ϵ) times.The rigorous statement on the query complexity of QAE is given as Theorem 5 in Appendix A 1.

Quantum random access memory
As we will see later, in the algorithm we propose, we need to load the data obtained outside the quantum algorithm onto the quantum register.A QRAM [65] is a device for such an aim.Specifically, we assume that, given a set of N real numbers X = {x 0 , ..., x N −1 }, we can implement the unitary operation on a two-register system like for any α 0 , ..., α Originally, a QRAM was proposed as a dedicated device for data loading, consisting of "atoms" that emit "photons" and enabling the parallel data loading in superposition in time O(log N ) [65].However, some issues on the implementability of such a device have been pointed out, for example, the feasibility of error correction [92,93].On the other hand, there are some proposals on quantum circuit-based QRAM [94][95][96][97][98], into which we can incorporate error correction.However, proposed circuit-based implementations need O(N ) gates and qubits, which reduces the feasibility for large N , although the circuit depth can be logarithmic in N .In any case, implementing a QRAM is highly challenging, and in this paper, we simply assume the availability of the operation (9).

III. QUANTUM ALGORITHM FOR VLASOV SIMULATION OF MASSIVE NEUTRINOS
Now, we present the quantum algorithm for the neutrino Vlasov simulation, including the procedure to estimate the neutrino power spectrum.

A. Discretizing the Vlasov equation
Let us start from discretization, which transforms Eq. ( 4) to a set of ODEs.
We set the range [0, L] for each of the position coordinates x = (x, y, z) and the range [−V, V ] for each of the velocity coordinates v = (u, v, w) with sufficiently large L, V ∈ R + .The region in the phase space we consider is thus As a common boundary condition, we impose the periodic condition for each position coordinate, and the Dirichlet condition f = 0 for each velocity coordinate.Then, we introduce n gr grid points in each coordinate: We set for x, and similarly for y and z, and for u, and similarly for v and w.The grid points in the 6D phase space are We hereafter label the grid points with vectors of indexes: Sometimes, we also use the label i ∈ [N gr ] 0 , which is related to i as Here, N gr := n 6 gr is the total number of the grid points in the phase space.For the later convenience, we assume that n gr = 2 mgr with some m gr ∈ N.
Although we take the same number of grid points in each of position and velocity coordinates just for simplicity, we can take different numbers and such a generalization is straightforward.In fact, in simulations in [24][25][26], the grid number in each position coordinate is different from that in each velocity coordinate.
Using this grid, we approximate the partial derivative with the central difference, that is, where e x = (1, 0, 0), and similarly for the derivatives with respect to y, z, u, v and w.Then, Eq. ( 4) is transformed into f (t) is a vector in R Ngr , which can be seen as . In light of this, we label each entry in f (t) with i ∈ I 6 , as well as i ∈ [N gr ] 0 in Eq. ( 14).
Although we define f (t) as a solution of Eq. ( 16) with some initial value, we expect that the ith entry of f (t) approximates the value of f at time t on the i-th grid point: A(t) is a N gr × N gr matrix and we again label its rows and columns with both i ∈ I 6 and i ∈ [N gr ] 0 .A(t) is written as A x is defined as where D per and E u are n gr × n gr matrices defined as and . . .
respectively.A y and A z are defined similarly.A u (t) is a time-varying N gr × N gr matrix and its (i, j) entry is ; if j = i + e u FCDM,x(t,xi x ,yi y ,zi z ) 2∆v ; if j = i − e u 0 ; otherwise where e u := (0, 0, 0, 1, 0, 0).A v (t) and A w (t) are defined similarly.
We now comment on what are conserved in the time evolution of the ODE system (16).First, ∥f (t)∥ the norm of the solution is constant in time, which is seen as follows.A is antisymmetric, as is easily checked from the definition.Thus, is Hermitian, and we can regard Eq. ( 16) as the Schrödinger equation with the Hamiltonian H.As a property of the Schrödinger equation, ∥f (t)∥ is conserved.This means that in the ODE system (16), there is no instability such that the solution explodes.This is a merit of the current discretization as Eq. ( 15), which leads to the antisymmetricity of A.
Second, f sum (t) := i∈I6 f i (t), the sum of all the entries in f (t) is approximately conserved.This is seen by noting that where (u ↔ v) (resp. (u ↔ w)) denotes the same as the first and second terms expect u is replaced with v (resp.w), and w are the sets of the indexes of the grid points on the boundaries in the velocity coordinates.That is, the time derivative of f sum is almost vanishing except for the contributions from the boundary grids.The conservation of f sum has the physical meaning that the number of the neutrino particles in the 6D simulation box is conserved, which is desirable in the simulation.The violation of the conservation corresponds to the escapes of the particles from the box, which is controllable by setting the box sufficiently large.
Although these conservation properties are desirable, there might be other properties the solution f (t) should have from the physical perspective.For example, since f (t) denotes the distribution density function in the phase space, it must be nonnegative, but, it is not clear that this is satisfied in the current approach.We leave incorporating some schemes guaranteeing the nonnegativity such as upwind difference into the quantum algorithm as future work.

B. Generating the solution-encoding quantum state
As seen above, we can regard Eq. ( 16) as the Schrödinger equation.If we rewrite it a la quantum mechanics, then it becomes where |f (t)⟩ is the quantum state encoding the value of f (t) in the amplitudes: Here, |i⟩ which can be also seen as |i⟩ under the correspondence between i and i in Eq. ( 14), since concatenating the binary representations of i x , ..., i w yields that of i.Now, let us consider generating the quantum state |f (T )⟩ encoding f at the terminal time T by Hamiltonian simulation.To apply this, we make the following assumption.
Assumption 1 (Piecewise time-constancy of F CDM and oracles to access it).Let T ∈ R + and n t ∈ N. Let t it = i t ∆ t for i t ∈ [n t ] 0 with ∆ t := T /n t .Then, for any i holds for any t ∈ [t it , t it+1 ).Here, x ix := (x ix , y iy , z iz ), and, for each CDM,x , F it,ix CDM,y and F it,ix CDM,z are some real numbers.Furthermore, for any i t ∈ [n t ] 0 , we are given accesses to the oracles O it FCDM,x , O it FCDM,y , and O it FCDM,z that act as for any i x ∈ I 3 , where The assumption that F CDM is piecewise constant in time matches the practical setting.From the N -body simulation, we obtain F CDM only on the discrete time points, since in the N -body simulation we discretize the time for time integration of the equation of motion, and interpolating them as Eq. ( 27) is a common approximation.We also note that, for this type of F CDM , the Hamiltonian H(t) is also piecewise constant in time: for any i t ∈ [n t ] 0 , holds for any t ∈ [t it , t it+1 ).
We also assume that we can prepare the quantum state encoding the initial value of f (t).
Assumption 2. For a given initial value f (0), we have an access to the oracle O f (0) that acts as We will discuss the implementation of the above oracles in Sec.III E.
Then, we can apply Hamiltonian simulation to generate |f (T )⟩.
Theorem 1. Suppose that Assumptions 1 and 2 hold.Then, for any ϵ ∈ (0, 1/2), we have a quantum circuit U f (T ),ϵ on the (2 lg N gr + 5)-qubit system that acts as where |f (T )⟩ is the quantum state that encodes the solution of Eq. ( 16) at time T with the initial value f (0) as Eq. ( 26), and |ψ gar ⟩ is some unnormalized quantum state satisfying ∥ |ψ gar ⟩ ∥ ≤ ϵ.In U f (T ),ϵ , (controlled) , arithmetic circuits, and their inverses are queried times, O f (0) is queried once, and qubits are used including ancillary ones.
The proof of this theorem is given in Appendix A 2. Although we postpone the details of construction of U f (T ),ϵ , we now give a rough outline.For the piecewise constant H(t), the solution of Eq. ( 25) is We generate this state via implementing exp (−i∆ t H it ) by block-encoding, which is possible with the oracles O it FCDM,x , etc.We can simplify the query complexity bound (32) by taking into account how to set L and V in practice.To solve the Vlasov equation precisely, we set the considered region V in the phase space so that it contains the region where the distribution function f takes a nonnegligible value.In other words, noting that f reflects the neutrino motion, we set V so that no particle escapes from it.For this purpose, it is sufficient to set L and V , the side lengths of V, as since in the simulation time T a particle traverses at most V T in the position coordinate and at most F max T in the velocity coordinate.Conversely, much larger values of L and V than these are too much.Thus, assuming Eq. ( 35), we can simplify Eq. (32) as as announced in Introduction.Besides, Eq. ( 33) becomes which means that the number of qubits used is only logarithmic.

C. Extracting the power spectra
Even though we can generate the quantum state |f (T )⟩ that encodes f (T ), extracting some information of interest from it is another issue.Here, we consider how to estimate the neutrino power spectrum from the quantum state.In the current setting of finite volume and discrete grid points, the power spectrum is given as [99] Here, is the discrete Fourier transform of the neutrino density perturbation at x ix the grid point in the position space.
is the neutrino density at the position space grid point, calculated by integrating f with respect to the velocity coordinates in the phase space.In Eq. ( 41), (i x , i v ) ∈ I 6 is made by concatenating i x , i v ∈ I 3 .
Leaving how to take the ensemble average in Eq. (38) to Sec.III D, we now consider how to estimate δν Theorem 2. Suppose that Assumptions 1 and 2 hold.Suppose that we are given the value of Let ϵ ∈ (0, 1/2), δ ∈ (0, 1), and i k ∈ I 3 \{(0, 0, 0)}.Then, there exists a quantum algorithm that, with probability at least 1−δ, outputs an ϵ-approximation of δν the solution of Eq. ( 16).In this algorithm, (controlled) , arithmetic circuits, and their inverses are queried times, and qubits are used.Proof.Supposing that we are given |f (T )⟩ = and regarding this as the quantum state on the system of six m gr -qubit registers, we consider the following unitary on this system: Here, Q mgr is the quantum Fourier transform on a m grqubit system, which acts as for any j ∈ [n gr ] 0 , and H mgr is the operation on the same system that is applying the Hadamard gate to each qubit.With m gr = lg n gr , Q mgr is implemented as a circuit consisting of O(log 2 n gr ) Hadamard gates and conditional rotation gates with depth of the same order [87], and H mgr is just a collection of O(log n gr ) Hadamard gates.We thus neglect their costs in the following discussion.We also note that the average density ρν (T ) is related to f (T ) as Using these along with Eqs. ( 39), (40), and (41), we obtain by a straightforward calculation that, for i k := (i kx , i ky , i kz ) ∈ I 3 , where Since the lefthand side is the squared amplitude of the computational basis state |Ω i k ⟩ in W |f (T )⟩, we can estimate this by QAE 5 .
In summary, we can get an ϵ-approximation of δν Input: Accuracy ϵ ∈ (0, 1/2), success probability 1 − δ ∈ (0, 1), the value of C. 1: Construct the unitary quantum circuit U f (T ),Cϵ/4 on a (2 lg Ngr + 5)-qubit system following Theorem 1. 2: Construct W ′ := I32N gr ⊗ W the unitary on the same system, where I32N gr is the identity operator on the first lg Ngr + 5 qubits and W is the operator in Eq. ( 46) that acts on the other lg Ngr qubits.3: Estimate the squared amplitude of The rest of the proof is on the accuracy and complexity of Algorithm 1.We postpone it to Appendix A 3.
To illustrate the outline of Algorithm 1, we present the overview diagram in FIG. 1.Using the oracles O it FCDM,x etc. that give the CDM gravity, we construct blockencodings of the time evolution operators exp(−i∆ t H it ).Combining these block-encodings and the initial-stateencoding oracle O f (0) , we construct the unitary U f (T ),ϵ to generate the state |f (T )⟩ encoding the solution f (T ) at time T with accuracy ϵ.Then, this unitary followed by W generate the quantum state, in which the squared amplitude of a specific computational basis state . Thus, we estimate this squared amplitude by QAE to get δν i k

2
, the squared amplitude of the specified Fourier mode of the neutrino density perturbation.Note that, in this QAE, we use U f (T ),Cϵ/4 , whose accuracy is Cϵ/4, to guarantee the accuracy ϵ for the estimation of δν Let us make some comments on C. In Theorem 2, we assume that we know the value of C in advance, even though it is defined with f (T ) the solution at the terminal time.In fact, this assumption is plausible.This is because, as explained above, ∥f (t)∥ is constant over time, and f sum (t) is also almost constant.After all, we can calculate C using the initial value f (0) instead of f (T ).The cost for this preliminary calculation will be discussed in Sec.III E.
Also, note that C is written as that is, the ratio of the squared average of f i (T ) over the grid points to the average of f i (T ) squared, which we expect is of order 1.By using Eq. ( 35) and C ∼ 1, we simplify Eq. ( 43) to which is the query complexity bound announced in Introduction.Besides, Eq. ( 45) becomes of order (37), the space complexity announced in the Introduction.

D. Extensions
Equipped with Algorithm 1 as a basic one, we now consider its extensions so that it matches complications in practice.

Integrated power spectrum
Although Algorithm 1 outputs an estimate of δν this might not be feasible, since the magnitude of δν i k 2 may be tiny.In fact, in the current discrete setting with n 3 gr grid points in the position coordinates, δ ν (T, x) consists of the n 3 gr Fourier components, and thus δν i k 2 the amplitude of each Fourier components is typically suppressed by a factor 1/n 3 gr .Therefore, to obtain a nonzero estimate of a value of such an order, we need to run QAE with O(1/n 3 gr ) accuracy and O(n 3 gr ) query complexity.After all, the query complexity scaling like Eq. ( 51) is not achieved.
However, this issue is not serious in practice.First, under the usual assumption of isotropy, P ν (k) does not depend on the direction of k but only its norm k = |k|.Furthermore, for the purpose of comparing the numerical simulation with observations, it often suffices to get the integrated power spectrum with some interval [k 1 , k 2 ], which indicates the total magnitude of the neutrino density perturbations of scales between k 1 and k 2 .In the current discrete setting, the corresponding quantity is where Thus, if we want a quantity like this, whose magnitude is larger than the Fourier-component wise amplitude, then the query complexity of our quantum algorithm might not blow up.
Eq. ( 53) can be estimated by Algorithm 1 with a slight modification.Instead of Eq. (A23), we estimate .This probability can be also estimated by QAE with query complexity in Eqs. ( 43) and (44).

Ensemble average
So far, we have considered how to obtain δν i k 2 for one realization, that is, the result of solving Eq. ( 16) with one initial value f (0).However, what we really want is its ensemble average in Eq. (38).In the classical Vlasov simulation, we estimate this ensemble average as that is, the sample average of δν,0 in the n IV different runs of the simulation, where we randomly take the different initial values according to the theoretical distribution of the primordial perturbation [24][25][26].
Also in quantum computing, we may do the same thing, by not running the quantum algorithm many times separately but utilizing quantum superposition.That is, we extend Assumption 2 so that we can generate the superposition of the initial values f (0) (0), ..., f (nIV−1) (0): where n IV = 2 mIV with m IV ∈ N for convenience.Note that the N -body simulation for CDM should be also run many times with the initial values taken randomly, which leads to the different F CDM (t, x) in the different runs.Thus, we also extend Assumption 1 so that we can use the oracle to access F CDM (t, x) in the specified run: where F it,ix,iIV CDM,x is F CDM,x (t it , x ix ) in the i IV th run of the N -body simulation for CDM, and so on.We postpone the implementation of these oracles to Sec.III E.
Equipped with these oracles, we can easily modify Algorithm 1 so that it outputs Eq. (56).That is, we add m IV qubits, and replace O f (0) with O ′ f (0) .We also use Õit FCDM,x instead of O it FCDM,x , and so on.By this change, U f (T ),Cϵ/4 in Algorithm 1 is modified to the operator U ′ f (T ),Cϵ/4 on the (2 lg N gr + m IV + 5)-qubit system that acts as where f (iIV) (T ) is the solution of Eq. ( 16) for the i IV th initial value f (iIV) (0), and |ξ gar ⟩ is some unnormalized quantum state satisfying ∥ |ξ gar ⟩ ∥ ≤ ϵ.Then, we see that the following quantity is equal to the right-hand side of Eq. ( 56), where Eq. ( 60) is the probability that we obtain and can be estimated by QAE.The numbers of queries to the replaced oracles are still of order ( 43) and (44).Note that in the above discussion, we have implicitly assumed that the value of C is the same for different f (iIV) (0).This in fact holds approximately for large N gr , since C is written with the averages of the values of f (0, x, v) and its square on the N gr grid points as Eq. ( 50), and the realizations of f (iIV) (0) for various i IV are generated based on the same distribution of the primordial perturbation.

E. Implementation of the oracles
Now, let us consider the implementations of the oracles used in our quantum algorithm.

Oracles to access FCDM
We first consider , the oracles to access F CDM (t, x) the gravitational force by CDM as Eq. ( 28).Since we now assume that the Nbody simulation for CDM gives the values of this on the n 3 gr grid points in the position coordinates, it seems that we need to resort to QRAM that stores those n 3 gr real numbers.In the recent Vlasov simulations with massive neutrino, the number of the grid points in the position coordinates is of order 10 9 [25].Because of the issues on QRAM mentioned in Sec.II C 5, realizing a QRAM with such a large number of entries seems challenging.Nevertheless, compared to the classical Vlasov simulation, in which the O(n 6 gr ) memory space is used, the QRAM size of O(n 3 gr ) indicates a large improvement with respect to scaling on n gr .
We should also note that calculating {F it,ix CDM,x , F it,ix CDM,y , F it,ix CDM,z } it,ix by the N -body simulation and preparing QRAMs that store them takes O(n t n 3 gr ) time, which has the worse scaling on n gr than the query complexity of our quantum algorithm in Eq. (43).However, as mentioned in Sec.II B, in the current classical computing, the N -body simulation for only CDM is less heavy than the Vlasov and N -body simulations including massive neutrino.Therefore, it is reasonable to consider only the application of quantum computing to the Vlasov simulation for neutrino, leaving the N -body simulation for CDM classical.Also note that, once QRAMs storing F it,ix CDM,x and so on are prepared, we can reuse it for the neutrino Vlasov simulations in the various settings, e.g., various neutrino masses, which relatively diminishes the cost for the N -body simulation for CDM.
We also comment that, the approximation that F CDM (t, x) is piecewise constant in time reduces the size of QRAMs to be prepared.Although we are given {F it,ix CDM,x , F it,ix CDM,y , F it,ix CDM,z } it,ix , whose number per vector element, say F CDM,x , is n t n 3 gr , our algorithm uses not one QRAM with n t n 3 gr entries per vector element, but n t QRAMs with n 3 gr entries, which are {F it,ix CDM,x } ix for one value of i t ∈ [n t ] 0 .This is because, as explained in Sec.III B, to get |f (T )⟩, we successively apply Hamiltonian simulations for the time intervals [t 0 , t 1 ), ..., [t nt−1 , t nt ), and in each of them we use only the constant value of F CDM (t, x) in the interval.Dividing the QRAM into smaller ones reduces the difficulty of implementation.In particular, for the circuitbased QRAM, we do not need to use n t n 3 gr qubits, but it suffices to reuse n 3 gr qubits, with different QRAMs implemented in the Hamiltonian simulation for the different time interval.We do not enjoy this reduction if F CDM (t, x) continuously varies in time and we use the time-dependent Hamiltonian simulation [100] As the last comment, we note that when we want to estimate Eq. ( 56), the average of δν,0 , the number of the entries in each QRAM increases to n 3 gr n IV .

Oracle to generate the state encoding the initial value
Next, we consider the oracle O f (0) in Eq. ( 30), which generates the quantum state |f (0)⟩ encoding the initial value f (0).
In the neutrino Vlasov simulation, the initial distribution function is set to a Fermi-Dirac distribution [26], Here, v b (x) is the neutrino bulk velocity at position x, and where k B is the Boltzmann constant and T ν is the neutrino temperature at the initial time.The initial density fluctuation and bulk velocity are calculated based on the theory for the primordial cosmological perturbation [26].Then, our aim is generating the following quantum state |f (0)⟩ encoding Eq. (61).It is written as where and C and D v b are normalization constants.We now assume that {δ ν (0, x ix )} ix are classically computed as a preparation of our quantum algorithm and stored in a QRAM.Given such a QRAM, we can generate the quantum state C ix∈I3 (1 + δ ν (0, x ix )) |i x ⟩ querying the QRAM O(log n 3 gr ) times by the method in [101].We also assume that {v b (x ix )} ix are precomputed and stored in a QRAM.There are some methods to generate a quantum state encoding a function given as an explicit formula in the amplitudes [102][103][104][105].By such a method, we can perform the operation gr ) times.Combining these circuits and QRAMs, we can generate |f (0)⟩ with O(polylog(n gr )) query complexity.Note that the sizes of the aforementioned QRAMs are O(n 3 gr ), which is of the same order as QRAMs storing F CDM .Also note that preparing those QRAMs takes a O(n 3 gr ) time and this does not increase the order of the preparation time, which already includes the O(n 3 gr ) time to prepare the QRAMs for F CDM .
In the above discussion, we have assumed that there is only one neutrino flavor with mass m ν .This is just for simplicity, and the extension to the more realistic multiflavor setting is straightforward as follows.The point is that we are now approximating the neutrino Vlasov equation as Eq. ( 4) and thus the different masses of the different flavors affect the solution only through the initial value, which depends on the mass.Since Eq. ( 4) is linear with respect to f , if we set the initial value to the linear combination of those for the various flavors, then the solution becomes the distribution function for the mixture of the flavors.In our quantum algorithm, the modification is only replacing |F FD,v b ⟩, which encodes the single-flavor distribution, with the quantum state encoding the initial distribution of the mixture.Since it is still given as an explicit formula, we can use the aforementioned methods for function-loading to the quantum state.
Lastly, we make a comment on the cost to calculate C, which has been postponed.The numerator of Eq. ( 50) is equal to and thus written with ρν (0), which is determined by the normalization of f .On the other hand, the denominator of Eq. ( 50) becomes 2 (66) for f (0, x, v) in the form of Eq. ( 61).The latter part of Eq. ( 66) is calculated as where we replace the sum with the integral at the first approximation, and at the second approximation we neglect the contribution to the integral from the region outside 67) is independent of i x and evaluated numerically.The remaining part is We can calculate this at the same time as generating {δ ν (0, x ix )} ix , keeping the order of the preparation time O(n 3 gr ).

IV. DEMONSTRATION
Lastly, we present an illustrative numerical model.It is impossible to run our algorithm on a current real quantum computer or a quantum circuit simulator.Instead, we demonstrate a core part of our algorithm, namely, the Hamiltonian simulation-based time evolution of the neutrino distribution function f (t, x, v).By considering a toy problem, we find the solution of Eq. ( 25), the discretized Vlasov equation, as Eq. ( 34) with the matrix exponentiation exp (−i∆ t H it ) calculated explicitly on a classical computer.This is what we aim at by the quantum algorithm of Hamiltonian simulation, and differs from common classical methods for time integration.
We consider the 1D setting with the position coordinate x and the velocity coordinate v.We give the CDM gravity by the following explicit function with real constants A and K.Then, under the discretization described in Sec.III A, we find f (T ) = f (T, x 0 , u 0 ), ..., f (T, x ngr−1 , u ngr−1 ) , the vector of the values of f (T, x, u) on the grid points.Since F CDM,x does not depend on time and thus the Hamiltonian H does not either, we can get f (T ) by applying a single operator exp (−iHT ) to f (0): The initial value is set as follows: f (0, x, u) is constant in x, and the Maxwell distribution is adopted in the u direction.That is, with σ v > 0.
We show the result for A = −1, K = π and σ v = 0.1 under the discretization with L = 2, V = 1 and n gr = 64.FIG. 2 shows the heatmaps of f (T, x, u) drawn with f (T ) obtained by Eq. ( 69) with T = 0, 0.1, and 0.2.From this f (T ), we calculate δ ν ix , the neutrino density perturbation at each grid point, as Eq. ( 40), and plot it as a function of x = i x L/n gr in FIG. 3.Because of the functional form of F CDM,x , we expect that the density is enhanced at x = 1 and suppressed at x = 0, 2, and FIG. 3 matches this expectation.Using this δ ν ix , we calculate | δν i k | 2 , the squared amplitude of the Fourier mode of the neutrino density perturbation, via Eq.( 39), and plot it as a function of k = 2πi k /L in FIG. 4. Since the neutrino distribution evolves under the force as a single sinusoidal function with wavenumber k = K, we expect that only the Fourier mode with that wavenumber is enhanced.In fact, this is observed in FIG. 4. In summary, these figures imply that our calculation is working in this demonstration.

V. SUMMARY
In this paper, we considered the quantum algorithm for simulations of LSS formation with massive neutrinos.The large-scale neutrino distribution is an important issue for both cosmology and particle physics.In order to follow the growth of density perturbations of neutrinos with a large velocity dispersion, it is desirable to solve directly the Vlasov equation in an efficient way rather than performing conventional N -body simulations.However, one needs to solve the PDE in a (6+1)-dimensional space, and thus it is a challenging task: Taking n gr grid points in each of six space coordinates and n t time grid points leads to O(n t n 6 gr ) and O(n 6 gr ) space complexity.We thus proposed a quantum algorithm for this task.First, by neglecting the gravity generated by the neutrinos that contribute only a tiny fraction in mass, and hence by taking into account only the gravity by CDM, we approximated the Vlasov equation in a linear form.We showed that the discretized Vlasov equation can be regarded as a Schrödinger equation.Then, by calculating the distribution of CDM using a fast N -body simulation in advance, we applied the Hamiltonian simulation to solve the equation.We proposed not only how to obtain the solution-encoding quantum state but also a way to extract the neutrino power spectrum, which is an important quantity for comparisons with cosmological observations and thus is of practical interest, from the quantum state by QAE.Our method outputs a ϵ-approximation of the power spectrum with O((n gr + n t )/ϵ) query complexity, using O log 5/2 (n gr /ϵ) qubits.Although our algorithm has some shortcomings such as the assumption on the availability of large-sized QRAMs, it is the first proposal of a quantum algorithm for the LSS simulation that outputs the quantity of practical interest with guaranteed accuracy.We thus believe that our algorithm is an important step in the application of quantum computing to the LSS simulation or, more broadly, challenging numerical tasks in cosmology and particle physics.
In future works, we will explore the possibility of extending our algorithm.Obviously, taking into account the neutrino self-gravity via, e.g., Carleman linearization, is one possible way.Another possibility is incorporating higher-order schemes for partial derivatives, for which we are now adopting the central difference as Eq. ( 15).Adopting higher-order schemes has two effects.On the one hand, it makes the coefficient matrix A(t) less sparse: If we adopt a n-th-order scheme, then the sparsity s scales as O(n).In the sparse-access setting, the query complexity of Hamiltonian simulation scales as O(s) on s, and thus so does our algorithm.In total, the query complexity of our algorithm scales as O(n) on n.On the other hand, by adopting a higher-order scheme, we may configure spatial grid points with larger interval, which means smaller n gr , keeping the accuracy of the solution.Since the query complexity of our algorithm scales as O(n gr ), in total, if an n-th-order scheme leads to n gr scaling better than n gr = O(1/n), then higherorder schemes would be beneficial.As a previous study on the effect of higher-order schemes in practice, we refer to [106], which performed a set of simulations with third-, fifth-, and seventh-order schemes and in fact observed the improvement of the accuracy.However, in general, the extent to which higher-order schemes improve accuracy is highly problem-dependent, and thus we leave incorporating higher-order schemes into our algorithm for future works.

c. Quantum amplitude estimation
The complexity of QAE is evaluated as the following theorem.
Theorem 5 ([64], Theorem 12,modified).Suppose that we are given an access to a quantum circuit A that acts on a quantum register as , where |ϕ⟩ and |ϕ ⊥ ⟩ are mutually orthogonal states and a ∈ (0, 1).Also, suppose that we have access to a quantum circuit V on the same system that acts as V |ϕ⟩ = − |ϕ⟩ and V |ϕ ⊥ ⟩ = |ϕ ⊥ ⟩.Then, for any ϵ ∈ R + and δ ∈ (0, 1), there exists a quantum algorithm that with probability at least 1 − δ outputs a ϵ-approximation of a calling A and times.
In this paper, we use QAE for the case that the target state |ϕ⟩ is a superposition of some set S of computational basis states and |ϕ ⊥ ⟩ is a superposition of the other ones.In this case, a is the probability that we get a bit string corresponding to any state in S when we measure |Φ⟩.
Note that, although the success probability in the original theorem in [64] is a constant, it is 1 − δ in Theorem 5 with a factor log(1/δ) added in the query complexity bound.This is due to taking the median of the results in multiple runs of the algorithm [107], which is based on the powering lemma in [108].
We also comment on the qubit number.The original version of QAE proposed in [64] is based on quantum phase estimation (QPE) [109], and uses O(log(1/ϵ)) qubits to output the estimate of a in addition to qubits used in A and V .Furthermore, it requires the controlled versions of A and V , for which we may need to use additional qubits and gates compared to the uncontrolled ones.Fortunately, after [64], many variants without QPE, which require neither additional qubits nor controlled oracles, have been proposed [110][111][112][113]. Thus, we hereafter consider that in QAE we use only qubits needed to operate A and V .
In this operation, the number of queries to (controlled) and that to O f (0) is 1.Then, combining these, we bound the total query number in Algorithm 1 as Eqs.( 43) and (44).
On the qubit number, we note that, to operate U f (T ),Cϵ/4 , we use qubits whose number if of order (45), as implied by Theorem 1.Since operating W and QAE do not require additional qubits, the number of qubits used in the whole of Algorithm 1 is also of order (45).

i k 2 from
|f (T )⟩.Formally, we have the following theorem.

FIG. 1 .
FIG. 1. Overview of the quantum circuit used in Algorithm 1.

i k 2 by the procedure shown in Algorithm 1 . 1
Algorithm Estimation of δν i k 2
with any t ∈ [0, T ), we have O in Eq. (18), then we immediately have those for H(t) in Eq. (23), since the latter is just the former multiplied by i.Thus, we hereafter consider O row .For A(t), whose sparsity is 12, r ik is given as6 cuits are queried O (1) times.Proof.If we have O A(t) row , O A(t) col , and O A(t) ent for A(t) it FCDM,x , O it FCDM,y , O it FCDM,z it, arithmetic circuits and their inverses are that in each V it , which is bounded as Eq.(A12), multiplied by n t , that is, O n t ∥H it ∥ max ∆ t + log Lastly, we prove the statement on the qubit number.Combining Theorems 3 and 4, we see that the number of the ancilla qubits used in V it isO log 5/2 ∥H it ∥ maxPlugging Eq. (A20) into this, we get Eq.(33).Adding the qubits on which the quantum state in Eq. (A17) is generated, whose number is 2lgN gr +5 = O(log n gr ), does not change the order.Continuation of the proof.Let us see that the stated accuracy is achieved.By U f (T ),Cϵ/4 , we get the following stateU f (T ),Cϵ/4 |0⟩ ⊗(2 lg Ngr+5) = |0⟩ ⊗(lg Ngr+5) |f (T )⟩ + |ϕ gar ⟩ , (A22) where |ϕ gar ⟩ is an unnormalized state with ∥ |ϕ gar ⟩ ∥ ≤ Cϵ/4.Thus, definingp := Ω ′ i k W ′ U f (T ),Cϵ/4 |0⟩ = 2ℜ ⟨Ω i k | W |f (T )⟩ Ω ′ i k W ′ |ϕ gar ⟩ + Ω ′ i k W ′ |ϕ gar ⟩Here, at the first inequality, we use| ⟨Ω i k | W |f (T )⟩ | ≤ 1 and ∥ Ω ′ i k W ′ |ϕ gar ⟩ ∥ ≤ ∥(W ′ ) † Ω ′ i k ∥ • ∥ |ϕ gar ⟩ ∥ ≤ Cϵ 4, which follows from the Cauchy-Schwarz inequality.pobtained in Step 3 in Algorithm 1 is an Cϵ/4-approximation of p, and combining this with Eq. (A24) yields Thus, the error of the output of Algorithm 1 is bounded by ϵ.Lastly, let us prove the statement on the query complexity.In the QAE in Step 3 in Algorithm 1, U f (T ),Cϵ/4 is queried