Ab initio nuclear structure via quantum adiabatic algorithm
Weijie Dua,b, James P. Varya, Xingbo Zhaob,c, and Wei Zuob,c

aDepartment of Physics and Astronomy, Iowa State University, Ames, Iowa 50010,
USA
bInstitute of Modern Physics, Chinese Academy of Sciences, Lanzhou 730000,
China
cSchool of Nuclear Science and Technology, University of Chinese Academy of
Sciences, Beijing 100049, China

## Abstract

#### Background:

Solving nuclear many-body problems with an ab initio approach is widely
recognized as a computationally challenging problem. Quantum computers offer a
promising path to address this challenge. There are urgent needs to develop
quantum algorithms for this purpose.

#### Objective:

In this work, we explore the application of the quantum algorithm of adiabatic
state preparation with quantum phase estimation in ab initio nuclear structure
theory. We focus on solving the low-lying spectra (including both the ground
and excited states) of simple nuclear systems.

#### Ideas:

The efficiency of this algorithm is hindered by the emergence of small energy
gaps (level crossings) during the adiabatic evolution. In order to improve the
efficiency, we introduce techniques to avoid level crossings: 1) by suitable
design of the reference Hamiltonian; 2) by insertions of perturbation terms to
modify the adiabatic path.

#### Results:

We illustrate this algorithm by solving the deuteron ground state energy and
the spectrum of the deuteron bounded in a harmonic oscillator trap
implementing the IBM Qiskit quantum simulator. The quantum results agree well
the classical results obtained by matrix diagonalization.

#### Outlook:

With our improvements to the efficiency, this algorithm provides a promising
tool for investigating the low-lying spectra of complex nuclei on future
quantum computers.

##  1 Introduction

Ab initio nuclear structure theories are powerful tools for investigating
nuclear properties [1, 3, 4, 2, 5]. One approach [2] resembles configuration
interaction (CI) methods in atomic physics and quantum chemistry since it
seeks the Hamiltonian eigenvalues and eigenvectors in a chosen basis space.
Such theories formulate the nuclear many-body problem in terms of the strong
inter-nucleon interactions within a Hamiltonian that acts on all the nucleons
in the nucleus. The three-dimensional harmonic oscillator (3DHO) basis is
commonly adopted. Within the basis representation, the Hamiltonian of the
nuclear system, incorporating the nucleon kinetic energies with the inter-
nucleon interactions, generates a sparse matrix of large dimension [4, 2, 6,
7, 8]. One of the main goals is to compute the low-lying spectrum to compare
with and/or predict experimental results.

However, the reach of the CI ab initio method (e.g., the no-core shell model
[2, 9, 10]) is restricted by available computational resources [6, 7, 8]. Due
to the complex nature of the nuclear many-body problem, the dimension of the
model space grows exponentially with increasing number of constituent
nucleons. The resources necessary for obtaining converged results with
realistic inter-nucleon interactions exceed those available on world-leading
supercomputers for all but nuclei of low atomic number, typically less than 20
nucleons. Quantum computers hold the promise of reducing the exponential
scaling requirements to polynomial [11]; and quantum computing techniques
would open up a fruitful path to ab initio solutions for nuclei spanning the
full table of nuclear isotopes.

For these reasons, there is a great need to develop quantum algorithms for
performing nuclear structure calculation on quantum computers. Several
heuristic and ground-breaking algorithms have been proposed in quantum
chemistry and condensed matter physics [12, 13], including the adiabatic state
preparation with quantum phase estimation (ASP-QPE) [14, 15, 16], the hybrid
quantum-classical methods such as the variational quantum eigensolver (VQE)
[17, 18, 19], the quantum approximate optimization algorithm [20], the quantum
imaginary time evolution [22, 21] and the quantum Lanczos methods [21, 23].
Among these, the VQE method has been adopted to treat simple nuclear systems
such as the deuteron [24, 25]. This successful hybrid method is flexible with
respect to circuit depth and can be executed on near-term quantum computers
[26]. However, it is still unclear if this method can be easily adopted to
excited states and to more complex nuclei (for one reason, the VQE method
entails a judicious design of the trial function and additional noisy high-
dimensional classical optimization [27]).

In this work, we explore an alternative idea of using the ASP-QPE algorithm to
calculate the low-lying spectrum of the nuclear system. In principle, if one
starts with the ground state of the reference Hamiltonian and slowly evolves
the state along the adiabatic path that connects to the target Hamiltonian,
then one is guaranteed to obtain the ground state of the target Hamiltonian
according to the well-celebrated adiabatic theorem [28, 29]. By the consequent
QPE algorithm [30, 31, 32], the ground state energy of the target Hamiltonian
can be obtained. The ASP-QPE algorithm can also be generalized to solve the
excited states and their eigenenergies [33], as we will show in this work. In
view of the fact that the computational complexity of the ASP-QPE algorithm
dramatically increases due to the emergence of small energy gap(s) (referred
to as “level crossing” in the text) during the adiabatic evolution, we
investigate two ideas to address this issue (cf. Refs. [16, 34, 35]): 1) to
design suitable reference Hamiltonians; and 2) to introduce additional
perturbation terms. These two methods modify the adiabatic path and avoid the
level crossing such that the efficiency of the ASP procedure as well as the
probability for the QPE to observe the desired eigenenergy are improved.

We illustrate the ASP-QPE algorithm by solving the deuteron ground state
energy. For the illustration of solving excited states, we apply the ASP-QPE
algorithm to investigate the spectrum of the deuteron bounded in an external
harmonic oscillator trap, where this system mimics the neutron-proton pair
that exists in the mean field of the nuclear environment.

We note that the ASP-QPE algorithm takes deep circuits with extra ancilla
qubits, which makes it hard to execute the algorithm in near-term quantum
platforms [26]. However, since the qubitization procedure is simple (we adopt
a compact mapping scheme to map the Hamiltonian eigenvalue problem onto the
quantum computer) and the ASP prepares deterministically the eigenstate of
interest from the input reference state, one can generalize the ASP-QPE to
complicated nuclei directly, without the complexity of designing the ansatz
trial function and solving high-dimensional optimization problems as in the
VQE approach. Also, since the target Hamiltonians are large sparse matrices
for complex nuclei, efficient quantum circuits can be designed for the
simulation [36, 37]: with the level-crossing problem addressed here, the
efficiency of the ASP-QPE algorithm can be improved, and we expect that the
algorithm will provide a promising candidate for ab initio nuclear structure
calculations on future quantum computers.

This paper is organized as follows. We present the theory in Sec. 2, where we
discuss various target Hamiltonians, our designs of the reference Hamiltonians
and the adiabatic path. In Sec. 3, we discuss the ASP-QPE algorithm, where we
also present our ideas to avoid the level crossing. We show the simulation
conditions in Sec. 4 of the model problem, where we also discuss the compact
mapping scheme, the preparation of the reference states, and the design of the
quantum circuit. We present the results and discussions in Sec. 5. We conclude
in Sec. 6, where we also present future plans.

##  2 Theory

