Quantum Algorithms for Computational Nuclear Physics revisited , particular case of second quantized formulation

No core Full Configurational Interaction (NCFCI) calculations of Nuclear Bonding energy are resource demanding, in particular, computational time scales exponentially with the nucleon number A. In contrast to that, usage of quantum computers would allow an efficient (in polynomial time) NCFCI calculation and speed-up for other beyond-Mean-Field (correlation energy including) methods. To initiate feasibility studies of given quantum algorithms, we present an introduction to preliminary classicalcomputer simulation for the case of spherical nuclei (and He in particular) within NCFCI with realistic chiral NNLO_opt potential.


,
(2) the creation/annihilation refers to w complete orthonormal (ON) sets of one-particle functions (defined on spatial and spin-1/2 domains) -for an ordinary nuclei w = 2, one set corresponding to protons, one to neutrons (e.g.we can reserve negative integers for the former and positive for latter).In the standard formulation, number of creation and annihilation operators for the given type of nucleon must be same in each term in (1).Creation and annihilation operators, as connected to two distinct types of fermions has to obey anticommutation relations (5) and commute with operators for the other nucleonic type

ˆˆˆˆ( ) ,
sign i j a i a j a i a j sign i j a j a i i j     () ˆ,0 While particular example tested consists of at most 2-body interactions (k = 2, the simplest non-trivial), the quantum algorithm presented can be generalized for any k (and any w).In (4) the complex coefficient is therefore identified with matrix element of n-body force operator in one-particle ON basis set.Due to the indistinguishability of identical nucleons and hermitian character of hamiltonian, coefficients V (n) must obey relations where P is an arbitrary permutation of identical nucleon indices.For presented case of k = 2 and w = 2 ON set of Hartree-Fock (HF, Mean-field, [9][10][11]) solutions for the ground state of double magic nuclei has been used, starting with 4 He.The HF equations has following form where Fockians ("(p)" for protons, "(n)" for neutrons ) in arbitrary ON basis sets   where V (1) = h (r) represents operator of kinetic energy, second line on right side ( 10) is an antisymmetrized matrix element of proton-proton (for r = p) or neutronneutron (r = n) interaction a operator and third line represent proton-neutron interaction energy operator matrix element (s = n if r = p and vice versa), the summations is over occupied (0s 1/2) one-particle states for given nucleon.Eq. ( 8) allows (1) to be expressed without V (1) = h (r) matrix elements.There is also possibility to express (1) by (slightly transformed) V (2) elements only (see Supplementary Information).The connection between notations in (1-7) and (8-10) is where  = -1 for r = p and  = +1 for r = n, for the HF ground state wave-function where for 4 He, kZ = k2 = (-1,-2), qA-Z = q2 = (1,2) for the two lowest lying m-levels for each nucleon type.For any realistic nucleon-nucleon potential there is no analytical solution for either HF eigen-problem or an general beyond-HF eigen-problem of the hamiltonian (1).The numerical solution is done by restriction to finite-dimensional one-particle subspaces (defined by index sets I (p) , I (n) ).For most bounded state properties HF solution should serve as a mere ON one-particle basis set generating procedure and more realistic A-body wavefunctions (respecting correlations between particle motions) should be obtained through diagonalization of (1) within (subspace) of product space of restricted oneparticle state spaces for each particle.Due to the Ritz a Please see Supplementary Information for the details.The two-body interaction contains also (second quantized form of) mass-polarization term (1/M)pi•pj (pi being linear momentum of i-th particle) due to the CM S energy subtraction.variational principle [12,13], eigenvalues obtained are upper bounds to eigenvalues of (1) with respect to nonrestricted, complete space.
In case of (1) being diagonalized within all possible (appropriately anti-symmetrized) products of oneparticle functions (Configurations) term Full Configuration Interaction (FCI) is used [11,14,15].In case only certain configurations are used, the method is called Limited Configuration Interaction and possible variants could be CIS, CISD, CISDT, …, for set of at most single, double, triple, … excitations from ground state (or from set of model configurations in MRCI variant) [16].FCI is size-consistent and therefore preferable over limited variants (even when calculated with respect to a greater one-particle sets), only FCI energy will converge with one-particle basis (limiting to completeness) to exact eigenvalue [16]  While the first two classes in fact avoid to calculate FCI energy itself, the only true solution is the last one.Proof-of-principle of a particular algorithm for a theoretical quantum computer device will be presented in the next section.

