Toward nuclear physics from lattice QCD on quantum computers

One of the ultimate missions of lattice QCD is to simulate atomic nuclei from the first principle of the strong interaction. This is an extremely hard task for the current computational technology, but might be reachable in coming quantum computing era. In this paper, we discuss the computational complexities of classical and quantum simulations of lattice QCD. It is shown that the quantum simulation scales better as a function of a nucleon number and thus will outperform for large nuclei.


Introduction
Quantum computing era is coming around the corner.When quantum computation is put to practical use, paradigm shifts would take place in computational studies of hadron and nuclear physics.One of the most expected fields is lattice quantum chromodynamics (QCD).Lattice QCD researchers have an ambition to study various important physics which suffer from the sign problem in classical computations, e.g., real-time dynamics [1][2][3][4][5][6][7][8][9][10][11][12][13] and properties of dense quark matters [14][15][16].The applicability of quantum computations is not limited to the sign problem, and it is important to open a new avenue to study unsolved problems in lattice QCD.
A nuclear many-body system is one of such challenging problems in lattice QCD.The nuclear many-body system can be, in principle, simulated by lattice QCD because it is just a composite state of quarks and gluons.The simulation is however computationally difficult even with state-ofthe-art supercomputers.As will be discussed later, the computational cost for the corresponding multi-nucleon correlation function is known to blow up quickly as a nucleon number increases.In addition, there exists a so-called signal-to-noise problem, where a correlation function suffers from an exponentially larger statistical error with an increased nucleon number [17], and this issue is related to the sign problem in finite density.Thanks to the significant progress over the decades, the current front-line reaches the simulations of a few-body systems: For instance, two-baryon interactions near the physical point [18][19][20][21][22] and three-nucleon interactions at heavy quark masses [22,23] are obtained from the HAL QCD method [24,25].See also recent calculations of two-baryon interactions at heavy quark masses [26][27][28] from the finite volume method [29].It is, however, still far out of reach of current computations to go beyond few-body systems to determine the properties of large nuclei or many-body forces directly from QCD.
In this paper, we discuss nuclear many-body problem in lattice QCD on quantum computers.
The most important question is whether the quantum computation overcomes the conventional one on classical computers.To answer this question, we estimate their computational complexities in Sec. 2. We also show some results of small-scale noiseless quantum simulation in Sec. 3. Since the Hilbert space of QCD is too large, we employ a toy model of nucleons.Finally, in Sec. 4, we conclude the paper with some future applications.

Computational cost comparison
We compare the computational cost of classical and quantum simulations of lattice QCD.To make the argument concrete, let us consider the simulations for the ground state mass of one nucleus.We estimate how the simulation cost scales as a function of the three-dimensional spatial volume V and the nucleon number A. We treat them as independent parameters, although they are not independent in the sense that the spatial volume must be taken larger than the physical size of the nucleus.We consider quark masses at the physical point and estimate the cost with the lattice spacing a fixed.We implicitly assume the Wilson fermion formalism although it is not essential for the scaling argument below.The Wilson fermion has an advantage in nuclear physics because flavor number is tunable.The general advantages of the Wilson fermion on quantum computing are described in the literature [30].