Table 1: The summary of the target Hamiltonians, the reference Hamiltonians, and the adiabatic paths in this work. See the text for details. | Natural deuteron | Trapped deuteron
---|---|---
Target Hamiltonian | $H_{\rm target,1}=T_{\rm rel}+V_{\rm NN}$ | $H_{\rm target,2}=T_{\rm rel}+V_{\rm NN}+{\alpha[V_{\rm HO}+V_{\rm depth}]}$
Reference Hamiltonian | $H_{\rm ref,1}=H_{\rm HO}+V^{(0)}_{\rm NN}+V_{\rm shift}$ | $H_{\rm ref,2}=H_{\rm HO}+V^{(0)}_{\rm NN}+\alpha V_{\rm depth}+V_{\rm shift}$
Adiabatic path | $H_{1}(t)=[1-f(t)]H_{\rm ref,1}+f(t)H_{\rm target,1}$ | $H_{2}(t)=[1-f(t)]H_{\rm ref,2}+f(t)H_{\rm target,2}$
Path modification | N.A. | $H^{\prime}_{2}(t)=H_{2}(t)+h_{i,j}(t)[1-h_{i,j}(t)]H^{i,j}_{\rm path}$

In this work, we apply the ASP-QPE algorithm to solve for the ground state
energy of the deuteron. We also generalize the method to solve the low-lying
spectrum of simple nuclear systems: we illustrate this exploration with the
trapped deuteron (deuteron system bounded in an external 3DHO trap). In order
to improve the efficiency of the algorithm, we discuss the ideas for 1)
adopting proper reference Hamiltonian for the ASP algorithm; and 2) modifying
the adiabatic path by additional perturbation terms. In Table 1, we summarize
the Hamiltonians and adiabatic paths applied in this work, where the detailed
discussions are presented in the following subsections.

###  2.1 The natural deuteron

We start from the intrinsic Hamiltonian of the deuteron:

$\displaystyle{H}_{\rm target,1}={T}_{\rm rel}+{V}_{\rm NN},$ (1)

where ${T}_{\rm rel}$ denotes the relative kinetic energy of the neutron-
proton system and $V_{\rm NN}$ is the nuclear interaction. In general, $V_{\rm
NN}$ can be taken from various model-dependent nuclear interactions (see,
e.g., Refs. [38, 39]), or those derived from first principles, e.g., chiral
effective field theory (see, e.g., Refs. [40, 41, 42]).

We follow Ref. [2, 9, 10, 43] and use the three-dimensional harmonic
oscillator (3DHO) basis throughout this work. We set the oscillator energy
$\hbar\omega=5$ MeV. For the deuteron, we take the orbital angular momenta
$l=0,\ 2$, the total spin $S=1$, the total angular momentum (coupled from the
orbital angular momentum and the total spin) $J=1$, and the magnetic
projection $M=0$. Together with the principle quantum number $n$, we have the
complete set of quantum numbers for expressing the 3DHO basis for the deuteron
structure calculation [2, 9, 10, 43]. The matrix element for the kinetic
energy reads [44]:

$\displaystyle\langle
n^{\prime}l^{\prime}S^{\prime}J^{\prime}M^{\prime}|{T}_{\rm rel}|nlSJM\rangle$
$\displaystyle=$
$\displaystyle\frac{\hbar\omega}{2}\delta^{l^{\prime}}_{l}\delta^{S^{\prime}}_{S}\delta^{J^{\prime}}_{J}\delta^{M^{\prime}}_{M}\Big{[}(2n^{\prime}+l^{\prime}+\frac{3}{2})\delta_{n}^{n^{\prime}}+\sqrt{(n^{\prime}+l^{\prime}+\frac{3}{2})(n^{\prime}+1)}\delta^{n^{\prime}}_{n-1}+\sqrt{(n+l+\frac{3}{2})(n+1)}\delta^{n^{\prime}}_{n+1}\Big{]}.$
(2)

It is also convenient to introduce the harmonic oscillator potential ${V}_{\rm
HO}$, which admits the following analytic form [44]:

$\displaystyle\langle
n^{\prime}l^{\prime}S^{\prime}J^{\prime}M^{\prime}|{V}_{\rm HO}|nlSJM\rangle$
$\displaystyle=$
$\displaystyle\frac{\hbar\omega}{2}\delta^{l^{\prime}}_{l}\delta^{S^{\prime}}_{S}\delta^{J^{\prime}}_{J}\delta^{M^{\prime}}_{M}\Big{[}(2n^{\prime}+l^{\prime}+\frac{3}{2})\delta_{n}^{n^{\prime}}-\sqrt{(n^{\prime}+l^{\prime}+\frac{3}{2})(n^{\prime}+1)}\delta^{n^{\prime}}_{n-1}-\sqrt{(n+l+\frac{3}{2})(n+1)}\delta^{n^{\prime}}_{n+1}\Big{]}.$
(3)

The Hamiltonian of the 3DHO then reads:

$\displaystyle H_{\rm HO}={T}_{\rm rel}+{V}_{\rm HO},$ (4)

and the corresponding matrix element is

$\displaystyle\langle
n^{\prime}l^{\prime}S^{\prime}J^{\prime}M^{\prime}|{H}_{\rm
HO}|nlSJM\rangle=(2n^{\prime}+l^{\prime}+\frac{3}{2})\hbar\omega\delta_{n}^{n^{\prime}}\delta^{l^{\prime}}_{l}\delta^{S^{\prime}}_{S}\delta^{J^{\prime}}_{J}\delta^{M^{\prime}}_{M}.$
(5)

We note that both the ${T}_{\rm rel}$ and ${V}_{\rm HO}$ matrices are of the
tri-diagonal form (sparse matrices), while the $H_{\rm HO}$ matrix is diagonal
in the 3DHO basis.

Finding a suitable choice of the reference Hamiltonians is important for
improving the computational efficiency of the ASP-QPE algorithm (cf. Ref. [16]
and references therein). We design a reference Hamiltonian $H_{\rm ref,1}$ to
be our first approximation to the target Hamiltonian $H_{\rm target,1}$. For
the case of the natural deuteron, we model the reference Hamiltonian $H_{\rm
ref,1}$ as a harmonic oscillator (with its mass being the reduced mass of the
two-nucleon system) present in a mean nuclear field ${V}^{(0)}_{\rm NN}$:

$\displaystyle H_{\rm ref,1}={H}_{\rm HO}+{V}^{(0)}_{\rm NN}+{V}_{\rm shift},$
(6)

where ${V}^{(0)}_{\rm NN}$ is taken to be the diagonal entries of $V_{\rm
NN}$, which serves as our approximation to the role of the mean nuclear field.
Note that Hartree-Fock calculations of nuclear structure theory [45, 5] are
conveniently used to provide the mean field. To lift the possible
(near-)degeneracy of the reference Hamiltonian, we also introduce a diagonal
shift Hamiltonian ${V}_{\rm shift}$ into $H_{\rm ref,1}$, which will be
defined below.

An adiabatic path connects $H_{\rm ref,1}$ to $H_{\rm target,1}$. We design
the adiabatic path as:

$\displaystyle{H}_{1}(t)=[1-f(t)]{H}_{\rm ref,1}+f(t){H}_{\rm target,1},$ (7)

