A Quantum Simulation Approach to Implementing Nuclear Density Functional Theory via Imaginary Time Evolution

The quantum imaginary time evolution (QITE) algorithm is a direct implementation of the classical imaginary time evolution algorithm on quantum computer. We implement the QITE algorithm for the case of nuclear Hartree-Fock equations in a formalism equivalent to nuclear density functional theory. We demonstrate the algorithm in the case of the helium-4 nucleus with a simplified effective interaction of the Skyrme kind and demonstrate that the QITE, as implemented on simulated quantum computer, gives identical results to the classical algorithm.


I. INTRODUCTION
Quantum computers have developed from hypothetical devices to real systems capable of harnessing quantum entanglement and the efficient representation of information to bring a new paradigm to computation and information processing.One major area of application of quantum computers is in the simulation of other quantum systems, where qubit encoding of wave functions of many-particle systems can be enacted particularly efficiently.These ideas have been used in quantum fermionic systems such as quantum chemistry [1][2][3][4], materials [5,6], and nuclear physics [7][8][9][10].In nuclear physics, general quantum algorithms such as the variational quantum eigensolver (VQE) and its extensions [11] have been widely used to find energy eigenvalues in systems such as the deuteron, as represented through an effective field theory [12,13], the Lipkin Meshkov Glick model [14][15][16][17][18][19], and the nuclear shell-model [17,[20][21][22][23].
In the present work, we are concerned with the development and implementation of a quantum computational algorithm of the imaginary time evolution type to solve the nuclear Hartree-Fock equations.We use a simplified but typical nuclear effective interaction of the Skyrme form that gives rise to the usual non-linear Schrödinger equation, which can also be cast as a density functional theory problem.Although Hartree-Fock solutions are achievable on classical computers, implementation on quantum computer is a kind of benchmark [4], and we intend that our algorithm can serve as a starting point for much more sophisticated models with highlycorrelated states, where quantum algorithms should be more useful.The procedure can be considered as a kind of state preparation algorithm for deeper quantum simulation, itself a topic of general interest [24], or as the starting point for quantum algorithms such as symmetry restoration as proposed for use with nuclear-like systems [25,26].We also note that the solution, on quantum computer, of the type of nonlinear Schrödinger equation presented here is of interest beyond nuclear physics [27].
In what follows, we give a brief summary of the nuclear physics problem, a description of the imaginary time evolution algorithm with our particular quantum implementation, and follow with results in the simplest case -a spherical 4 He nucleus in the absence of a Coulomb interaction and where there is a single (s 1/2 ) single-particle wave function to determine.

II. NUCLEAR MODEL
For a nucleus consisting of N strongly interacting nucleons, the time-independent Schrödinger equation governing the N -body wavefunction Φ({x i }) is where x = {⃗ r, S, I 3 } includes all of space, spin, and isospin coordinates.This equation cannot be solved exactly in general and reasonable approximations are usually required.One common method for describing the nuclear structure and low energy dynamics [28,29] is the mean field approximation.By treating each of the N interacting particles as a single particle in the field created by the remaining N − 1 particles, the N -body problem is essentially reduced to a self-consistent one-body problem.Accompanying the mean-field approximation is the form of the nuclear interaction, which needs to be appropriate for use in the mean-field limit; either in the form of a renormalized realistic interaction, or a phenomenological one.We use here the latter form, of the Skyrme kind [30,31].
The Skyrme interaction starts with a potential consisting of a two-body term and a three-body term ( The three-body term is assumed to be a zero-range force arXiv:2308.15425v2 [nucl-th] 23 Jan 2024 while the two-body term is assumed to be short ranged, with matrix elements in momentum space given by where P σ , ⃗ σ, and ⃗ k are the spin-exchange operator, Pauli spin matrices, and relative wave vector respectively (denoting v (2) ijk ).In configuration space equation (4) becomes where ⃗ k L,R = 1 2 i ⃗ ∇ i − ⃗ ∇ j are relative wave vector operators, with the subscripts L, R denoting leftmultiplying and right-multiplying respectively.
It can be shown that [31], for Hartree-Fock calculations, in which the nuclear wave function is assumed to be a Slater determinant of single particle states, the three-body term is equivalent to the density-dependent two-body interaction where ρ (⃗ r) is the single-particle density.
Considering the case of a nucleus with A = 2Z (e.g.

4
He, 16 O), the further simplifications of the absence of a Coulomb field, and with a simplified Skyrme force in which the t 1 , t 2 , and W 0 terms are neglected, we have a simplified model often used for exploratory studies [32,33].The potential is then reduced to an effective two-body interaction With a Slater determinant wave function, the expectation value of the Hamiltonian (i.e. the energy) is By minimizing E via the variational principle, while requiring that the single particle (SP) states are normalized, we obtain the Hartree-Fock equations [32] − (10) The Lagrange multipliers ε j can be identified as the SP energies.
For spherically-symmetric nuclei in the absence of spin-orbit interaction, the SP wave functions φ nlmmsq (⃗ r, s, I 3 ) can be factorised into radial, angular, spin, and isospin parts, in which the angular and spin parts remain uncoupled: The spherical harmonics Y lm (θ, ϕ) form an orthonormal basis and describe the angular behaviour of φ nlmmsq (⃗ r, s, I 3 ), while the spinors χ ms (s) and χ q (I 3 ) represent the spin and isospin dependencies of φ nlmmsq (⃗ r, s, I 3 ) respectively.Equation ( 10) is then reduced to the radial equation where is the Hartree-Fock Hamiltonian.The density ρ (r) is now given by with the factor of 4 in the first line arising from spin and isospin degeneracies.

III. IMAGINARY TIME EVOLUTION
In order to solve the nonlinear Schrödinger equation (10), an iterative method is usually employed.Here, we make use of the imaginary time evolution method.

The time-dependent Schrödinger equation (TDSE) of the Hartree-Fock Hamiltonian Ĥl
has the formal solution Under imaginary time (t → −iτ ), equation ( 16) becomes where N is a normalization operator to renormalize the state after the application of the non-unitary imaginary time evolution operator.When τ → ∞, |ψ (r, τ )⟩ converges to |u 0l ⟩, the ground state of Ĥl HF provided that the initial state |ψ (r, 0)⟩ is not orthogonal to it.[34].
In principle, the imaginary time evolution (ITE) operator Û (τ ) = exp − Ĥl HF τ /ℏ could be applied once, with a large enough imaginary time τ , for a good approximation of the ground state.However, the true form of Ĥl HF is unknown due to its density dependence, and the ITE has to be separated into k total steps, each with an imaginary time step of ∆τ = τ k total .This is equivalent to writing equation (17) as where Û is updated after each step using a newly obtained |ψ (r)⟩.Providing the imaginary time step ∆τ ≪ 1 is small enough, Û can be approximated by

A. Quantum Imaginary Time Algorithm
There have been multiple attempts at implementing ITE on quantum devices [35][36][37][38][39].One of the main obstacles is the non-unitarity of the operator Û as quantum circuits can only handle unitary operators.In the original QITE [40] algorithm, this is overcome by sectioning the qubit chain and performing approximated Trotter evolution [41][42][43] section by section.A later proposal [44] makes use of a randomized qDRIFT algorithm [45] which significantly reduces the circuit depth.In this paper, we use the idea of a duality computer [46,47] to implement non-unitary Hermitian gates, with the aid of some auxiliary qubits.
In order to encode the unknown wave function solution to the HF equation, we have chosen the 3D isotropic oscillator basis (with oscillator length 1 b ) [48][49][50][51] as our computational basis, where α nln ′ are the expansion coefficients and R b,l+ 3 2 n ′ (r) are the oscillator radial wavefunctions.In this basis, the matrix elements of the density are given by the sum of integrals of the product of four basis wavefunctions, These integrals are computed and tabulated classically at the beginning of the calculation, for evaluation of the density at each step of ITE, and a matrix representation of the time-evolution operator ( 19) is constructed.
Using N qubits, the first n ′ = 2 N coefficients of expansion α 0ln ′ can be represented as a state vector |ψ l ⟩ such that where The ITE operator Û , previously obtained classically, is a 2 N ×2 N Hermitian matrix and can be decomposed into a sum of products of Pauli matrices (and identity matrices) acting on individual qubits [52] Û = where i ′ 4 [q] is the q th digit from the right of i ′ when expressed in quaternary and are the identity and Pauli matrices acting on the (q − 1) th qubit.Using i ′ = 45 as an example, since 45 10 = 231 4 , the corresponding gate is The coefficients β i ′ can then be stored in an ancillary state |ϕ a ⟩ using 2N ancillary qubits such that where B = In the larger Hilbert space spanned by |ϕ a ⟩and |ψ l ⟩ the outcome of the non-unitary Û can be obtained in a particular subspace.This is achieved by applying a series of controlled Pauli gates ÔP = and 2N Hadamard gates The result of the operations is given by Thus we obtain our expected value and result state in the first 2 N entries of the statevector.

B. N = 1 case
The details of the N = 1 case help clarify the general algorithm.We thus write out in full the N = 1 case using three (one target and two ancillary) qubits, in which the ground state |ψ l ⟩ is expanded in the two lowest oscillator states.
The ITE operator, Û , is then a 2×2 Hermitian matrix and its Pauli decomposition is Comparing with equation ( 23), we see that and Hence, For a general target state the expected output is Applying Ô from equation ( 28), we obtain Combining the contribution using the Hadamard gates on the ancillary qubits as described in equation ( 30), the qubits will be in the final state where the coefficients of the states |000⟩ and |001⟩ returns the required results from equation (36), up to a normalization constant 1 2B .

C. Quantum Circuits
For a state |ψ l ⟩ expanded in the first 2 N basis states, the QITE algorithm can be performed by a 3N qubit circuit.The N = 1 and N = 2 quantum circuits, which perform one iteration of imaginary time evolution, are For N = 2, the subcircuit SP (2) is (QC5) , with the angles given by tan In general the statevector is not accessible through measurements as only the relative amplitudes can be obtained through the probabilties.In our algorithm we exploit the fact that all of our components are real and the relative phases can either be 0 or π.After direct measurements for deducing the relative amplitudes of the states, N more circuits are composed by appending a Hadamard gate to the N th target qubit.This mixes the amplitudes.An increase in probablity indicates a relative phase of 0 between the two states, and a decrease indicates a relative phase of π.

A. N = 1 case
We applied the approach above on a 3 qubit simulator (N = 1, 2 expansion states) to calculate the ground state energy of 4 He [53,54], where all the nucleons are in the s state (l = 0) with spin and isospin degeneracy = 4, using the parameter values t 0 = −1090.0MeV fm 3 and t 3 = 17288.0MeV fm 6 [55].
We choose a time step of ∆τ ℏ = 0.005 MeV −1 and use an optimized oscillator length of 1 b = 1.5284 fm.Then we perform an ITE, using a classical algorithm, from an initial trial state with equal amplitudes of the first two oscillator states, The wave function starts achieving self-consistency (up to 3 significant figures) after 20 iterations.Further evolution (40 iterations) of the state gives a 4 He ground state energy of −29.1567MeV.The precision of the energy here indicates the self-consistency of the value over the last 3 iterations.In this case further iteration only changes digits from the fifth decimal place.The quantum imaginary time evolution procedure was coded in Qiskit [56] and simulated on a classical computer using the QASM backend with 10000 shots per measurement.The same starting wave function, oscillator length, and number of time steps were used as in the classical state for comparison.Fig. 2 shows the evolution of the wave function and potential as a function of iteration for this implementation of the quantum imaginary time evolution algorithm.The final ground state energy obtained was −29.14 MeV.Selfconsistency is only achieved up to 2 decimal places after 40 iterations.
Fig. 3 shows the final ground state wave functions after 40 iterations from classical and quantum ITE.The per-iteration results and the final energy are thus in close agreement with the classical implementation, as expected.

B. N = 2 case
We repeated the above procedures on a 6 qubit simulator (N = 2, 4 expansion states), without changing any parameters apart from the oscillator lengths and the initial trial state.The initial state is chosen to have equal amplitudes of the first four oscillator states, while the oscillator lengths are optimized to get the lowest ground state energies.The values we use is Compared to the N = 1 case, more iterations are needed for the wave function to achieve selfconsistency.Under the classical algorithm, the wave  the precision of the ground state energy is lower.Improvement of this can be done by using more number of shots per measurement and evolving the states for more iterations.

C. Pauli Coefficients
To gain some understanding of the operation of the quantum imaginary time algorithm, we examine the coefficients, β i of the Pauli operators in the expansion of the unitary operator Û (23).The results are shown in Fig. 6 and 7.  We see the changing over iteration number of the structure of the imaginary time operator as the effective density-dependent Hamiltonian changes, until the operator converges to a fixed form.The relative importance of different terms may allow for a future optimised implementation which eliminates the smallest contributions, but we have not attempted that here.

D. Quantum Advantage
The present work uses a small basis expansion as a proof of concept, where the quantum advantage is not apparent.In this section we compare the computational cost in two different aspects, the number of bits/qubits used (circuit width) and the number of operations/gates used.
With 2 N expansion states, the classical algorithm requires 2 N bits to set up, whereas its quantum counterpart uses only 3N qubits for the presented implementation of QITE.This improvement is exponential and allows one to carry out calculations with a much larger basis expansion size in the eventual availability of larger quantum computers.In our case, we have used a technique in which general non-unitary operators can be used.In theory only one ancillary qubit is needed to implement the circuit.In our algorithm the extra ancillary qubits are traded for a reduced circuit depth.Specialisation to the QITE operator can reduce the ancillary qubit cost further.
In terms of the size of the algorithm, the classical ITE algorithm requires O 2 2N operations per time step.One common measure of the complexity of a quantum ciruit is the number of multi-qubit gates in it.As all quantum gates can be decomposed into 1qubit gates and the 2-qubit CNOT gates, this is equivalent to counting the number of CNOT gates.
Our quantum circuit (QC1, QC2, QC6, QC7) for carrying out an imaginary time step can be decomposed into 4 parts, the state preparation sub-circuits SP (N ) and SP (2N ), the ITE operator ÔP , the contribution combiner ÔH , and the measurement subcircuits.
ÔP consists of 3N doubly-controlled Pauli gates (CCX, CCY, and CCZ), each of which contains five CNOT gates in their respective decomposition [57].Hence ÔP contributes O (N ) CNOT gates to the circuit.ÔH and measurement sub-circuits consist only of the 1-qubit Hadamard gate H and do not add any CNOT gate to the circuit.
The state preparation sub-circuit SP (2N ) is where the bottleneck of the algorithm lies.In our algorithm, we adopt a simple but inefficient way to perform quantum state preparation.In this implementation, the number of CNOT gates scales as O 2 2N .This is acceptable in the small basis expansion we showed.In general this amount of CNOT gates is required, as it corresponds to the 2 2N +1 − 2 independent coefficients of the 2N -qubit statevector [58].There has been attempts to prepare states more efficiently via different approaches [59][24] [60].A recent work using tensor networks proved to achieve linear efficiency for up to 250 qubits [61].With improved state preparation algorithms, our work will prove to be an exponential improvement of its classical counterpart.

V. CONCLUSION
We have presented an implementation of the quantum imaginary time evolution method to solve the Hartree-Fock equations in the case of the helium-4 nucleus.The method is demonstrated to be equivalent to the classical imaginary time evolution algorithm, with a resource scaling superior to the classical case.
The Hartree-Fock technique finds approximate solutions of the nuclear many-body problem that are nearly uncorrelated (except for Pauli exclusion), while much of the promise of quantum computing lies in its ability to simulate highly-correlated states.Hence, the utility of the quantum computation method here is limited to the scaling of the Hilbert space size with the number of qubits.On the other hand, the Hartree-Fock solution is the standard starting point for building correlated wave functions through e.g.densitymatrix methods [62], and our algorithm can thus serve as the starting point for studies in which the ability of quantum computers to efficiently access highly correlated states is exploited.Work along this line is under investigation.
circles and filled circles indicate control activation if the control qubit is in the state |0⟩ and |1⟩ respectively.The subcircuit SP (N) for a Nqubit real state preparation is .a N-qubit real state there are 2 N coefficients (2 N − 1 independent ones), denoted {c i }, i = 0, 1, . . ., 2 N − 1.For the target state |ψ l ⟩ preparation subcircuits SP (N ), c n ′ = α 0ln ′ (as in equation 22).For the ancillary state |ϕ a ⟩ preparation subcircuits SP (2N ), c i = β i (as in equation 23).The angles of rotation, ϑ ij , for the state preparation, are given by Fig 1 shows the evolution of the wave function and the potential as a function of iteration.Even iteration numbers only are shown for clarity.

FIG. 1 :
FIG. 1: Classical imaginary time evolution for N = 1.Upper panel: Wave function as a function of iteration.Lower panel: Potential as a function of iteration.

FIG. 2 :
FIG. 2: Quantum imaginary time evolution for N = 1 carried out using a quantum simulator.Upper panel: Wave function as a function of iteration.Lower panel: Potential as a function of iteration.

FIG. 6 :
FIG. 6: Per-iteration development of the β coefficients for the terms in the Pauli expansion of the unitary operator in the quantum imaginary time evolution algorithm for N = 1 case.Upper panel: Full view.Lower panel: Zoomed in view.

FIG. 7 :
FIG. 7: Per-iteration development of the β coefficients for the terms in the Pauli expansion of the unitary operator in the quantum imaginary time evolution algorithm for N = 2 case.