While the first applications of Abrams -Lloyd algorithm [4] have been tailored for computational (ab initio)
quantum chemistry [4-8] (e.g.FCI for electronic structure for fixed nuclei [26,38], FCI for electrons and nuclei in molecule treated on equal footing via the NOMO approach [30]), the possibility to use Abrams -Lloyd algorithm for computational physics problem has been also previously suggested (e.g.[31]).
Here the main steps of the algorithm will be briefly reminded.The algorithm aim is to estimate hamiltonian H eigenvalue E for given multiparticle stationary state based on representation of evolution operator U = exp(-i t H).

Inputs for algorithm
The inputs for algorithm consists of:

Bounds Emin, Emax to estimated eigenvalue, such t h at
Emin < E < Emax and desired accuracy in number of bits, m (the usual choices would be m = 17 to 27).

Estimate (initial guess)
for any combination of  and m, h has a lower bound, therefore the condition S > Scrit is necessary to assure that pm > 0.5.By eventual repetition of quantum computer algorithm, the success probability can be then amplified as close to 1.00 as desired.
When the algorithm is just simulated for smaller number of one-particle basis states on a classical computer, the desired result and corresponding probability is well known and number of algorithm repetitions r to achieve probabilities of, e.g., 0.99, 0.9999 or 0.999999 can be computed c [30].However, in b Even though this expression is a lower bound to probability, in most cases, the true probability isn't much higher than bound.The possible difference is due to the "overlap" of probability distribution peaks associated to different eigenvalues of studied hamiltonian H.In a typical case, the gap between eigenvalue in question and its neighbors is by orders larger than read-out accuracy (2 -m ) and (13) can be considered as an equation.Possible initial guess should have to be representable efficiently, i.e. through at most polynomially (in oneparticle basis set size) many unitary gates.The most straightforward choice would be Hartree-Fock (Mean Field) guess.However, for a multireference states (common for excited states but in case of s mall gap between highest occupied one-particle state and lowest unoccupied one or for open-shell systems possible for ground state as well) this choice might lead to S < Scrit and more sophisticated initial guess would be recommended.First way would be beyond-Mean field eigenvector calculated with less computationally demanding method than FCI -e.g.CISD, CISDTQ, CCSD, MBPT or FCI [11,16] with excitations restricted to given, limited number of unoccupied (virtual) spin orbitals (Fockian eigenvectors) within same one-particle basis set generated by Hartree-Fock calculation.Second, different one-particle basis set (but connected with the one used in quantum algorithm via finite-dimensional unitary transformation) might be generated from MCSCF or as Natural Orbitals (NO) through diagonalization of one-particle density matrix from Many Body Perturbation Theory (MBPT) of some higher order [11,16].Third, purely quantum approach would be the use of Adiabatic State Preparation (ASP) method [35], [4], which facilitates transition from ground state (represented in quantum register) of one hamiltonian (here Hartree-Fock hamiltonian) to ground state (represented in quantum register) of another hamiltonian (e.g.FCI hamiltonian (1)) under preposition that smooth transformation of hamiltonian in question do exist and that ground state and lowest excited state will neither exchange nor degenerate inside the process (the lowest gap between them then limits number of steps needed).
where overall success probability from majority vote on all r independent algorithm runs is 1-,  = 10 -d (a typical choice for d = 2, 4 or 6 [30]).Wider spread of a register readouts allows for r smaller than upper bound on right side of (16) [30].