where the scalar function $f(t)\in[0,1]$ is defined within the evolution time
duration $t\in[t_{i},t_{f}]$ and it parameterizes the path. For simplicity, we
take

$\displaystyle f(t)=\frac{t-t_{i}}{t_{f}-t_{i}},$ (8)

while other designs of the adiabatic path and choices of the scalar function
are also applicable (cf. Ref. [46]). It is clear that

$\displaystyle{H}_{1}(t_{i})={H}_{\rm ref,1},\ {H}_{1}(t_{f})={H}_{\rm
target,1}.$ (9)

It is worth noting that we select ${H}_{\rm ref,1}$ to be diagonal. The
corresponding eigenstates are the standard bases in the 3DHO representation
(see Eq. (20) below for an example). These eigenstates are taken as the
reference states, which are input to the ASP algorithm.

###  2.2 The trapped deuteron

In order to show our idea of solving the low-lying spectrum of the nuclear
system (including both the ground and excited states) and that of improving
the efficiency of the ASP-QPE algorithm, we also work on the model problem in
which the deuteron is bounded in an external harmonic oscillator trap (note
that the physical, free-space, deuteron has only one bound state). The target
Hamiltonian in the relative coordinates (we factorize and henceforth omit the
center-of-mass motion as illustrated in the problems in Ref. [43]) is

$\displaystyle H_{\rm target,2}=T_{\rm rel}+V_{\rm NN}+{\alpha(V_{\rm
HO}+V_{\rm depth})},$ (10)

where $V_{\rm trap}=\alpha(V_{\rm HO}+V_{\rm depth})$ encapsulates the role of
the external trap on the relative variables with $\alpha$ defining the trap
strength (set to be 0.5). $\alpha V_{\rm depth}$ is a diagonal Hamiltonian
defining the depth of the trap (the diagonal entries of $\alpha V_{\rm depth}$
are constant and set to be -12 MeV).

As our first approximation to this target system, we design the reference
Hamiltonian to be

$\displaystyle H_{\rm ref,2}=H_{\rm HO}+V^{(0)}_{\rm NN}+\alpha V_{\rm
depth}+V_{\rm shift},$ (11)

which represents a harmonic oscillator in the presence of a model for a mean
field $V^{(0)}_{\rm NN}+\alpha V_{\rm depth}$. $V_{\rm shift}$ is again
applied to lift the possible (near-)degeneracy of $H_{\rm ref,2}$ and, again,
it will be defined below. We construct $H_{\rm ref,2}$ to be diagonal such
that its eigenstates are the standard bases in the 3DHO representation.

Similar to the natural deuteron case [Eq. (7)], the adiabatic path for solving
the eigenstates of the trapped deuteron is designed to be

$\displaystyle{H}_{2}(t)=[1-f(t)]{H}_{\rm ref,2}+f(t){H}_{\rm target,2}.$ (12)

where ${H}_{2}(t_{i})=H_{\rm ref,2}$ and ${H}_{2}(t_{f})=H_{\rm target,2}$.

We remark that the proper design of ${H}_{\rm ref,2}$ is important to avoid
the level crossing(s) during that adiabatic evolution (cf. Ref. [16] and
references therein): we find that it would be desirable to include also the
diagonal entries of $V_{\rm trap}$ in ${H}_{\rm ref,2}$ such that no crossing
exists. However, we deliberately exclude these entries in Eq. (11) to retain
the level crossing, such that we can illustrate the alternative idea to treat
the crossing problem.

This idea originates with an insight from many-body theory. In particular,
given that the crossing occurs between the $i^{\rm th}$ and $j^{\rm th}$
levels111The $i^{\rm th}$ and $j^{\rm th}$ levels are the neighboring levels
and we take $j=i+1$. of the instantaneous Hamiltonian ${H}_{2}(t)$ during the
period $[t_{a},t_{b}]$, we propose to introduce an extra perturbation term to
modify the adiabatic path ${H}_{2}(t)$ for this period. Serving as a temporary
perturbation, we require this term to strengthen the level interplay between
the $i^{\rm th}$ and $j^{\rm th}$ levels, such that the crossing is avoided.
As such, we propose the following form of the perturbation term:

$\displaystyle P^{i,j}_{\rm path}(t)=h_{i,j}(t)[1-h_{i,j}(t)]{H}^{i,j}_{\rm
path}.$ (13)

The modified adiabatic path is

$\displaystyle H^{\prime}_{2}(t)=H_{2}(t)+h_{i,j}(t)[1-h_{i,j}(t)]H^{i,j}_{\rm
path}.$ (14)

We require that $P^{i,j}_{\rm path}(t)$ varies slowly and smoothly within
$[t_{a},t_{b}]$, and it vanishes when $t\leq t_{a}$ and $t\geq t_{b}$. For
simplicity, we choose $h_{i,j}(t)={(t-t_{a})}/{(t_{b}-t_{a})}$, with
$t\in[t_{a},t_{b}]$. The term $h_{i,j}(t)[1-h_{i,j}(t)]$ then peaks at
$t=(t_{a}+t_{b})/2$ with the maximum $1/4$. In particular, we have constructed
the $H^{i,j}_{\rm path}$ to be of non-zero matrix elements only in the
subspace $\mathcal{S}$ spanned by the reference states that later cross.
Recall we design $H_{\rm ref,2}$ to be diagonal, and the corresponding
reference states (eigenstates of $H_{\rm ref,2}$) are the standard bases.
Therefore, the subspace $\mathcal{S}$, and the corresponding submatrix, can be
easily determined. For simplicity, we take only the off-diagonal elements in
the submatrix to be non-zero, while setting the diagonal elements to be zeros.
For the remaining entries of ${H}^{i,j}_{\rm path}$, we set them to be zeros
(see in Sec. 5.2 for a detailed example).

The impact of $P^{i,j}_{\rm path}(t)$ on the $i^{\rm th}$ and $j^{\rm th}$
levels is at the second order of perturbation, while its influence on the rest
of the spectrum will be of higher order. We also note that additional path
modifications of the form $P^{i,j}_{\rm path}(t)$ would be necessary when
there are multiple level crossings. We comment that this method will likely
fail when the spectrum (or subspectrum) of the target Hamiltonian is near-
degenerate, since it is not practical to insert the perturbation terms (to
modify the adiabatic path) at the end of the evolution.

##  3 The algorithm

In this section, we show the methodology to compute the deuteron ground state
energy via the ASP-QPE algorithm on the quantum computer. The calculation of
the spectrum of the trapped deuteron is analogous. We also discuss on the
application of the path-modification procedure to improve the computational
efficiency.

###  3.1 The deuteron ground state energy

As illustrated in Fig. 1, the calculation of the deuteron ground state energy
is accomplished in three steps: 1) preparing the deuteron ground state by the
ASP scheme with the qubit(s) in the quantum register $S$ (state); 2) applying
the QPE algorithm to encode the eigenenergy into the relative phase shift on
the probing qubit(s) in the register $R$ (readout); and 3) applying the
inverse Quantum Fourier Transformation ($QFT^{{\dagger}}$) [50] to access the
phase shift and, in turn, the eigenenergy.