Classical simulation
The conventional lattice QCD calculation is based on the path integral formalism.The lattice QCD action is given by where ψ and ψ are Grassmann quark fields and U is a gauge link field.The quark fields have internal spinor, color, and flavor spaces and the link field has an internal color space and four directions, but the indices are omitted for notational simplicity.The path integral is given by and the expectation value of the observable O is calculated by They are defined on a four-dimensional lattice V × N τ .
For the calculation of a nucleus with a mass number A, one typically considers a two-point correlation function of nucleus composite fields.The nucleus composite field is defined from A baryon fields as The coefficient C N is a matrix of spin and flavor such that N represents the nucleus of interest.The baryon field is defined from three quark fields as The coefficient C B is a matrix of spinor, color, and isospin such that B represents a color-singlet proton or neutron with spin 1/2.The ground state mass M can be obtained from the temporal correlation In actual calculations, one usually repeats the calculation of correlation functions with various choices for x 1 , • • • , x A (or similar choices in the momentum space) with zero-momentum projection for the total system.As a simplest choice, we consider the spatial summation for the center of mass coordinate of the sink field in the correlation function, and count the corresponding computational cost.Spatial summation usually leads to the implicit improvement of statistical fluctuations with larger V , but we do not consider such effect in the estimate.
The computational cost is estimated as follows.
Quark propagator The quark propagators are computed by conjugate gradient (CG) based algorithms.The convergence of the CG based algorithms is generally dependent on a condition number but asymptotically independent of V .This is because the condition number is given by the difference between the maximum and minimum eigenvalues and they are bounded by lattice cutoff and quark mass, respectively.When the convergence is independent of V , the computational cost for one propagator is proportional to the matrix size V .Gauge configuration The path integral is numerically evaluated by the Hybrid Monte Carlo (HMC) method.As the volume becomes larger, the time step of the evolution of molecular dynamics (MD) should be taken finer, which increases the computational cost with a scaling ) in the case of a second-order MD integrator as leap-frog or Omelyan [31].The scaling can be slightly improved with higher-order MD integrator, but the difference is marginal in the following discussions.Since the computational cost of each time step scales as O(V ) due to the calculation of propagators, the total cost for generating one gauge configuration is 4 ) [31,32].Correlation function The computation of a correlation function involves that of quark propagators and that of contractions.The cost for the propagators with all The contractions consist of two kinds of calculations: Wick contraction and color/spinor contraction.The former corresponds to quark permutation and thus the cost of the straightforward calculation scales factorially with A, t 3a ≃ O(( 32 A!) 2 ).On the other hand, the latter corresponds to the contraction for internal degrees of freedom, whose cost scales t 3b = O(e cA ) [33], where c is a constant which depends on the details of the computational setup and algorithm.We also consider the cost for the zero-momentum projection, t 3c = O(V ).Therefore, the cost of straightforward contraction calculation amounts to be In the last decade, significant algorithmic development to reduce the cost of contraction have been achieved [33][34][35][36][37][38], which are crucial to perform lattice QCD studies for few-baryon systems.If one further considers the scaling with a much larger A, the most significant cost is the factorial factor associated with the permutation, and the algorithm which can reduce this cost to t 3a = O(A 3 ) [39] would be most efficient.Therefore, we take the cost of the contraction as t 3 = O(V A 3 e cA ).Statistical error It is known that statistical error becomes exponentially worse in a larger A system [17].Therefore, it is necessary to increase the number of configurations to keep a relative statistical error O(1) by the scaling of N conf = O(e c ′ A ), where c ′ is a number which linearly depends on the Euclidean time τ in the correlation function.We take τ to be independent of V and A, since the energy gap between the ground state and the first excited state is rather insensitive to V and A, so c ′ becomes constant.
Combining all the estimates above, the total cost scales as in the classical simulation, where 0 < c 1 (≡ c ′ ) < c 2 (≡ c + c ′ ) are some constants which depend on details of computations.