Abrams
is applied on each of m qubits of a register (boxes with H symbol in Fig. 1) creating equal weights superposition of binary numbers in a.
Then sequence of conditioned applications of quantum gate representation of U on register b is done (the lowest written wire in Fig. 1), with control qubit being from a register (note the corresponding full black circles in the middle of Fig. 1).This creates entangled state ((2) in [31]) where all x-th powers of U are entangled to corresponding binary representations of x in a register (19, first line).
The pre-last step of PEA is the inverse Quantum Fourier Transform (QFT, represented by a box with ˆI QFT U in Fig. 1, the best-known versions has computational time cost of O(m log m) [39]) which in fact transform state vector in concatenated a+b quantum register (denoted |core> as in [31]) from time domain (time is discretely sampled in the a register) to the energy domain, where All sources of quantum computer powersuperpositioning (Hadamard gate application), parallelism (ease to create all 2 m components in ( 19) in m steps), entanglement (the entangled state in ( 19)) and destructive interference (the inverse QFT) are exploited in PEA and this would lead to an exponential speed up for FCI energy calculation when compared to classical computer-based procedure… under an important preposition that both initial eigenvector guess and evolution operator U application (in general it is a l-qubit gate and should be decomposed into elementary one and two-qubit gates from universal quantum logic gate set [3]) will be represented efficiently (i.e. in number of time-steps (elementary gates) polynomial in one-particle basis set cardinality N).The latter preposition and U gate decomposition will be discussed in following section 2.2.4 and the initial sate preparation in 2.2.5