Figure 1: (Color online) Illustration of the ASP-QPE algorithm. $|\Psi_{\rm
ref}\rangle$ denotes the reference state, which is typically chosen as the
eigenstate of the reference Hamiltonian. $|\Psi_{\rm ASP}\rangle$ is the trial
eigenstate prepared by the ASP algorithm, which is then input to the QPE
algorithm to evaluate the corresponding eigenenergy.

In the first step, we apply ASP to prepare the trial eigenstate of the
deuteron. To make the ASP algorithm efficient, proper choice of the reference
Hamiltonian is crucial. In particular, one wants 1) to make the reference
state easy to prepare (we design the reference Hamiltonian to be diagonal such
that the reference states are the standard bases); 2) to avoid closely spacing
reference states (we introduce the shift Hamiltonian to lift the
(near-)degeneracy); and 3) to design the reference Hamiltonian such that it
resembles the features of the target Hamiltonian (e.g., in quantum chemistry,
it is advantageous to apply the Hartree-Fock Hamiltonian as the approximation
to the full CI Hamiltonian [15]; we model our reference Hamiltonian as the
harmonic oscillator as an approximation to the role of the mean nuclear field
that arises from the other nucleons).

We adiabatically evolve the ground state of $H_{\rm ref,1}$ along the path
[Eq. (7)]. According to the adiabatic theorem [16, 29, 28], we are guaranteed
to obtain the deuteron ground state at the end of the evolution. This
adiabatic evolution can be solved as [47, 48, 46, 49]:

$\displaystyle|\psi;t_{f}\rangle={U}(t_{f};t_{i})|\psi;t_{i}\rangle=\hat{T}\Bigg{\\{}\exp[-i\int_{t_{i}}^{t_{f}}{H}_{1}(t^{\prime})dt^{\prime}\Big{]}\Bigg{\\}}|\psi;t_{i}\rangle,$
(15)

where ${U}(t_{f};t_{i})$ is the time-evolution operator that propagates the
system from $t_{i}$ to $t_{f}$, while $\hat{T}$ is the time-ordering operator.
$|\psi;t_{i}\rangle=|\Psi_{\rm ref}\rangle$ denotes the reference state
(input), while $|\psi;t_{f}\rangle=|\Psi_{\rm ASP}\rangle$ is the output of
the ASP algorithm (i.e., the trial eigenstate of the deuteron).

We can approximate the complete adiabatic evolution by a sequence of discrete
adiabatic steps (time discretization). That is, we approximate the time-
evolution operator ${U}(t_{f};t_{i})$ as a product of $n$ unitaries [47, 48,
46, 49]:

$\displaystyle{U}(t_{f};t_{i})\approx{U}(t_{n};t_{n-1}){U}(t_{n-1};t_{n-2})\cdots{U}(t_{2};t_{1}){U}(t_{1};t_{0}).$
(16)

The step length $\delta t=(t_{f}-t_{i})/n$ is set to be small such that
$||H_{1}(t)||\cdot\delta t\ll 1$ (we take $\delta t=0.01$ MeV-1 in this work).
For consistency, we set $t_{n}=t_{f}$ and $t_{0}=t_{i}$. The unitary
${U}(t_{k};t_{k-1})$ for the adiabatic step at $t_{k}$ can be further
approximated by the first-order Trotter decomposition [50] based on the
adiabatic path [Eq. (7)]:

$\displaystyle{U}(t_{k};t_{k-1})=e^{-i{H}_{1}(t_{k})\delta
t}=e^{-i[1-f(t_{k})]{H}_{\rm ref,1}\delta t}e^{-if(t_{k})H_{\rm
target,1}\delta t}+\mathcal{O}(\delta t^{2}).$ (17)

In the second step, we apply the QPE algorithm [50] to the trial ground state
$|\Psi_{\rm ASP}\rangle$, which, presumably, overlaps significantly with the
exact deuteron ground state. In practice, we construct the unitary operator
$U_{\rm QPE}=e^{i2\pi\frac{H_{np}+\varphi}{\tau}}$ to encode the deuteron
ground state energy $E$ into the relative phase $\phi$:

$\displaystyle{U}_{\rm QPE}|\Psi_{\rm ASP}\rangle=e^{i2\pi\frac{H_{\rm
target,1}+\varphi}{\tau}}|\Psi_{\rm
ASP}\rangle=e^{i2\pi\frac{E+\varphi}{\tau}}|\Psi_{\rm
ASP}\rangle=e^{i2\pi\phi}|\Psi_{\rm ASP}\rangle.$ (18)

Through repeated controlled actions of powers of ${U}_{\rm QPE}$, the relative
phase is encoded into the bit-string value on the qubits in the $R$ register.
We note that the iterative phase estimation algorithm [14] can improve the
precision using small number of qubits in the $R$ register.

In the last step, we apply the inverse quantum Fourier transformation (QFT†)
[50] to the $R$ register to obtain an approximation of $\phi$ written to $R$
in binary. By classical approximation methods, we obtain an initial estimation
of the ground state energy as $E=\tau\phi-\varphi$. In practice, we choose
$\tau$ and $\varphi$ based on the initial estimation of $E$ such that the
relative phase is $0\leq\phi(\approx\frac{1}{2})<1$.

###  3.2 The trapped deuteron and the modified adiabatic path

In principle, the ASP-QPE approach can also be generalized to solve the low-
lying spectrum of the nuclear systems as well [33]. Analogous to solving the
deuteron ground state energy, we calculate the low-lying spectrum of $H_{\rm
target,2}$, i.e., the deuteron bounded in a harmonic oscillator trap. In
particular, one prepares the trial eigenstates (both the bound and the excited
states) of the trapped deuteron by adiabatically evolving respective reference
states (i.e., the eigenstates of the reference Hamiltonian $H_{\rm ref,2}$)
along the adiabatic path $H_{2}(t)$. With the subsequent implementation of the
QPE algorithm, the eigenenergies of the low-lying states can be obtained.

We note that the level-crossing problem poses challenges to the ASP-QPE
algorithm [15, 28, 34, 35, 51] in solving the spectrum efficiently. While the
proper design of the reference Hamiltonian is helpful, we can also modify the
adiabatic path by introducing perturbation terms to avoid the crossing. As
discussed in Sec. 2.2, such modifications are only effective during the period
of level crossings. They strengthen the interplay between the crossing levels
such that the level crossings are avoided [52]. The efficiency of the ASP-QPE
algorithm can then be improved. We illustrate this idea by solving the
spectrum of the trapped deuteron, where we insert one additional perturbation
term to the adiabatic path $H_{2}(t)$ in order to avoid the level crossing.

##  4 Simulation conditions

In this exploratory work, we retain a limited 4-dimensional model space for
the deuteron channel. For this restricted model space, the principle and
angular momentum quantum numbers of the 3DHO bases are $nl=0s,\ 1s,\ 0d,\
2s$.222We follow the spectroscopic notation in Ref. [5]. The lower case Roman
letters are conventional for basis states: the quantum number of the angular
momentum $l=0$ is denoted by $s$, while $l=2$ is denoted by the $d$. In order
to recover the deuteron ground state energy, we take $V_{\rm NN}$ to be the
effective nuclear interaction obtained from the Daejeon16 inter-nucleon
interaction [53] via the Okubo-Lee-Suzuki similarity transformation [43]:

$\displaystyle V_{\rm
NN}=\begin{pmatrix}-5.72084&-3.51148&-0.164215&-0.533112\\\
-3.51148&-8.72219&-0.0348415&-5.73651\\\
-0.164215&-0.0348415&-8.66007&-0.0436011\\\
-0.533112&-5.736510&-0.0436011&-13.8506\end{pmatrix}.$ (19)

Note that, unlike the effective interaction above, the nuclear interaction for
the applications of complicated nuclei are large and sparse matrices [8]. To
lift the (near-)degeneracy of the reference states, we introduce the diagonal
matrix $V_{\rm shift}$ with entries to be $+3$ MeV ($-3$ MeV) on the first
(second) diagonal entry and alternating again for the third and fourth
diagonal entries.

With the setup in this work, the eigenstates of the diagonal reference
Hamiltonians, $H_{\rm ref,1}$ and $H_{\rm ref,2}$, are the standard bases in
the 3DHO representation:

$\displaystyle|1\rangle=(1,0,0,0)^{T},\ |2\rangle=(0,1,0,0)^{T},\
|3\rangle=(0,0,0,1)^{T},\ |4\rangle=(0,0,1,0)^{T}.$ (20)

The ordering of basis states is, in principle, arbitrary and our chosen
ordering does not produce ordered diagonal entries. Another choice of ordering
should not affect the results of this work.

Analogous to the compact encoding scheme shown in Ref. [49], we map the 3DHO
bases to distinct qubit configurations formed by a sequence of binaries, where
each binary corresponds to one qubit in the quantum register $S$. In our model
problem, the mapping from the standard bases (reference states) to the qubit
configurations is

$\displaystyle|1\rangle\mapsto|00\rangle,\ |2\rangle\mapsto|01\rangle,\
|3\rangle\mapsto|11\rangle,\ |4\rangle\mapsto|10\rangle.$ (21)

With this mapping scheme, the reference states can be easily prepared on the
$S$ register by the action of the Pauli-X gate(s). For more complicated
applications (e.g., heavier nuclei), it takes only $\lceil\log_{2}N_{d}\rceil$
qubits to encode the information of $N_{d}$-dimensional model space.

The quantum circuit of the ASP algorithm is designed following Ref. [49]. In
particular, the unitary ${U}(t_{k};t_{k-1})$ of a discretized adiabatic step
in Eq. (16) is converted into a series of unitaries: some are time-invariant
unitaries that apply in other adiabatic steps, while others are diagonal
unitaries with entries being explicit functions of $t_{k}$. We adopt standard
procedures to decompose each unitary into quantum gates [54, 50]. Combining
these gates in order, we obtain the circuit of ${U}(t_{k};t_{k-1})$, which is
of fixed structure with gate parameters directly parameterized by $t_{k}$.
Indeed, the circuit of ${U}(t_{k};t_{k-1})$ serves as an elementary module:
the complete circuit for the ASP algorithm can be automatically constructed by
sequentially combining these modules (with the parameter $t_{k}$ updated in
each module) on the quantum computer.

With the input reference state in the $S$ register, the corresponding trial
eigenstate of the target Hamiltonian can be obtained by the operation of the
circuit. This trial eigenstate is then input to the QPE algorithm to evaluate
the eigenenergy, where the circuit of the QPE algorithm is designed based on
the target Hamiltonian following standard techniques [54, 50]. We can also
interrupt the adiabatic evolution at intermediate time and design the QPE part
of the circuit based on the instantaneous Hamiltonian. In this way, we can
prepare the trial eigenstates of the instantaneous Hamiltonian and measure the
corresponding eigenenergies during the evolution.

We remark that we do not seek optimal design of the circuit in this initial
work. Our foci are 1) to explore the feasibility of applying the ASP-QPE
algorithm to solve the low-lying spectrum of nuclear systems; and 2) to
explore methods to improve the efficiency of the algorithm. In future work,
the circuit design will be improved. We implement the ASP-QPE algorithm on the
IBM Qiskit quantum simulator [55, 56] and do not take into account the noise
in simulations.

##  5 Results and discussions

###  5.1 The deuteron ground state energy

Figure 2: (Color online) The evolution of eigenenergies of the instantaneous
Hamiltonian $H_{1}(t)$. Solid lines present the results of classical
calculations. The ground state energies at selected moments computed via the
ASP-QPE algorithm implemented with the IBM Qiskit quantum simulator are
presented as black dots. The inset panel shows the squared overlap
$|\langle\Psi_{\rm ASP}(t)|\Psi_{\rm exact}(t)\rangle|^{2}$ of the ground
states of $H_{1}(t)$ at intermediate times during the adiabatic evolution. See
the text for details.

In Fig. 2, we present the classical and quantum calculations of the
eigenenergies of the instantaneous Hamiltonian during the adiabatic evolution.
The classic results of all the four states are obtained by exact matrix
diagonalization, while the quantum results are computed via the ASP-QPE
algorithm for only the ground state.

Starting from $t=0$, where the y-intercepts show the eigenvalue of $H_{\rm
ref}$, we find that the energy gap between the ground state (red line) and the
first excited state (blue line) is easily visible throughout the evolution.
This gap implies good efficiency can be achieved when computing the deuteron
ground state energy via the ASP-QPE algorithm [29, 28].

In the inset panel, we show the squared overlap $|\langle\Psi_{\rm
ASP}(t)|\Psi_{\rm exact}(t)\rangle|^{2}$ between the trial ground state,
$|\Psi_{\rm ASP}(t)\rangle$, prepared by the ASP procedure and the exact
ground state, $|\Psi_{\rm exact}(t)\rangle$, solved by diagonalizing the
instantaneous Hamiltonian $H_{1}(t)$ for every adiabatic step. Both the states
and the squared overlap are computed on the classical computer. The squared
overlap values indicate the quality of the trial ground states at intermediate
times: better quality of the trial ground states ensures higher probability
for the QPE algorithm to observe the ground state energies during the
evolution. We find that the ASP procedure well prepares the instantaneous
ground state along the adiabatic path: the squared overlap $|\langle\Psi_{\rm
ASP}(t)|\Psi_{\rm exact}(t)\rangle|^{2}$ is above 93$\%$ throughout the
evolution. We also checked that these values approach 100$\%$ given moderately
longer evolution time (e.g., at $T=10$ MeV-1, versus $T=1$ MeV-1 used here,
the squared overlap values are all beyond 99.5$\%$).

Implementing the IBM Qiskit quantum simulator [55, 56], we compute the ground
state energies at selected moments along the adiabatic path $H_{1}(t)$ based
on the ASP-QPE algorithm. These results are shown as the black dots in Fig. 2.
Note that the errors of these quantum results are mainly from the statistical
variance during the measurement, which are smaller than the size of the
markers. Therefore, we do not show the error bar explicitly. We find the
quantum results (black dots) agree well with the classical results (red line).
At $t=T=1$ MeV-1, the deuteron ground state energy is recovered. To obtain
more precise result of the deuteron ground state energy making use of limited
number of qubits, a recursive phase estimation algorithm [14] will be useful.