Quantum simulation
The Hamiltonian formalism is favored for quantum simulation.The lattice QCD Hamiltonian is the sum of the quark and gluon parts, The quark Hamiltonian is given by the quark creation and annihilation operators, ψ † , ψ, and the spatial link operator, Û .The gluon Hamiltonian is given by the chromoelectric operator, Ê, and the spatial link operator, Û .These operators can be constructed by sequential manipulations of fundamental gates.The gate costs are O(1) for the gluon operators and O(log V ) for the quark operators [40], but we neglect the logarithmic correction on the following estimation.The Hilbert space is the direct product of the quark and gluon parts.The quark Hilbert space is 2 D quark dimensional with D quark = V (site)×4(spinor)×3(color)×2(flavor) and the gluon Hilbert space is d D gluon dimensional with D gluon = 3V (link).The number of gluon modes d is infinite in principle, but truncated by a finite number in numerical simulation.Although the Hilbert space can be reduced by gauge fixing or solving the Gauss law, we do not consider such a reduction here.
We adopt the quantum adiabatic algorithm [41].The algorithm consists of three steps: state preparation, adiabatic evolution, and measurement.First, we prepare the ground state of a solvable Hamiltonian Ĥ0 for given proton and neutron numbers.The choice for Ĥ0 is not unique.One where m is the quark mass, w is a positive constant, and plaq is the summation over plaquettes.
The potential term v(x) is introduced to resolve the ground state degeneracy and non-existence of a quantum phase transition is assumed between the ground states of Ĥ0 and Ĥ 1 .We define the nucleus operator N from A baryon operators, whereas each baryon operator B(x) is defined from three quark operators, in the same way as Eqs.( 4) and ( 5) except that the argument τ does not exist.The nucleus operator N is constructed such that the ground state of Ĥ0 is given by where |Φ 0 ⟩ is the ground state of Ĥ0 for A = 0, i.e., the vacuum.The ground state of the Hamiltonian Ĥ is obtained by the evolution equation where ĥ(s) interpolates Ĥ0 at s = 0 and Ĥ at s = S.The simplest choice is although the interpolation function is not unique in general.As long as Ĥ0 conserves and down quark numbers, the evolution equation ( 11) does not change proton and neutron numbers.We can obtain the ground state |Ψ A ⟩ with the same proton and neutron numbers as the initial state |Φ A ⟩.
The ground state mass is defined by the energy difference from the vacuum, The computational cost is estimated as follows.
State preparation For the initial Hamiltonian (9), the vacuum |Φ 0 ⟩ is constructed by making all the spatial links be unity and all the antiquarks occupy over the whole lattice.The manipulation of N † is given by making A baryons occupy the sites from lower energy of v(x).The sum of the cost is T 1 = O(V ) + O(A).Adiabatic evolution For the validity of the adiabatic evolution (11), S is large enough that The lower bound is sensitive to the smallest value of the spectral gap ∆(s) between the ground state and the first excited state of ĥ(s).The gap around s = 0 can be controlled and taken arbitrarily large by the choice of Ĥ0 , while the gap around s = S is given by the physical somewhere and the algorithm breaks down.Also, if the evolution passes through a quantum phase transition, the gap goes to exponentially small and the convergence is exponentially slow.Fortunately, lattice QCD has no quantum phase transition between weak coupling limit and strong coupling limit.
energy gap of Ĥ.If ĥ(s) can be optimized, the physical energy gap will be a bottleneck.Two kinds of physical excitation are possible: the internal excitation and the central motion of a nucleus.The internal excitation is independent of V and its energy ∆ int is typically a few MeV.The central motion depends on V and the minimum kinetic energy is given by in an isotropic cube V = L 3 .An example value is GeV, and A = 12.When V and A are large, ∆ cen is smaller than ∆ int .Thus ∆ cen gives the upper bound of asymptotic scaling.The matrix norm || Ĥ − Ĥ0 || of Eq. ( 14) scales as O(V ).
Therefore the lower bound of S is given by O(V ∆ −2 cen ).Since the circuit depth of the evolution ( 11) is proportional to S and one manipulation of ĥ(s) is O(V ), the cost is estimated as If the central motion excitation is absent, the convergence is bounded by ∆ int and the scaling is improved as ).When Ĥ0 is chosen to be translationally invariant like the free Dirac Hamiltonian 2 , the evolution equation (11) conserves the total momentum.We can exclude the central motion excitation by setting the initial total momentum to zero.Since the general construction of translationally-invariant Ĥ0 , other than the free Dirac Hamiltonian, is not known, further studies on this point are desirable in future.Measurement The energy is measured as the expectation value of the Hamiltonian over the volume, so the cost is T 3 = O(V ).Quantum fluctuation One circuit is a sequential execution of the above steps and the circuit is executed N shot times to take the average E = ⟨ Ĥ⟩.The statistical average is accompanied by the statistical error δE = (⟨ Ĥ2 ⟩ − ⟨ Ĥ⟩ 2 )/(N shot − 1).If the theory is located at a quantum critical point, quantum fluctuation diverges, so δE diverges.The circuit must be executed exponentially many times to obtain an order-one relative statistical error, δE/E = O(1), as the volume increases.In nuclear many-body problems, however, we do not encounter a quantum critical point (except for the liquid-gas critical point of nuclear matter).The relative statistical error is insensitive to V and A. Thus N shot = O(1).