Computational complexity
The U, exponential of shifted and scaled hamiltonian (1), can be rewritten in a simplified form as where hX are hermitian combinations of terms from expansion (1), e.g. for some X = X0, ( Furthermore, in accordance to the  While the direct mapping is qubit cost uneconomic (there are exponentially many more unphysical components in (25) than actively used in quantum computation!), the gate implementation is very straightforward and computational time is at most polynomial in Nm (if initial state preparation is efficient).

Computational time can be estimated by
where 2 m arise from PEA circuit (Fig. 1), Te = q,

Preliminary classical CI calculations indicate quantum
Abrams-Lloyd algorithm feasible for 4 He ground state FCI calculations, (at least for smaller basis set size) with the simplest single-determinantal initial eigenvector guess (identical to Hartree-Fock solution).For excited states, linear combination of a smaller set of configuration should be used as an initial eigenvector guess.As will be listed in following section, the feasibility studies should carry on for 4 He and later to be extended to larger nuclei.Different quantum algorithms than Abrams-Lloyd should be investigated as well.

3 .
Elements of second quantized hamiltonian (1) The quantum algorithm is probabilistic and provides any of the two closest binary approximates Ebin or Ebin+2 -m , such that the correct result E lies in between them (Ebin < E < Ebin+2 -m , or E = Ebin+2 -m , where   (0;1)) with success probability [3,29] pm overlap S has been defined in second point of input list and factor h depends on reminder  and accuracy m (bits),

-Fig. 1 :
Fig. 1: Phase Estimation Algorithm (PEA) quantum logical circuit.Figure adopted from [31] parametrization of unitary matrix U in question as U = exp(-i t (H -Emin)) (with t = 2/(Emax -Emin) and H being FCI hamiltonain matrix) is supposed from now, Ej is exact FCI energy (minus Emin) for (j+1)-th state and cj is the coefficient of expansion of initial guess is a measurement of all m qubits in a register, denoted by " " in Fig. 1 -to obtain m bits f1, f2, …, fm binary representation of energy through (17).Since in previous step, the quantum core is in an entangled state, where, however, amplitude of component entangled to desired binary representation of a correct energy E(0) should dominate, this step will lead to correct result with probability pm given by (13).At the end of the algorithm, the quantum register a is in state corresponding to measurement outcome and eigenvector in question is encoded inside the quantum register b as 1 b  (with the probability pm it is the eigenvector we searched for).As wave-function is not an observable, there is no simple way how to read it from quantum computer qubit-based memory.However, we could exploit it for some subsequent quantum algorithm -e.g.recover as much as possible information through Quantum Tomography [40-41].

, 1 components
 is number of terms in operator sum in(21) and p is projection from set {1,2,…,2} on set {1,2,…,}.For k,p(k) value determination, please see[44],[3].As discussed in[37], the order  should be optimized.Too low order will provide an inaccurate approximation to the evolution operator U, yet too high order will lead to greater accumulated error (in classical computer simulation arising from rounding, in quantum computation from gate application error and decoherence).The optimization can be formulated (for an ideal quantum computer) as a minimization of q product (since the total number of exponential terms used to represent unitary operator U by elementary quantum gates is equal to q and  is a fixed constant, for minimization technical details, see formula (296) in[37] and Fig.7there).Then, each creation and annihilation operator in (22) or in any hX term is represented via quantum gates acting on register a storing the eigenvector of U. In case of the direct mapping the b quantum register qubits simply store the occupation numbers of m-levels (one-particle state characterized in Harmonic Oscillator scheme by fixed n, l, j and m quantum numbers) and since nucleons are fermions the binary character of qubits is well fit for this representation (occupation numbers are in Let us denote the total number of m-levels used in suggested calculation as Nm.As an example, wave-function representation for 4 He nuclei is is occupation number for j-th m-level for protons, nj n for neutrons, the sum runs over 2^(2Nm)non-zero ci, or are further relevant (since we suppose no-sea approximation, j nj p = j nj n = 2).The size of b quantum register l = 2Nm.In the direct mapping case, creation and annihilation operators are represented by strings of Pauli matrices (Jordan-Wigner operator algebra mapping [45] employs O(Nm) gates, while more economic Bravyi-Kitaev mapping [46] needs just O(log Nm), please see section 2.5.2.1 in [31]) acting on individual qubits of b register.
for number of terms in Trotterization of U = exp(i t H),  stands for number of terms in second quantized hamiltonian (1), (21), Np for the number of them to be executable in parallel and N1 for the number of elementary gates needed to represent exponential of one single term hX in (21).For general case of at most k-

23) the expansion should be partitioned into largest commuting hermitian blocks to minimize number of terms for the next step. In this step, Trotter-Suzuki formula [3,37,42,43] is applied, with general -th order form
, (

the quantum algorithm should be quicker than the classical one. However, since classical FCI diagonalization is done in j-scheme and FCI hamiltonian sparsity is exploited, the denominator in (29) should decrease to 2 or 1.5, leading to critical A in interval (3; 8), where upper bounds are estimated with k = 3 and  = 7. But it can be concluded that quantum FCI will be beneficial, if not from 4 He, from 16 O with certainty. 3 Numerical tests In so far only classical FCI and limited CI preliminary calculations has been done for 4 He, In our calculations we use the nucleon-nucleon optimized chiral potential NNLO_opt. The coupling constants of this interaction were determined by a new optimization method in order to minimize the effects of the three-nucleon force [50], the results are summarized in the following table,Table 1 . Selected 4 He ground state (0 + ) CI energies (ECI) for NNLO_opt potential with  = 16 HO basis set parameter, Hartree-Fock (EHF) and second order Many Body Perturbation Theory (EMBPT), N denote maximal major shell number e and |c0| 2 = S the overlap between CI wavefunction and Hartree-Fock ground state
d Note

that Nm might not be considered as independent to A. At least Nm > A/2 for a non-zero number of virtual one-particle states in FCI calculation. However, the Nm for fixed error in energy might not need to grow more than linearly in A, that is Nm
= A+ for some small positive  and .Number of mshells, Nm, is in relation to number of major shells, n, as Nm ~ O(n 3 ), but number of j-shells as Nj ~ O(n 2 ), therefore Nj ~ O(Nm 2/3 ), further explaining change of denominator in (29).e N = 2n+l, where n starts from zero, e.g.N = 2 means collection of j-shells 0s1/2, 0p3/2, 0p1/2, 1s1/2, 0d5/2