###  5.2 The deuteron in the harmonic oscillator trap

Figure 3: (Color online) Evolution of the spectra via the normal adiabatic
path ${H}_{2}(t)$ (short dot-dashed lines) and via the modified adiabatic path
$H^{\prime}_{2}(t)$ (long dashed lines). The spectra at selected moments in
time computed via the ASP-QPE algorithm with the modified adiabatic path
$H_{2}^{\prime}(t)$ implementing the IBM Qiskit quantum simulator are shown as
colored dots. With the normal adiabatic path, the level crossing occurs at
$t\approx 2.9$ MeV-1 between the $2^{nd}$ and $3^{rd}$ excited states, with
the minimal gap size being $\Delta E\approx 0.02$ MeV.

The ASP-QPE algorithm can be applied to solve the eigenenergies of the excited
states in addition to the ground state energy. We illustrate this scheme by
solving the spectrum of the deuteron system bounded in an external harmonic
oscillator trap.

We again apply the ASP algorithm to prepare the trial eigenstates of the
system: each trial eigenstate (output) is evolved from the corresponding
reference state (input) via the adiabatic path. The straightforward choice of
the path is $H_{2}(t)$. By exact diagonalization of the instantaneous
Hamiltonian, we compute the spectrum of $H_{2}(t)$ as functions of time. The
results are shown as the short dot-dashed lines in Fig. 3. Note that the
y-intercepts correspond to the eigenvalues of ${H}_{\rm ref,2}$, which are
indeed the diagonal entries of ${H}_{\rm ref,2}$ [Eq. (11)]. We find that the
minimal energy gap between the 2nd and 3rd excited states occurs at
$t_{c}\approx 2.9$ MeV-1: such a level crossing causes a transition between
these two states and hinders the efficiency of the ASP algorithm [15, 28, 34,
35, 51, 33].

In order to avoid the crossing, we modify the adiabatic path. In particular,
we introduce an extra perturbation term to the adiabatic path $H_{2}(t)$
during the period $[t_{a},t_{b}]$ when the levels approach (here we take
$t_{a}=0$ and $t_{b}=6$ MeV-1). We note that these crossing levels evolve from
the reference states $|3\rangle=(0,0,0,1)^{T}$ and $|4\rangle=(0,0,1,0)^{T}$,
which span the subspace $\mathcal{S}=\\{|3\rangle,|4\rangle\\}$. Therefore, we
design the perturbation term to be $P_{\rm
path}^{3,4}=h_{3,4}(t)[1-h_{3,4}(t)]H_{\rm path}^{3,4}$, where we take the
matrix elements $[H_{\rm path}^{3,4}]_{34}=[H_{\rm path}^{3,4}]_{43}=-2$ MeV,
while the remaining elements of the matrix $[H_{\rm path}^{3,4}]_{4\times 4}$
are all zeros.

The modified adiabatic path $H^{\prime}_{2}(t)$ connects the same reference
and target Hamiltonians as ${H}_{2}(t)$. We compute the spectrum of
$H^{\prime}_{2}(t)$ by exact diagonalization during the evolution. The results
are shown as the long dashed lines in Fig. 3. The additional perturbation
avoids the level crossing by strengthening the interplay between between the
2nd and 3rd excited states. We remark that we have also applied this path-
modification scheme to the same model problem but with an enlarged
(64-dimensional) model space; we verified the feasibility of this scheme to
avoid multiple level crossings during the evolution.

Figure 4: (Color online) Evolution of the squared overlap $|\langle\Psi_{\rm
ASP}(t)|\Psi_{\rm exact}(t)\rangle|^{2}$ for each state based on the normal
path $H_{2}(t)$ (short dot-dashed lines) and the modified path
$H^{\prime}_{2}(t)$ (long dashed lines) during the evolution. The magenta long
dashed/short dot-dashed lines overlap with their respective gray lines, while
all the blue lines overlap with all the red lines.

Based on either the normal path $H_{2}(t)$ or the modified path
$H^{\prime}_{2}(t)$, we compute for every intermediate moment $t$: 1) the
$k^{th}$ ($k=1,\ 2,\ 3,\ 4$) trial eigenstate $|\Psi_{\rm ASP;k}(t)\rangle$
evolved from the $k^{th}$ reference state by the ASP procedure; and 2) the
$k^{th}$ exact eigenstate $|\Psi_{\rm exact;k}(t)\rangle|^{2}$ obtained by
matrix diagonalization. We then compute the squared overlap $|\langle\Psi_{\rm
ASP;k}(t)|\Psi_{\rm exact;k}(t)\rangle|^{2}$ for each state. Note that these
calculations are all done on a classical computer. These values characterize
the quality of the trial eigenstates prepared by the ASP algorithm.

We show the squared overlap values in Fig. 4. For the evolution via the normal
path, we note that the squared overlap values of the ground state and of the
$1^{st}$ excited state approach 1. However, the squared overlap values of
$2^{nd}$ and $3^{rd}$ excited states drop to about zero starting from the
crossing point at $t\approx 2.9$ MeV-1 (seen from the overlapping magenta and
gray dot-dashed lines): a significant crossover transition takes place between
these two states with the limited total evolution time. The small energy gap
and the resulting transition hinder the efficiency of the ASP-QPE algorithm.

For comparison, we present the evolution of the squared overlap value of each
state computed based on the modified path $H^{\prime}_{2}(t)$ in Fig. 4. We
find that the squared overlap values of the ground state and the $1^{st}$
excited state overlap with those computed based on the norm path $H_{2}(t)$
(they all nearly approach 1). We also note that a significant crossover
transition is prevented since the level crossing is avoided: 1) the squared
overlap values of the $2^{nd}$ and $3^{rd}$ excited states remain above 0.75
during the evolution; 2) the squared overlap values of these two states are
restored to more than 0.97 at the end of the evolution. We also checked to see
that, with longer evolution time (e.g., $T\geq 50$ MeV-1), the squared overlap
values for all the states computed based on the modified path
$H^{\prime}_{2}(t)$ are more than 0.99 throughout the adiabatic evolution; for
achieving the same quality, it would take $T\gg\frac{1}{\Delta E^{2}}\approx
2500$ MeV-1 [29] for the adiabatic evolution via the normal path $H_{2}(t)$.
Thus, the efficiency of the ASP algorithm can be significantly improved by the
path modification scheme.

The trial eigenstate $|\Psi_{\rm ASP;k}(t)\rangle$ prepared based on the
modified adiabatic path $H^{\prime}_{2}(t)$ is input to the QPE algorithm for
evaluating the eigenenergy. For the total evolution time $T=10$ MeV-1 applied
in this work, the trial eigenstates overlap significantly with the
corresponding exact eigenstates in the intermediate times and at the end of
the evolution. Therefore, one observes the eigenenergies of the corresponding
trial eigenstates with large probabilities.

Implementing the IBM Qiskit quantum simulator [55, 56], we compute the spectra
of $H^{\prime}_{2}(t)$ implementing the ASP-QPE algorithm at selected times
during the evolution. The results are presented as the colored dots in Fig. 3.
We find that the results obtained via the quantum computation agree well with
the classical results. We do not show the error bar (which is mainly from the
statistical variance in measurement and of the size smaller than the marker
size). At the end of the evolution, the ASP-QPE algorithm successfully
recovers the eigenenergies of the target Hamiltonian.