VWH[FLWHGVWDWH JURXQGVWDWH
x/a ρ1 (0)ρ 2 (x) Fig. 1 Exact diagonalization results of the energy eigenvalues E (left) and the ground-state correlation function ⟨ρ 1 (0)ρ 2 (x)⟩ (right).The parameters are ma = 1 and λa = 0.5 and the lattice The total cost is given by in the quantum adiabatic algorithm.The quantum scaling ( 17) is at least exponentially better than the classical scaling (7) as a function of A. If the central motion excitation can be excluded, Eq. ( 17) will be improved to

Model demonstration
Unfortunately, we cannot numerically check the asymptotic scaling (17) in lattice QCD.It is far beyond the reach of current technology.Here we adopt a simple toy model, which is a much simplified version of lattice effective field theory [42], to demonstrate the quantum simulation described in the previous section.
Let us consider nonrelativistic fermions on a one-dimensional torus.There are four spin-isospin species, i.e., protons and neutrons with up and down spins.The nuclear force is approximated by contact attractive interaction.The Hamiltonian is with the spin-isospin index i = 1, 2, 3, 4. The Hamiltonian (18) preserves the fermion numbers A i = x ⟨ρ i (x)⟩ and the total fermion number A = i A i .The energy spectrum of each fermion number sector can be exactly computed by matrix diagonalization, as shown in Fig. 1.We considered the four sectors . As the total fermion number increases from A = 2 to A = 4, the energy gap between the ground state and the first excited state decreases and approximately reproduces ∆ = O(A −1 ).At the same time, the difference between the ground state energy and the non-interacting one Am increases due to the attractive interaction.The wave function gets more localized, as can be seen from the two-point correlation function ⟨ρ 1 (0)ρ 2 (x)⟩ in Fig. 1.In the nuclear physics language, up to four nucleons can occupy the first nuclear shell and a helium (A = 4) is the most tightly bound.(Note however that we need more analysis to know whether these are really bound states.)In this model, the ground state mass is equal to the ground state For A ≤ 4 in Fig. 1, the solvable Hamiltonian can be chosen as The corresponding initial state is , where |0⟩ is the empty vacuum.The final Hamiltonian (18) is encoded by the Jordan-Wigner transformation as Because of a periodic boundary condition, a fermionic parity is multiplied to the boundary terms, X i (La) = (−1) A i −1 X i (0) and Y i (La) = (−1) A i −1 Y i (0).The evolution operator is discretized by the second-order Suzuki-Trotter decomposition with a small step δ, and This can be implemented by one-qubit and two-qubit gates, and so no ancilla qubit is necessary.
The number of qubits is D qubit = 4V = 4L.We numerically calculated the evolution equation ( 11) on a noiseless simulator.In Fig. 2, we plot the S-dependence of the ground state mass.The calculation converges at S/a ≃ 15 for A = 2 and at S/a ≃ 40 for A = 4.As explained in Sec.2.2, the convergence is expected to scale as a function of the physical energy gap of the final Hamiltonian (18).From the exact results in Fig. 1, the ratio is {∆(A = 4)/∆(A = 2)} 2 ≃ (0.52/0.87) 2 ≃ 0.4.This is comparable with S(A = 2)/S(A = 4) ≃ 0.4.This demonstrates the expected relation between the convergence and the physical energy gap.The data are averaged over 10 3 circuit executions and the standard deviation is shown as statistical error.The statistical error is insensitive to A as expected.

Summary and far future applications
We have shown that the quantum simulation of lattice QCD has better dependence on the atomic number A than the conventional classical simulation.While the classical simulation suffers from crucial difficulties that the computational cost blows up exponentially as a function of A, it can be avoided by quantum simulations with the Hamiltonian formalism.As for the dependence on the volume V , the quantum scaling (17) can be worse than the classical scaling (7).In practice, however, V and A are not independent.The required size of V will be proportional to A. If the proportionality is assumed, the quantum scaling O(A 16 3 ) is always better than the classical scaling O(e c 2 A A 4 ).
We have assumed fixed lattice spacing and discussed the complexity scaling of A and V because A and V are relevant parameters for nuclear many-body systems.Of course, one eventually needs to take the continuum limit a → 0. In the classical simulation, the computational cost for large A systems on a fixed physical volume grows at least a −4 by the naive counting of fourdimensional lattice points, and the actual scaling could be worse as O(a −4 ) to O(a −7 ) considering additional effects such as a larger autocorrelation in Monte Carlo sampling [32].In the quantum simulation, a similar naive counting of lattice points leads to a −4 scaling.However, the required value of d also increases so that more ultraviolet modes are taken into account.Consequently, the computational cost will rapidly grow [43].This is a crucial issue in the Hamiltonian formalism of lattice QCD and must be resolved to realize quantum advantage.
The scaling argument suggests that quantum computers are supposed to take the place of classical supercomputers at some point in the future, although the practical implementation will take a long time.The required number of qubits is huge; e.g., D qubit = D quark + D gluon log 2 d = 119,537,664 for V = 128 3 , even if the drastic approximation log 2 d = 11 is adopted [44].The problem of quantum error must be resolved, too.If the implementation is achieved after a long journey, the first target would be the ground states of nuclei discussed in this paper, and then low-lying excited states.The calculations provide not only the masses but also the state vectors.
They can be used as the state preparation for the real-time simulation of nuclear reactions, such as nucleus-nucleus scattering and nuclear fusion and fission.Quantum computers could work as emulators of nuclear experiments.

Fig. 2
Fig. 2 Quantum adiabatic calculation for the ground state energy.The broken lines are the exact values and the error bars are statistical errors.The time step is δ = 0.5/a and other conditions are the same as in Fig. 1.