Finally, we notice that it is important to know where and which neighboring
levels cross in order to implement the ASP-QPE algorithm efficiently. In
quantum computations, we can achieve this by energy measurement at
intermediate times along the adiabatic path. While the decreasing energy gap
is one of the signs of the level crossing, multiple eigenvalues observed with
comparable probabilities on one trial eigenstate, which persists even with
longer evolution time, will also be indicative.

In many applications, we are interested in a few low-lying states of the
nuclear system. Once the crossing(s) is (are) detected, we can implement the
path-modification scheme to avoid the crossing(s) and work to improve the
computational efficiency. We expect that the ASP-QPE algorithm, together with
the path modification scheme, will be applicable for solving the (sub)spectra
with moderately separated states. However, the algorithm will be less
efficient as the spectrum becomes more dense (e.g., with near-degenerate
states of the target Hamiltonian), where the energy gaps between the levels of
the target Hamiltonian are small.

##  6 Summary and outlook

Solving nuclear many-body problems using an ab initio approach is a
computationally hard problem. For complex nuclei with increasing atomic
number, this task will be intractable even for world-leading supercomputers.
Quantum computers potentially open a promising path to address such nuclear
many-body problems, and there is hence a need to develop the quantum
algorithms for ab initio nuclear structure theory. In this work, we
investigate the application of the quantum algorithm of the adiabatic state
preparation with quantum phase estimation (ASP-QPE) in solving ab initio
nuclear structure theory.

We first apply the ASP-QPE algorithm to solve the ground state energy of the
natural deuteron. We apply the ASP algorithm to prepare the trial ground
state. Proper design of the reference Hamiltonian is important for the
algorithm to be efficient. In this work, we model the reference Hamiltonian as
a harmonic oscillator (its mass being the reduced mass of the two-nucleon
system) existing in the mean nuclear field. We design the adiabatic path to
connect the reference Hamiltonian and the deuteron Hamiltonian. By
adiabatically evolving the ground state of the reference Hamiltonian along the
path, we prepare the trial ground states of the instantaneous Hamiltonians.
These trial states are input to the QPE algorithm for evaluating the
instantaneous ground state energies. We adopt a compact mapping scheme to map
the model problem onto the quantum register. The quantum circuit of the ASP
algorithm is designed following the idea of Ref. [49] and the circuit of the
QPE algorithm is designed based on the instantaneous Hamiltonians. We
implement the IBM Qiskit quantum simulator to simulate the circuits and the
quantum results agree well with the classical results solved by exact matrix
diagonalization. At the end of the evolution, we recover the deuteron
eigenenergy.

We then generalize the ASP-QPE algorithm to solve the low-lying spectra (both
the ground and excited states) of simple nuclear systems. We note that the
efficiency of the algorithm is hindered if there is a small/vanishing energy
gap (level crossing). While suitable design of the reference Hamiltonian is
helpful to address this issue, we also introduce the idea of path modification
by insertions of intermediate perturbation terms. In particular, we design
these perturbation terms such that they are only effective during the periods
of level crossing and require these terms to strengthen the interplay between
the approaching levels. The level crossing can be avoided by this path-
modification scheme, via which process we manage to improve the efficiency of
the ASP-QPE algorithm.

As an illustration, we solve the spectrum of the deuteron bounded in an
external harmonic oscillator trap. The application of the ASP-QPE algorithm is
analogous to that of the natural deuteron, except for the modified path to
avoid the level crossing. Implementing the IBM Qiskit quantum simulator, we
solve the eigenenergies of the instantaneous Hamiltonians at selected times.
All of the quantum results agree well with the exact results obtained by
matrix diagonalization. At the end of the evolution, the spectrum of the
bounded deuteron is recovered.

In the future, we plan to develop the ASP-QPE algorithm to solve the low-lying
spectra of complicated nuclei. Adopting the compact encoding scheme, it takes
only $\lceil\log_{2}N_{d}\rceil$ qubits to encode the $N_{d}$-dimensional
model problem. Meanwhile, proper design of the reference Hamiltonian and the
path-modification scheme help to improve the efficiency. Since all the
Hamiltonians involved in constructing the adiabatic path and in the phase
estimation procedure are (or can be divided into) large sparse matrices with
elements explicitly computable, an efficient quantum circuit implementing the
ASP-QPE algorithm can be designed [36, 37, 57] and we expect that quantum
advantage can be achieved.

## Acknowledgments

This work was supported in part by the US Department of Energy (DOE) under
Grants No. DE-FG02-87ER40371 and No. DE-SC0018223 (SciDAC-4/NUCLEI). A portion
of the computational resources were provided by the National Energy Research
Scientific Computing Center (NERSC), which is supported by the US DOE Office
of Science. X.Z. and W.Z. are supported by the Strategic Priority Research
Program of Chinese Academy of Sciences, Grant No. XDB34000000. X.Z. is
supported by new faculty startup funding from the Institute of Modern Physics,
Chinese Academy of Sciences and by Key Research Program of Frontier Sciences,
CAS, Grant No. ZDBS-LY-7020. W.Z. are supported by the National Natural
Science Foundation of China (Grants No. 11975282 and No. 11435014) and the 973
Program of China (Grant No. 2013CB834405).

## References

  * [1] S. C. Pieper and R. B. Wiringa, Ann. Rev. Nucl. Part. Sci. 51 (2001), 53-90 doi:10.1146/annurev.nucl.51.101701.132506 [arXiv:nucl-th/0103005 [nucl-th]].
  * [2] B. R. Barrett, P. Navratil and J. P. Vary, Prog. Part. Nucl. Phys. 69, 131-181 (2013) doi:10.1016/j.ppnp.2012.10.003
  * [3] G. Hagen, T. Papenbrock, M. Hjorth-Jensen and D. J. Dean, Rept. Prog. Phys. 77 (2014) no.9, 096302 doi:10.1088/0034-4885/77/9/096302 [arXiv:1312.7872 [nucl-th]].
  * [4] E. Caurier, G. Martinez-Pinedo, F. Nowacki, A. Poves and A. P. Zuker, Rev. Mod. Phys. 77, 427-488 (2005) doi:10.1103/RevModPhys.77.427 [arXiv:nucl-th/0402046 [nucl-th]].
  * [5] J. Suhonen, From Nucleons to Nucleus: Concepts of Microscopic Nuclear Theory (Springer-Verlag, Berlin, 2007)
  * [6] J. P. Vary, P. Maris, E. Ng, C. Yang and M. Sosonkina, J. Phys. Conf. Ser. 180 (2009), 012083 doi:10.1088/1742-6596/180/1/012083 [arXiv:0907.0209 [nucl-th]].
  * [7] P. Maris, J. P. Vary and A. M. Shirokov, Phys. Rev. C 79 (2009), 014308 doi:10.1103/PhysRevC.79.014308 [arXiv:0808.3420 [nucl-th]].
  * [8] J. P. Vary, R. Basili, W. Du, M. Lockner, P. Maris, D. Oryspayev, S. Pal, S. Sarker, H. M. Aktulga, E. Ng, M. Shao and C. Yang, Proceedings of the International Conference ‘Nuclear Theory in the Supercomputing Era-2016’ (NTSE-2016), Khabarovsk, Russia, September 19-23, 2016. Eds. A. M. Shirokov and A. I. Mazur. Pacific National University, Khabarovsk, Russia, 2018, p. 15. [arXiv:1803.04101 [nucl-th]].
  * [9] P. Navratil, J. P. Vary and B. R. Barrett, Phys. Rev. Lett. 84, 5728-5731 (2000) doi:10.1103/PhysRevLett.84.5728 [arXiv:nucl-th/0004058 [nucl-th]].
  * [10] P. Navratil, J. P. Vary and B. R. Barrett, Phys. Rev. C 62, 054311 (2000) doi:10.1103/PhysRevC.62.054311
  * [11] R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982). doi:10.1007/BF02650179
  * [12] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, and Alán Aspuru-Guzik, Chem. Rev. 119, 10856 (2019).
  * [13] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020).
  * [14] A. Aspuru-Guzik, A. Dutoi, P. Love and M. Head-Gordon, Science 309, 1704 (2005).
  * [15] J. D. Whitfield, J. Biamonte, A. Aspuru-Guzik, Molecular Physics, 109, 735 (2011) [arXiv:1001.3855 [quant-ph]].
  * [16] T. Albash and D. A. Lidar, Rev. Mod. Phys. 90 015002 (2018).
  * [17] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nat. Commun. 5, 4213 (2014).
  * [18] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. 18, 023023 (2016).
  * [19] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nat. Commun. 10, 3007 (2019).
  * [20] E. Farhi, J. Goldstone, and S. Gutmann, arXiv:1411.4028 (2014).
  * [21] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. Kin-Lic Chan, Nat. Phys. 16, 205–210 (2020).
  * [22] S. McArdle, T. Jones, S. Endo, Y. Li, S. Benjamin, and X. Yuan, npj Quantum Information 5, 75 (2019).
  * [23] K. Yeter-Aydeniz, R. C. Pooser, and G. Siopsis, npj Quantum Information 6, 63 (2020).
  * [24] E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean, P. Lougovski, Phys. Rev. Lett. 120, 210501 (2018).
  * [25] O. Di Matteo, A. McCoy, P. Gysbers, T. Miyagi, R. M. Woloshyn and P. Navrátil, [arXiv:2008.05012 [quant-ph]].
  * [26] J. Preskill, Quantum 2, 79 (2018), arXiv:1801.00862 [quant-ph].
  * [27] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nat. Commun. 9, 4812 (2018).
  * [28] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren and D. Preda, Science 292, 472 (2001), quant-ph/0104129.
  * [29] A. Messiah, Quantum mechanics: Volume II (North-Holland Publishing Company Amsterdam, 1962).
  * [30] A. Y. Kitaev, arXiv:quant-ph/9511026.
  * [31] D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 79, 2586-2589 (1997) doi:10.1103/PhysRevLett.79.2586 [arXiv:quant-ph/9703054 [quant-ph]].
  * [32] D. S. Abrams and S. Lloyd, Phys. Rev. Lett. 83, 5162-5165 (1999) doi:10.1103/PhysRevLett.83.5162 [arXiv:quant-ph/9807070 [quant-ph]].
  * [33] J. Liu, Y. Xin, arXiv: 2004.13234.
  * [34] E. Farhi, J. Goldstone, S. Gutmann, arXiv: quant-ph/0208135.
  * [35] E. Farhi, J. Goldstone, D. Gosset, S. Gutmann, H. B. Meyer, and P. Shor, arXiv: 0909.4766 [quant-ph].
  * [36] D. Aharonov, A. Ta-Shma, arXiv:quant-ph/0301023.
  * [37] Stephen P. Jordan and Pawel Wocjan, Phys. Rev. A 80, 062301 (2009).
  * [38] R. Machleidt, Adv. Nucl. Phys. 19, 189 (1989).
  * [39] R. Machleidt et al., Phys. Rep. 149, 1 (1987).
  * [40] D. R. Entem1, and R. Machleidt, Phys. Rev. C 68 041001(R) (2003).
  * [41] E. Epelbaum, H. W. Hammer and U. G. Meissner, Rev. Mod. Phys. 81, 1773-1825 (2009) doi:10.1103/RevModPhys.81.1773 [arXiv:0811.1338 [nucl-th]].
  * [42] E. Epelbaum, H. Krebs and U. G. Meißner, Eur. Phys. J. A 51, no.5, 53 (2015) doi:10.1140/epja/i2015-15053-8 [arXiv:1412.0142 [nucl-th]].
  * [43] J. P. Vary, R. Basili, W. Du, M. Lockner, P. Maris, S. Pal and S. Sarker, Phys. Rev. C 98, no.6, 065502 (2018) doi:10.1103/PhysRevC.98.065502 [arXiv:1809.00276 [nucl-th]].
  * [44] W. Du, P. Yin, Y. Li, G. Chen, W. Zuo, X. Zhao and J. P. Vary, Phys. Rev. C 97, no.6, 064620 (2018) doi:10.1103/PhysRevC.97.064620 [arXiv:1804.01156 [nucl-th]].
  * [45] J. P. Delaroche, M. Girod, J. Libert, H. Goutte, S. Hilaire, S. Peru, N. Pillet and G. F. Bertsch, Phys. Rev. C 81, 014303 (2010) doi:10.1103/PhysRevC.81.014303 [arXiv:0910.2940 [nucl-th]].
  * [46] B. F. Schiffer, J. Tura, and J. I. Cirac, arXiv:2103.01226 (2021).
  * [47] L.-A. Wu, M. S. Byrd, and D. A. Lidar, Phys. Rev. Lett. 89, 057904 (2002).
  * [48] Y. Sun, J. Zhang, M. S. Byrd, and L.-A. Wu, New J. Phys. 22 053012 (2020).
  * [49] W. Du, J. P. Vary, X. Zhao and W. Zuo, [arXiv:2006.01369 [nucl-th]].
  * [50] M. Nielsen and I. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, UK, 2001).
  * [51] N. G. Dickson, and M. H. S. Amin, Phys. Rev. Lett. 106, 050502 (2011).
  * [52] C. Cohen-Tannoudji, B. Diu, and F. Laloe, Quantum Mechanics, vol 1 (Wiley-VCH, 1992).
  * [53] Y. Kim, I. J. Shin, A. M. Shirokov, M. Sosonkina, P. Maris and J. P. Vary, [arXiv:1910.04367 [nucl-th]].
  * [54] C. Li, R. Roberts, X. Yin, arXiv:1210.7366 [quant-ph].
  * [55] A. C. Santos, Rev. Bras. Ens. Fis. 39(1), e1301 (2017). arXiv:1610.06980 [quant-ph].
  * [56] A. W. Cross, L. S. Bishop, J. A. Smolin, and J. M. Gambetta, arXiv preprint arXiv:1707.03429 (2017).
  * [57] D. W. Berry, G. Ahokas, R. Cleve, B. C. Sanders, Communications in Mathematical Physics 270, 359 (2007).