# Nuclear shell-model simulation in digital quantum computers

A. Pérez-Obiol<EMAIL_ADDRESS>Barcelona Supercomputing Center, 08034
Barcelona, Spain  A. M. Romero<EMAIL_ADDRESS>J. Menéndez
<EMAIL_ADDRESS>A. Rios<EMAIL_ADDRESS>Departament de Física
Quàntica i Astrofísica (FQA), Universitat de Barcelona (UB), c. Martí i
Franqués, 1, 08028, Barcelona, Spain
Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB), c.
Martí i Franqués, 1, 08028 Barcelona, Spain  A. García-Sáez
<EMAIL_ADDRESS>Barcelona Supercomputing Center, 08034 Barcelona, Spain
Qilimanjaro Quantum Tech, 08007 Barcelona, Spain  B. Juliá-Díaz
<EMAIL_ADDRESS>Departament de Física Quàntica i Astrofísica (FQA),
Universitat de Barcelona (UB), c. Martí i Franqués, 1, 08028 Barcelona, Spain
Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (UB), c.
Martí i Franqués, 1, 08028 Barcelona, Spain

###### Abstract

The nuclear shell model is one of the prime many-body methods to study the
structure of atomic nuclei, but it is hampered by an exponential scaling on
the basis size as the number of particles increases. We present a shell-model
quantum circuit design strategy to find nuclear ground states by exploiting an
adaptive variational quantum eigensolver algorithm. Our circuit implementation
is in excellent agreement with classical shell-model simulations for a dozen
of light and medium-mass nuclei, including neon and calcium isotopes. We
quantify the circuit depth, width and number of gates to encode realistic
shell-model wavefunctions. Our strategy also addresses explicitly energy
measurements and the required number of circuits to perform them. Our
simulated circuits approach the benchmark results exponentially with a
polynomial scaling in quantum resources for each nucleus. This work paves the
way for quantum computing shell-model studies across the nuclear chart and our
quantum resource quantification may be used in configuration-interaction
calculations of other fermionic systems.

††preprint: APS/123-QED

Atomic nuclei are complex many-body systems formed by protons and neutrons
(collectively denoted as nucleons) bound by the strong nuclear force. Nuclei
exhibit captivating properties such as the coexistence of spherical and
deformed shapes at low energies Taniuchi  _et al._ (2019); Butler  _et al._
(2019); Tsunoda  _et al._ (2020), strong short-range correlations between
pairs of nucleons Schmidt  _et al._ (2020), or decay modes driven by the
strong Mukha  _et al._ (2006), weak Hinke  _et al._ (2012) or electromagnetic
Walz  _et al._ (2015) forces. Furthermore, nuclear decays are crucial to
understand the origin of heavy elements in the universe Cowan  _et al._
(2021), and experiments using nuclei aim to answer fundamental physics
questions such as which is the nature of dark matter Aalbers  _et al._ (2023),
why matter dominates over antimatter in the universe Engel  _et al._ (2013),
or whether neutrinos are their own antiparticles Avignone  _et al._ (2008).

The nuclear shell model, also known as the configuration interaction method,
is one of the leading many-body approaches to study the structure of nuclei.
The shell model is grounded in the idea that, in a similar fashion to
electrons in an atom, nucleons occupy orbitals organized in _shells_ of
different energies Mayer (1949); Haxel  _et al._ (1949). Nuclear states are
then obtained by computationally intensive diagonalizations of the nuclear
Hamiltonian in a many-body configuration space comprising one or several
shells. In spite of impressive progress in recent decades Brown and Wildenthal
(1988); Caurier  _et al._ (2005); Otsuka  _et al._ (2020); Stroberg  _et al._
(2019), the exponential scaling of the many-body Hilbert space with the number
of nucleons ultimately prevents the application of the shell model across the
entire nuclear chart, particularly in heavy nuclei.

Quantum computing promises to circumvent limitations associated to any
exponentially-scaling many-body system using the principle of superposition of
qubit states Arute  _et al._ (2019). In the current noisy intermediate-scale
quantum (NISQ) device era Preskill (2018), _variational quantum eigensolvers_
(VQE) Peruzzo  _et al._ (2014); McClean  _et al._ (2016) are among the most
successful algorithms Bharti  _et al._ (2022) exploiting the benefits of
quantum computing to deal with complex many-body problems in physics Cerezo
_et al._ (2021); Tilly  _et al._ (2022) and chemistry Anand  _et al._ (2022);
McArdle  _et al._ (2020); Haidar  _et al._ (2022).

Quantum many-body systems that have been used as VQE testbeds include the
Fermi-Hubbard Cade  _et al._ (2020), Ising Cervera-Lierta (2018) and Lipkin-
Meshkov-Glick models Cervia  _et al._ (2021); Harsha  _et al._ (2018); Faba
_et al._ (2022); Wahlen-Strothman  _et al._ (2017); Robin and Savage (2023),
superfluid systems Lacroix (2020); Ruiz Guzman and Lacroix (2022), hadrons
Qian  _et al._ (2022) or molecules Grimsley  _et al._ (2019); Sapova and
Fedorov (2022); Feniou  _et al._ (2023).

Figure 1:  The shell model and quantum encoding. Panel (a): Segrè chart
covering the $p$, $sd$ and part of the $pf$ shell. Solid lines indicate
neutron and proton magic numbers. Open circles show the isotopes studied in
this work. Panel (b): schematic representation of the $p$-, $sd$\- and
$pf$-shell configuration spaces. The number on top of every single-particle
state is the qubit label for the implementation in a quantum device under a
Jordan-Wigner mapping. Panel (c): number of many-body configurations,
$\dim_{\rm{mb}}$, in the $M$-basis as a function of the number of active
neutrons in the configuration space, $N_{\rm{CI}}$. We show results for the
isotopic chains of He and Be in the $p$ shell; O, F, Ne, and Al in the $sd$
shell; and Ca, Ti, Cr, and Zn in the $pf$ shell. Isotopes beyond the middle of
the shell are not shown since the number of configurations is symmetric. Bold
marker lines highlight nuclei studied in this work.

In general, a VQE implementation requires a series of well-defined stages
Tilly  _et al._ (2022), involving a) a mapping between physical degrees of
freedom (eg fermionic operators) and the qubits in a quantum computer; b) the
preparation of an initial reference state; c) a (potentially iterative)
variational optimization; d) a measurement strategy for expectation values of
operators (most importantly, the Hamiltonian); and e) an error mitigation
scheme. Previous nuclear shell-model studies have only partially tackled these
problems Dumitrescu  _et al._ (2018); Lu  _et al._ (2019); Stetcu  _et al._
(2022a); Kiss  _et al._ (2022). The aim of this article is to present a
circuit design strategy that explicitly addresses all these aspects to solve
the nuclear shell model in a quantum computer. We also quantify the necessary
circuit resources, such as depths and widths, to achieve precise predictions
for nuclear masses. We do this in a set of test nuclei across different
nuclear shells. To this end, we perform (classical) baseline simulations on
the corresponding circuit architectures and benchmark the results against
diagonalizable shell-model simulations as well as independent ADAPT-VQE
simulations without an explicit circuit implementation.

## Results

### Nuclear shell model

The nuclear shell model Brown and Wildenthal (1988); Caurier  _et al._ (2005);
Otsuka  _et al._ (2020); Stroberg  _et al._ (2019) considers nuclei composed
by an inert core of nucleons, which do not explicitly contribute to the
dynamics, and a set of valence protons and neutrons interacting in a
relatively small configuration space. This space is usually bounded by two
_magic numbers_ , which denote special configurations of protons or neutrons
leading to particularly stable nuclei. Magic numbers thus define shells with
large energy gaps between them. Configuration spaces used in shell-model
calculations usually comprise one or two shells. Panel (a) of Fig. 1 shows the
light to mid-mass region of the isotope chart. We highlight areas where the
$p$, $sd$ and $pf$ shell-model calculations are routinely employed.

Since the nuclear force is rotationally invariant and nucleons are fermions,
it is useful to work in a single-particle basis with states with quantum
numbers $n\,l_{j}$, where $n$ is the principal quantum number, $l$ the orbital
angular momentum and $j$ the total angular momentum Shalit and Talmi (1963).
This basis also includes $m$ third-component projections of $j$ degenerate in
energy. The nuclear Hamiltonian is also to a very good approximation the same
for neutrons and protons, so it is customary to define, additionally, the
isospin quantum number $t=1/2$, with third component $t_{z}$ discerning
protons and neutrons Talmi (1993). Many-body nuclear states have good total
angular momentum $J$ and isospin $T$, with respective third components $M$ and
$T_{z}$ given by the sum of the third components of all nucleons in the
nucleus Varshalovich  _et al._ (1988).

The nuclear Hamiltonian in a given configuration space can be written as

$H_{\rm{eff}}=\sum_{i}\varepsilon_{i}a_{i}^{{\dagger}}a_{i}+\frac{1}{4}\sum_{ijkl}\bar{v}_{ijkl}a_{i}^{{\dagger}}a_{j}^{{\dagger}}a_{l}a_{k}\,,$
(1)

where $\varepsilon_{i}$ is the energy of the single-particle state $i$ and
$\bar{v}_{ijkl}=v_{ijkl}-v_{ijlk}$ are antisymmetrized two-body matrix
elements. $a_{i}$ and $a_{i}^{{\dagger}}$ are fermionic annihilation and
creation operators associated to each single-particle state, $i$. The matrix
elements $\bar{v}_{ijkl}$ can be obtained Hjorth-Jensen  _et al._ (1995);
Stroberg  _et al._ (2019) from an effective field theory of the underlying
theory of the nuclear force, quantum chromodynamics Epelbaum  _et al._ (2009).
Here, instead, we use standard phenomenological Hamiltonians, with components
adjusted to better reproduce key properties of selected nuclei Poves and Zuker
(1981). We choose the Cohen-Kurath interaction in the $p$ shell Cohen and
Kurath (1965), USDB in the $sd$ shell Brown and Richter (2006) and KB3G in the
$pf$ shell Poves  _et al._ (2001).

A suitable many-body basis, also referred to as Fock space, for shell-model
calculations is provided by the so-called $M-$scheme Talmi (1993), in which
the Slater determinant states are chosen to have a well-defined $M$.
$T_{z}=(N-Z)/2$ is also well defined because the number of neutrons $N$ and
protons $Z$ is fixed. Nuclear states are thus expanded in this basis,

$|JM\,TT_{z}\rangle=\sum_{\alpha}c_{\alpha}|\alpha,MT_{z}\rangle,$ (2)

and nuclear wavefunctions and their corresponding energies are eigenvectors
and eigenvalues of the Hamiltonian matrix in the basis of Slater determinants.
The $c_{\alpha}$ coefficients are obtained through diagonalization employing
state-of-the-art nuclear shell-model codes Caurier and Nowacki (1999); Shimizu
_et al._ (2019); Brown and Rae (2014); Johnson  _et al._ (2018) and ensure
that eigenstates have good $J$ and $T$ quantum numbers.

However, this framework faces a steep computational bottleneck in terms of the
maximum size of the Hamiltonian matrix from which the lowest eigenvalues and
eigenvectors can be calculated. The dimension of the single-particle basis of
a nuclear shell consisting of several orbitals $nl_{j}$ is

$\displaystyle\dim_{\rm{sp}}=\sum_{j}(2j+1),$ (3)

where the sum runs over the $j$ values in a given configuration space, see
panel (b) of Fig. 1 for details. The corresponding number of Slater
determinants grows combinatorially as

$\dim_{\text{mb}}=\binom{\dim_{\text{sp}}}{N_{\text{CI}}}\times\binom{\dim_{\text{sp}}}{Z_{\text{CI}}},$
(4)

where $N_{\text{CI}}$ ($Z_{\text{CI}}$) is the number of active neutrons
(protons) in the configuration space. Let us consider the $sd$ shell,
comprising the $1s_{1/2}$, $0d_{3/2}$ and $0d_{5/2}$ orbitals for both protons
and neutrons, and the $pf$ shell, comprising the $0f_{7/2}$, $0f_{5/2}$,
$1p_{3/2}$ and $1p_{1/2}$ orbitals. There are $12$ ($20$) single-particle
states in the $sd$ ($pf$) shell, so that it can describe the isotopic chains
of $12$ ($20$) elements with up to $12$ ($20$) valence neutrons, as shown in
panel (a) of Fig. 1. Panel (c) illustrates the exponential scaling of the
number of many-body configurations, $\dim_{\rm{mb}}$, present for isotopes of
elements in different shells. The number of basis states needed to describe
two isotopes of the same element, or two elements with the same $N$ in the
same shell, can differ by three or more orders of magnitude.

In practical calculations, this number may be reduced by about an order of
magnitude due to symmetry considerations, leading to a reduced number of
Slater determinants, $N_{\text{SD}}$ Caurier  _et al._ (2005). However, the
scaling in either $\dim_{\text{mb}}$ or $N_{\text{SD}}$ ultimately places a
limit in the computational resources needed to study heavy nuclei with the
nuclear shell model. This refers to both the number of operations per second,
or CPU time, and the memory to store all configurations. In fact, the shell-
model history is closely tied to that of computation, as larger-scale
calculations became feasible with the advances in computational power and
refined techniques in CPUs and GPUs Caurier  _et al._ (2005); Brown and
Wildenthal (1988); Otsuka  _et al._ (2020); Stroberg  _et al._ (2019).

### Variational algorithm

Here, we implement the nuclear shell model in a quantum computer following a
standard Jordan-Wigner (JW) mapping Seeley  _et al._ (2012); Romero  _et al._
(2022); Kiss  _et al._ (2022); Stetcu  _et al._ (2022a). We associate each
qubit with a single-particle state in the configuration space, which can
either be empty (projection $0$) or occupied (projection $1$). Panel (b) of
Fig. 1 shows the mapping between single-particle states and qubits for the $p$
(bottom), $sd$ (central) and $pf$ shells (top panel). From a memory-storage
perspective, a shell-model VQE under the JW mapping only requires as many
qubits as single-particle states in the configuration space. In other words,
the number of qubits remains constant for all nuclei described within a given
shell. If a VQE can be used to diagonalize the problem and is robust against
errors, the approach may provide access to much larger configuration spaces,
currently unattainable in classical computers.

A VQE uses the Rayleigh-Ritz variational principle Ritz (1909); Rayleigh
(1870) to calculate the ground-state of a Hamiltonian starting from an initial
ansatz. Our algorithm of choice is ADAPT-VQE Grimsley  _et al._ (2019); Romero
_et al._ (2022); Tang  _et al._ (2021); Haidar  _et al._ (2022); Feniou  _et
al._ (2023), which iteratively builds a wavefunction of the form

$|\psi(\bm{\theta})\rangle=\prod_{k=1}^{n}e^{i\theta_{k}A_{k}}|{\rm
ref}\rangle,$ (5)

where $|{\rm ref}\rangle$ is an initial (reference) state of the quantum
system, $k$ is the iteration (or layer) index, $A_{k}$ are particle-hole
excitation operators, and $\bm{\theta}=\\{\theta_{i},i=1,\dots,n\\}$ are a set
of variational parameters. We stress that the adapted wavefunction in Eq. (5)
is free of Trotter-Suzuki approximation errors Barkoutsos  _et al._ (2018);
Childs  _et al._ (2021). This ansatz does not require decomposing an
exponential map of a sum of excitation operators, as would be the case in
algorithms such as UCC-VQE Anand  _et al._ (2022); Kiss  _et al._ (2022).

The minimization of the energy of this wavefunction with respect to the
parameters $\bm{\theta}$,

$E=\min_{\bm{\theta}}\frac{\langle\psi(\bm{\theta})|H_{\rm{eff}}|\psi(\bm{\theta})\rangle}{\langle\psi(\bm{\theta})|\psi(\bm{\theta})\rangle},$
(6)

can be performed classically Pellow-Jarman  _et al._ (2021) and yields an
approximate ground-state energy. Here, we use the BFGS optimiser with a
gradient tolerance set to $10^{-6}$ at every iteration. At each layer $k$ of
the iterative procedure, the ansatz grows by one parametrized unitary,
$|\psi(\bm{\theta})\rangle\to e^{i\theta_{k}A_{k}}|\psi(\bm{\theta})\rangle$.
The new operator $A_{k}$ is selected according to the largest energy gradient
computed as

$\left.\frac{\partial
E^{(n)}}{\partial\theta_{k}}\right|_{\theta_{k}=0}=\left.i\langle\psi(\bm{\theta})|[H_{\rm{eff}},A_{k}]|\psi(\bm{\theta})\rangle\right|_{\theta_{k}=0}.$
(7)

Thus, at every layer, the wavefunction adapts to the new information acquired
in the previous optimization. The set of parameters $\bm{\theta}$ are obtained
anew for every layer, so an updated state has no ties to former states. The
adaptive character of ADAPT-VQE should lead to implementations with shallower
circuits Grimsley  _et al._ (2019); Feniou  _et al._ (2023).

A crucial point for the optimal convergence towards the target state is the
choice of excitation operators $A_{k}$. These are predefined in an operator
pool, prior to the start of the simulation. Since our interest lies in the
nuclear shell model, with a Hamiltonian of the form in Eq. (1), we use a pool
of two-body fermionic excitation operators

$T_{rs}^{pq}=i(a_{p}^{{\dagger}}a_{q}^{{\dagger}}a_{r}a_{s}-a_{r}^{{\dagger}}a_{s}^{{\dagger}}a_{p}a_{q}),$
(8)

where $p,q,r$ and $s$ are single-particle labels with quantum numbers
$n,l,j,m$ and $t_{z}$. The same operator may be selected more than once
throughout the iterative process, but not on consecutive iterations. We apply
symmetry considerations when building the Slater determinant basis for the
nuclear ground state, and only consider excitation operators which conserve
the total angular momentum and isospin projection $M$ and $T_{z}$. This
iterative procedure continues until convergence, defined when all the gradient
norms in Eq. (7) vanish and/or when the energy is close enough to a known
solution from, for instance, classical diagonalization benchmarks. While one
could consider more complex operators, involving triple or quadruple particle-
hole excitations Stetcu  _et al._ (2022a); Kiss  _et al._ (2022), our
simulations indicate that, for the wide set of nuclei studied in this work,
full shell-model correlations can be captured at the two-body level with a
commensurate number of ansatz layers, of at most a few hundred.

### Circuit design strategy

The main aim of this paper is to determine the optimal architecture of quantum
circuits that can implement a nuclear shell-model VQE. We explore all the
necessary stages of a VQE, from the encoding to the energy measurement in the
Methods section. Ultimately, the circuit design strategy that we propose
provides

an approximation-free implementation of ADAPT-VQE, in a one-to-one
correspondence with the method Grimsley  _et al._ (2019); Tang  _et al._
(2021). Having access to the circuit structure across the full VQE
minimization process, including energy measurements, is a key step forward in
discussing the scalability of nuclear shell-model simulations in quantum
devices, and it is particularly critical to estimate the necessary resources
for nuclear shell-model simulations with a real quantum advantage, that is, in
isotopes or regions of the chart where current classical devices cannot be
employed.

We benchmark our circuit implementation with circuit-free ADAPT-VQE
simulations Romero  _et al._ (2022). The latter implement the full algorithm
using regular matrix calculus, expressing statevectors, Hamiltonians and pool
operators as sparse matrices in the Fock basis. With the circuit for the
ansatz built and optimized, we simulate the energy measurement protocol, to
test the circuits for the changes of basis needed to extract energies in an
actual quantum computer.

The state preparation protocol is the most resource-intensive part of the
algorithm and we provide indications of the resource costs in the Simulations
subsection. We can also quantify and optimize the scaling of the energy
measurements. The nuclear shell-model Hamiltonian in Eq. (1) consists of one
and two-body operators, which can be expressed in terms of Pauli strings (see
the Methods section). The one-body part of the Hamiltonian is diagonal and can
be measured directly. We divide the two-body part in three different kinds of
terms, depending on the number of repeated indices. Table 1 lists the number
of circuits needed to measure the expectation value of each part of the
Hamiltonian for the _p_ , _sd_ and _pf_ shells. Our design strategy indicates
that $100$ circuits should suffice to compute any isotope in the $p$ shell and
semi-magic nuclei in the $sd$ shell. Open-shell isotopes require a factor of
$4-6$ more circuits than their semi-magic counterparts in a given shell.

In a quantum computer implementation, an energy calculation will be affected
by statistical errors. Across a whole ADAPT-VQE simulation, the total number
of circuits to be measured for each layer will be the product of three terms,
$N_{s}\times N_{tot}\times N_{fc}$. The number of shots, $N_{s}$, is of
statistical nature and, as discussed in the Methods section in the context of
Eq. (19), it will be sensitive to error mitigation schemes. $N_{tot}$ is the
number of different energy measurement circuits. We estimate this number and
show the results in Table 1. Finally, $N_{fc}$ is the number of function calls
from the classical optimizer, which we analyze in the Supplementary
Information.

shell | $N_{qb}$ | $N_{\text{h}}$ | $N_{\text{hh}}$ | $N_{tot}$
---|---|---|---|---
_p_ | 6 | 2 | 10 (9) | 13 (12)
12 | 4 | 109 (44) | 114 (49)
_sd_ | 12 | 8 | 203 (86) | 212 (95)
24 | 16 | 1389 (518) | 1406 (535)
_pf_ | 20 | 20 | 1507 (570) | 1528 (591)
40 | 40 | 10572 (3459) | 10613 (3500)
Table 1: Number of circuits needed to measure the expectation value of the
nuclear shell-model for the _p_ , _sd_ and _pf_ shells. $N_{qb}$ indicates the
number of qubits for only neutrons or protons (top row for each shell) or both
nucleon types (bottom). $N_{\text{h}}$ and $N_{\text{hh}}$ are the number of
single- and double-hopping terms in the Hamiltonian (related to $h_{ijki}$ and
$h_{ijkl}$, respectively), defining the number of circuits needed to measure
these parts. The last column lists the total number of circuits,
$N_{\text{h}}+N_{\text{hh}}+1$, accounting also for the single circuit needed
to measure $\langle n_{i}\rangle$ and $\langle h_{ijij}^{(l)}\rangle$. The
values in parenthesis correspond to the minimum number of groups containing
$h_{ijkl}$ terms that commute with each other and thus can be measured with
the same circuit.

### Simulations

shell | $N_{qb}$ | $N_{\text{SD}}$ | nucleus | $N_{{\rm layers}}$ |  $\varepsilon_{E}$ bound |  $N_{{\rm C}}$ (bound)
---|---|---|---|---|---|---
_p_ | 6 | 5 |  6Be | 2 | $10^{-8}$ | 42 (80)
12 | 10 |  6Li | 9 | $10^{-7}$ | 92 (176)
53 |  8Be | 48 | $10^{-7}$ | 68 (176)
51 |  10Be | 48 | $10^{-7}$ | 62 (176)
21 |  13C | 19 | $10^{-7}$ | 77 (176)
_sd_ | 12 | 14 |  18O | 5 | $10^{-6}$ | 99 (176)
37 |  19O | 32 | $10^{-6}$ | 85 (176)
81 |  20O | 70 | $10^{-6}$ | 98 (176)
142 |  22O | 117 | $10^{-6}$ | 93 (176)
24 | 640 |  20Ne | 167 | $2\times 10^{-2}$ | 137 (368)
4206 |  22Ne | 236 | $2\times 10^{-2}$ | 137 (368)
7562 |  24Ne | 345 | $2\times 10^{-2}$ | 138 (368)
_pf_ | 20 | 30 |  42Ca | 9 | $10^{-8}$ | 116 (304)
565 |  44Ca | 132 | $10^{-2}$ | 153 (304)
3952 |  46Ca | 124 | $10^{-2}$ | 139 (304)
12022 |  48Ca | 101 | $10^{-2}$ | 137 (304)
17276 |  50Ca | 221 | $10^{-2}$ | 130 (304)
Table 2:  Ansatz and circuit depth for a given energy bound. Number of ansatz
layers ($N_{{\rm layers}}$) and relative-error ($\varepsilon_{E}$) upper
bounds for the ground-state energy of all nuclei simulated in this work,
organized according to their configuration space ($p$, $sd$, and $pf$ shells),
number of qubits $N_{qb}$, and of many-body configurations (Slater
determinants) $N_{SD}$. The last column reports the average number of CNOT
gates per layer $N_{\rm C}$ together with its upper bound, $16(N_{qb}-2)$ (see
Methods). For nuclei with $N_{{\rm layers}}>100$, the average only accounts
for the first $100$ layers.

The systems we explore include nuclei across different shells, with even and
odd numbers of protons and neutrons (see panel (a) of Fig. 1). We find that
circuit-free and circuit-full simulations employing the same parameter
minimization algorithm agree to numerical accuracy.

We estimate the required depth of a circuit by imposing bounds on the relative
error of the ground-state energy, $\displaystyle\varepsilon_{E}=\frac{\lvert
E-E_{\textrm{SM}}\lvert}{E_{\textrm{SM}}}$, where $E_{\textrm{SM}}$ is the
corresponding classical shell-model diagonalization result. Table 2 lists the
number of ADAPT-VQE layers needed in an ansatz state to achieve a given value
of $\varepsilon_{E}$ for a series of nuclei across the $p$, $sd$ and $pf$
shells. All energies tend to converge to the benchmark values, albeit with
different rates. Semi-magic nuclei close to the closed shell typically
converge rapidly, with less than $10$ ADAPT-VQE layers. In contrast, the most
costly nuclei simulated in this work, neon isotopes, require a few hundred
ADAPT-VQE layers to reach a ground-state energy error of $2\%$. Nonetheless,
we stress that the optimizations do not get stuck in barren plateaus. A key
advantage of our circuit design strategy is that it allows us to quantify the
associated quantum circuit resources. We take the number of CNOT gates
required in the state preparation, $N_{\text{CNOT}}$, as a quantitative
indicator of circuit resources.

Figure 2 shows the evolution of $\varepsilon_{E}$ (top panel) and
$N_{\text{CNOT}}$ (bottom) as a function of the number of ADAPT-VQE layers for
four representative isotopes across different nuclear shells. Simulations for
all nuclei show that $\varepsilon_{E}$ decreases exponentially as the number
of layers in the ansatz increases, while the number of CNOT gates grows
linearly or polynomially.

This number depends on the particular operators chosen by the ADAPT-VQE
minimization, but it is at most $16\,(N_{qb}-1)$ per ansatz layer (see Methods
section). In contrast, the average number of CNOT gates per ansatz layer found
by ADAPT-VQE simulations is

roughly half of the corresponding upper bounds, see Table 2. As an example,
finding the ground-state energy of 22O with an error of few percent, requires
about $20$ ansatz layers and $\approx 2000$ CNOT gates. We provide more
details for all the nuclei studied in this work in the Supplementary
Information.

Figure 2:  Energy relative error and circuit complexity as a function of
ADAPT-VQE layers. Evolution of the relative error for the ground-state energy,
$\varepsilon_{\rm E}$, (top panel) and number of CNOT gates in the ansatz
circuit (bottom) as a function of the number of ansatz layers for simulations
of 8Be, 22O, 20Ne and 48Ca. As the algorithm adaptively iterates, errors decay
exponentially while the number of CNOT gates increases linearly or
polynomially.

Figure 2 and Table 2 demonstrate that ADAPT-VQE converges exponentially as the
number of layers, or equivalently CNOT gates, is increased.

Our results are either commensurate or competitive compared to previous
estimates of circuit depth based on UCC-VQE on the $p$ shell and on two oxygen
isotopes on the $sd$ shell Stetcu  _et al._ (2022a); Kiss  _et al._ (2022).
For 8Be, Stetcu _et al._ require $112$ variational parameters to reach
$\varepsilon_{E}\approx 1\%$ even after including triple and quadruple
excitation operators Stetcu  _et al._ (2022a). Our implementation of ADAPT-
VQE, with two-body excitation operators only, requires $48$ parameters to
reach $\varepsilon_{E}=10^{-7}$. In 22O, the UCC-VQE ansatz leads to
$\varepsilon_{E}\approx 3\%$ with $35$ parameters Stetcu  _et al._ (2022a),
whereas Fig. 3 indicates that ADAPT-VQE reaches a similar level of accuracy
with about $20$ layers. For 6Li, we find that $9$ layers suffice to get a
converged result up to $10^{-7}$, in contrast to the observations of Ref. Kiss
_et al._ (2022), where an alternative ADAPT-VQE implementation reaches only
$\varepsilon_{E}\approx 10^{-3}$. A difference between previous
implementations and our work is that we let our classical minimizer reach
bottom precision at each ADAPT-VQE layer, whereas Kiss _et al._ employ $10$
minimization steps per layer (with the SPSA optimiser) Kiss  _et al._ (2022).
Moreover, UCC-VQE shell-model implementations have so far relied on Hartree-
Fock reference states, which may not be optimal starting points for VQEs
Romero  _et al._ (2022); Stetcu  _et al._ (2022b). Either way, it appears that
ADAPT-VQE shell-model simulations outperform their UCC-VQE counterparts in
terms of layers, an observation that is in line with findings in quantum
chemistry Haidar  _et al._ (2022). We note, however, that an unbiased
comparison of quantum hardware efficiency between different methods requires a
one-to-one quantification of the resources in each approach, including
explicitly energy measurement overheads.

Figure 3:  Quality of the wavefunction and entanglement entropy as a function
of ADAPT-VQE layers. Evolution of the relative error for the ground-state
energy, $\varepsilon_{\rm E}$, the infidelity, $I$, and the average relative
error of single-orbital entropies, $\varepsilon_{S^{(1)}}$ for 20O as a
function of the number of ansatz layers (top panel). Evolution of
$S_{i}^{(1)}$ for the same nucleus and $i$ orbitals $0d_{3/2}$, $1s_{1/2}$ and
$0d_{5/2}$, where the dotted lines indicate the entropies for the exact
solution (bottom panel). The maximum $S_{k}^{(1)}$ is $1$, very close to the
value of the $0d_{5/2}$ orbitals.

ADAPT-VQE predicts the ground-state energy of the nucleus, but one also has
access to the nuclear wavefunction $|\psi(\bm{\theta})\rangle$, although
reconstructing it from quantum hardware may require costly quantum tomography.
One can quantify the quality with respect to a given benchmark wavefunction,
$|\psi_{\text{b}}\rangle$, by employing the infidelity
$I=1-|\langle\psi_{\text{b}}|\psi(\bm{\theta})\rangle|^{2}$. We take the
classical shell model as a benchmark, and the better the level of agreement
between both wavefunctions, the closer $I$ is to $0$. We also use the single-
orbital entanglement entropy,
$S_{i}=-(1-\gamma_{i})\log_{2}(1-\gamma_{i})-\gamma_{i}\log_{2}\gamma_{i}$,
with $\gamma_{i}=\langle\psi(\bm{\theta})\lvert
a^{\dagger}_{i}a_{i}\lvert\psi(\bm{\theta})\rangle$ , bound between $0$ and
$1$, to evaluate the importance of quantum correlations in the ansatz Gigena
and Rossignoli (2015); Robin  _et al._ (2021); Johnson and Gorton (2023);
Bulgac  _et al._ (2023); Pazy (2023); Bulgac (2022).

These two indicators provide quantitative complementary information on the
quality of the wavefunction and the variational process.

Focusing on the test case example of 20O, the top panel of Fig. 3 shows the
infidelity $I$ of the ground state with respect to the shell-model
wavefunction (dashed line). The panel also shows the average of relative
errors of each single-particle state entanglement entropy,
$\varepsilon_{S^{(1)}}=\frac{1}{N_{qb}}\sum_{i}\varepsilon_{S_{i}^{(1)}}$
(dotted line). These two quantities follow closely $\varepsilon_{E}$ along the
iterative process. We observe a few sudden drops in the relative error for the
energy, which correlate with similar drops in $I$ and $\varepsilon_{S^{(1)}}$.
This indicates that, at certain points in the optimization, ADAPT-VQE
entangles parts of the nucleus relatively faster than others. Overall, the
curves suggest that the ADAPT-VQE ansatz captures efficiently the entanglement
structure of the many-body wavefunction. A more extensive analysis of the
infidelity is provided in the Supplementary Information. The bottom panel of
Fig. 3 provides a closer inspection to the entanglement structure of this
nucleus. Based on previous studies Robin  _et al._ (2021); Stetcu  _et al._
(2022a); Johnson and Gorton (2023), we expect nuclear-structure features to
correlate with single-particle states entanglement properties. The panel shows
the quantum simulated single-orbital entropies of the $12$ single-particle
states as a function of the number of ansatz layers, compared to the classical
shell-model entropies (horizontal dotted lines). We clearly distinguish the
emergence of three subshells in the entropy. The most entangled qubits are
those in the lowest-energy orbital, $0{d_{5/2}}$, reaching almost the maximal
value. These are followed by the $1{s_{1/2}}$ and the $0{d_{3/2}}$ states,
which are correspondingly less entangled (and occupied). The entropies
saturate to the shell-model value relatively quickly, within about 20 layers.
We take this as an indication that ADAPT-VQE captures early on the most
important correlations of the nucleus, which are subsequently refined by the
variational process.

## Discussion

In this work, we provide a detailed framework for a quantum hardware
implementation of ADAPT-VQE tailored to nuclear shell-model calculations. The
algorithm requires as many qubits as the number of single-particle states, a
relatively small number ($\approx 50$) even for valence spaces demanding
currently unavailable classical computational resources. We benchmark our
results with calculations using a circuit-free, regular matrix implementation
of the algorithm.

Our simulations do not become stuck in local minima or barren plateaus. We
find that the majority of the resources in the quantum circuit are dedicated
to the construction of the parametrized ansatz wave function. Each additional
parameter in the ansatz increases the circuit depth linearly with the number
of qubits. In contrast, the preparation of the reference state and the
implementation of the basis changes to measure Hamiltonian expectation values
are comparatively small parts of the total circuit depth. We quantify (see
Methods) the number of circuits needed to measure energies in the different
isotopes. Our proposed energy-measuring circuits are not substantially deeper
than the corresponding circuit encoding the wave function.

We calculate the ground state of selected nuclei in the _p-_ , _sd-_ and _pf-_
shell valence spaces, using up to 24 qubits. For all these systems, our
simulations indicate that the relative error in the ground-state energy and
the infidelity decrease exponentially as the number of layers in the ansatz
increases (see Supplementary Information). While the number of parameters
needed to reach a certain precision depends on the nucleus, our results
indicate that at most $150$ CNOT gates per ADAPT-VQE layer are necessary to
get ground-state energies accurate at the percent level. This suggests that a
circuit implementation of the shell model with ADAPT-VQE may be a suitable way
forward for quantum computing simulations of nuclei. Nevertheless, the number
of layers and CNOTs shown in Table 2 do not demonstrate an exponential quantum
advantage Lee  _et al._ (2023) with respect to the classical computation cost.
This is indeed seen more clearly in Fig. 4, which shows the number of total
CNOTs needed to obtain an energy relative error of $2\%$, as a function of the
number of Slater determinants for all nuclei studied in this work. Figure 4
indicates that up to nuclear masses $A\simeq 50$ the number of CNOT gates
scales roughly as the number of Slater determinants.

Figure 4:  Correlation between number of CNOTs and Slater determinants. Total
number of CNOTs $N_{\textrm{CNOT}}$ needed to obtain a ground-state relative
energy error of $2\%$ as a function of the number of Slater determinants
$N_{\textrm{SD}}$ in the many-body basis for all nuclei considered in this
work. The observed trend does not indicate a quantum exponential advantage
over classical methods.

Our study opens several potential avenues for further exploration. First,
different fermionic encodings may reduce the number of CNOT gates, which are
subject to noise errors that can limit realistic implementations in quantum
devices. A preliminary analysis using the Bravyi-Kitaev basis Seeley  _et al._
(2012) (instead of a JW transformation) suggests a $\approx 10\%$ reduction in
the number of CNOT gates of ADAPT-VQE after $100$ iterations in the $sd$ and
$pf$ shells. Other options of fermionic mappings such as Gray code encoding Di
Matteo  _et al._ (2021); Siwach and Arumugam (2021) should also be explored.
Second, the present work is an ideal testbed for the implementation of quantum
information tools for the study of nuclear structure. Our calculated single-
particle state entropies reveal the entanglement structure of nuclei, in close
analogy to the occupation probabilities of the orbitals obtained in classical
diagonalization schemes. Other correlation measures, such as quantum discord
Faba  _et al._ (2021a, b, 2022), will be the subject of future work.
Furthermore, one should elucidate more clearly the sharp differences between
the UCC and ADAPT ansatz VQEs. On the one hand, the choice of initial states,
at the mean-field level Kiss  _et al._ (2022); Stetcu  _et al._ (2022a) or
mixing many-body configurations, may improve the overall performance Romero
_et al._ (2022); Ruiz Guzman and Lacroix (2022) of the minimization process.
On the other, understanding why the ordering in the choice of operators is so
relevant may provide further insights into nuclear many-body correlations. A
better understanding on these issues is key to find optimal algorithms and
circuit designs for the nuclear shell model that avoid the exponential scaling
of resources and can be realistically implemented in NISQ devices.

We note that there are promising alternative algorithms for nuclear shell-
model calculations based on the Lanczos method Kirby  _et al._ (2023).

## Methods

We simulate circuits for several $p$-, $sd$\- and $pf$-shell nuclei using the
statevector simulator qibo Efthymiou  _et al._ (2021), together with the
qibojit package, which harnesses multi-core parallelization based on JIT
(just-in-time) compilation and the numba compiler Lam  _et al._ (2015). qibo
has been found to be specially efficient when compared to other simulators for
similar fermionic quantum-circuit simulations Pérez-Obiol  _et al._ (2022). At
each layer, we execute the quantum circuit to extract a statevector
$|\psi_{n}\rangle$ of dimension $2^{N_{qb}}$. This extraction is limited by
classical computer resources, which in turn provide stringent mass limits for
our classical circuit simulations. For instance, simulating open-shell nuclei
in the $pf$ shell valence space, requires state-vectors with $2^{40}$ complex
coefficients, demanding $8$ TB of memory in single-precision format. When
dealing with $20$ or more qubits, we use GPUs and the cupy compiler Okuta  _et
al._ (2017) to accelerate computations.

| Fermion Operators | Qubit Operators
---|---|---
$n_{p}$ | $a_{p}^{{\dagger}}a_{p}$ | $\dfrac{1}{2}(1-Z_{p})\rule{0.0pt}{15.0694pt}\rule[-9.04166pt]{0.0pt}{0.0pt}$
$h_{pqrs}$ | $\begin{aligned} a_{p}^{{\dagger}}a_{q}^{{\dagger}}a_{r}a_{s}\\\ +a_{r}^{{\dagger}}a_{s}^{{\dagger}}a_{p}a_{q}\end{aligned}$ | $\begin{aligned} \dfrac{1}{8}P_{rs}^{pq}\,(&-X_{p}X_{q}X_{r}X_{s}+X_{p}X_{q}Y_{r}Y_{s}\\\ &-X_{p}Y_{q}X_{r}Y_{s}-X_{p}Y_{q}Y_{r}X_{s}\\\ &-Y_{p}Y_{q}Y_{r}Y_{s}+Y_{p}Y_{q}X_{r}X_{s}\\\ &-Y_{p}X_{q}Y_{r}X_{s}-Y_{p}X_{q}X_{r}Y_{s})\end{aligned}\rule{0.0pt}{36.5971pt}$
$T_{rs}^{pq}$ | $\begin{aligned} i(a_{p}^{{\dagger}}a_{q}^{{\dagger}}a_{r}a_{s}\\\ -a_{r}^{{\dagger}}a_{s}^{{\dagger}}a_{p}a_{q})\end{aligned}$ | $\begin{aligned} \dfrac{1}{8}P_{rs}^{pq}\,(&-X_{p}Y_{q}Y_{r}Y_{s}-Y_{p}X_{q}Y_{r}Y_{s}\\\ &+Y_{p}Y_{q}X_{r}Y_{s}+Y_{p}Y_{q}Y_{r}X_{s}\\\ &+Y_{p}X_{q}X_{r}X_{s}+X_{p}Y_{q}X_{r}X_{s}\\\ &-X_{p}X_{q}Y_{r}X_{s}-X_{p}X_{q}X_{r}Y_{s})\end{aligned}\rule{0.0pt}{36.5971pt}$
$h_{pq}$ | $a_{p}^{{\dagger}}a_{q}+a_{q}^{{\dagger}}a_{p}$ | $\dfrac{1}{2}\left(\prod_{n=p+1}^{q-1}Z_{n}\right)\left(X_{p}X_{q}+Y_{q}Y_{p}\right)\rule{0.0pt}{15.0694pt}\rule[-9.04166pt]{0.0pt}{0.0pt}$
$T_{pq}$ | $i(a_{p}^{{\dagger}}a_{q}-a_{q}^{{\dagger}}a_{p})$ | $\dfrac{1}{2}\left(\prod_{n=p+1}^{q-1}Z_{n}\right)\left(Y_{p}X_{q}-X_{q}Y_{p}\right)\rule{0.0pt}{15.0694pt}\rule[-9.04166pt]{0.0pt}{0.0pt}$
Table 3: Jordan-Wigner transformation for the main operators appearing in the
Hamiltonian and in our ADAPT-VQE operator pool. Indices run over $p<q$ and
$r<s$, assuming that all are different. If two indices are repeated, then
$h_{pqpr}=-n_{p}h_{qr}$ and $T_{pq}^{pr}=n_{p}T_{qr}$, with $q<r$. We note
that $h_{pqpq}=-2n_{p}n_{q}$ and $T_{pq}^{pq}=0$.

Next, we describe the five different stages Tilly  _et al._ (2022) of our VQE
circuit design strategy.

### Mapping

We consider the JW mapping Jordan and Wigner (1993); Seeley  _et al._ (2012),
which transforms nucleonic creation and annihilation operators as

$a_{i}^{{\dagger}}=\left(\prod_{k=0}^{i-1}Z_{k}\right)\sigma^{-}_{i},~{}~{}a_{i}=\left(\prod_{k=0}^{i-1}Z_{k}\right)\sigma^{+}_{i},$
(9)

where $\sigma^{\pm}_{j}=\frac{1}{2}(X_{j}\pm iY_{j})$ and $X_{j}$, $Y_{j}$,
$Z_{j}$ are the usual Pauli matrices applied to qubit $j$. Using these
relations we can express any fermionic operator in terms of Pauli strings.
Table 3 lists the expressions for the two types of (self-adjoint) terms
appearing in the nuclear shell-model Hamiltonian $H_{\rm{eff}}$ in Eq. (1). We
use an auxiliary operator

$P_{rs}^{pq}\equiv\left(\prod_{m=p+1,m\notin[r,s]}^{q-1}Z_{m}\right)\left(\prod_{n=r+1,n\notin[p,q]}^{s-1}Z_{n}\right).$
(10)

Table 3 also indicates the JW transformation for the pool operators
$T_{rs}^{pq}$, and for single-excitation operators which appear when indices
are repeated in either $h_{pqrs}$ or $T_{rs}^{pq}$. In this context, the most
important features of an operator are the numbers and lengths of the Pauli
strings they contain. These ultimately determine the efficiency in the circuit
implementation of ADAPT-VQE. The two operators $h_{pqrs}$ and $T_{rs}^{pq}$
contain $8$ Pauli strings, each of length
$L_{pqrs}=n_{2}+n_{4}-n_{1}-n_{3}+2$, where $n_{1}$, $n_{2}$, $n_{3}$ and
$n_{4}$ are the indices $p$, $q$, $r$ and $s$ sorted in ascending order. For
example, if $(p,q,r,s)=(2,8,5,7)$, then $(n_{1},n_{2},n_{3},n_{4})=(2,5,7,8)$
and $L_{2857}=6$. If two indices are repeated, the expressions simplify to
$h_{pqpr}$ and $T_{pq}^{pr}$, as indicated in Table 3. These consist of two
Pauli strings of length $L_{pqr}^{(1)}=r-q+1$ and two other strings of length
$L_{pqr}^{(2)}=r-q+2$.

Figure 5: Examples of main circuit blocks, separated by dashed boxes, in
ADAPT-VQE for the simulation of 6Be. Left: preparation of the reference state
defined in Eq. (11) and Eq. (12). Middle: implementation of
$e^{-i\frac{\theta}{2}X_{2}X_{3}Y_{4}Z_{5}}$ using the CNOT staircase
algorithm, one out of the many unitaries in the variational part of ADAPT-VQE.
Right: circuit of the basis change $M_{0123}$ needed to diagonalize
$h_{0123}$. The subcircuit in qubits $q_{2}$ and $q_{3}$ containing two CNOTs
and a Hadamard gate $H$ corresponds to the basis change $M_{23}$.

### Initial state preparation

To provide a minimal starting point to the simulations, we choose the lowest-
energy Slater determinant as a reference state. Under the JW mapping, Slater
determinants are mapped to the computational basis by flipping the qubits
corresponding to the occupied orbitals using $X$ gates. Considering for
example the case of 6Be, an isotope in the $p$ shell (panel (b) of Fig. 1) and
for our interaction of choice, the lowest-energy Slater determinant is

$|0,3\rangle=a_{0}^{{\dagger}}a_{3}^{{\dagger}}\,|\rm{vac}\rangle,$ (11)

where $|\rm{vac}\rangle$ is the vacuum state with no particles in the valence
space. After a JW mapping, the state is translated into the computational
basis as

$|100100\rangle=X_{0}X_{3}|000000\rangle.$ (12)

The leftmost block of Fig. 5 shows the corresponding circuit.

This choice of initial state preparation is minimal in terms of circuit
resources: it has unit depth independently of the number of orbitals in the
valence space and it does not involve any two-qubit gates. For a given valence
neutron and proton number, $N_{\text{CI}}$ and $Z_{\text{CI}}$, finding the
lowest energy Slater determinant requires at most

$N_{\text{SD}}$ operations. This task can be performed relatively quickly in a
classical computer, and is a one-off pre-processing overhead that we do not
incorporate in the circuit resources discussed below.

Figure 6: Circuit to prepare the 18O ground state. $X$ gates prepare the
reference state and FSWAP gates change the basis so that pool-operator
exponentials act on adjacent qubits. Multiqubit gates in boxes are defined as
$U_{rs}^{pq}(\theta)\equiv e^{i\theta T_{rs}^{pq}}$ and
$\theta_{0}=-0.157263$, $\theta_{1}=-0.437238$, $\theta_{2}=0.604663$,
$\theta_{3}=0.214431$, $\theta_{4}=-0.785469$.

### Variational optimization

The variational ansatz is parametrized as in Eq. (5), with pool operators
$A_{k}=T_{rs}^{pq}$ given in Table 3 after the JW transformation. We convert
the pool operators $T_{rs}^{pq}$ to Pauli strings using the OpenFermion
package McClean  _et al._ (2020), and for the circuits for the unitaries
$e^{i\theta T_{rs}^{pq}}$ we follow the staircase algorithm of Fig. 5. In the
simulated circuits we only use single-qubit and CNOT gates.

All Pauli strings in these sums commute with each other, so each term in
$T_{rs}^{pq}$ can be exponentiated separately and there is no need for a
Trotter-Suzuki approximation. This results in the expression

$\begin{split}e^{i\theta
T_{rs}^{pq}}=&e^{-i\theta^{\prime}P_{rs}^{pq}X_{p}Y_{q}Y_{r}Y_{s}}e^{-i\theta^{\prime}P_{rs}^{pq}Y_{p}X_{q}Y_{r}Y_{s}}\\\
\times&e^{i\theta^{\prime}P_{rs}^{pq}Y_{p}Y_{q}X_{r}Y_{s}}e^{i\theta^{\prime}P_{rs}^{pq}Y_{p}Y_{q}Y_{r}X_{s}}\\\
\times&e^{i\theta^{\prime}P_{rs}^{pq}Y_{p}X_{q}X_{r}X_{s}}e^{i\theta^{\prime}P_{rs}^{pq}X_{p}Y_{q}X_{r}X_{s}}\\\
\times&e^{-i\theta^{\prime}P_{rs}^{pq}X_{p}X_{q}Y_{r}X_{s}}e^{-i\theta^{\prime}P_{rs}^{pq}X_{p}X_{q}X_{r}Y_{s}},\end{split}$
(13)

with $\theta^{\prime}=\theta/8$ and $P_{rs}^{pq}$ given in Eq. (10). The
exponential of a single Pauli string is particularly easy to implement with
the staircase algorithm Sawaya  _et al._ (2020). If the Pauli string contains
only $Z$ matrices, the circuit contains two cascades of CNOTs and a $Z$
rotation, $R_{z}(\theta)\equiv e^{-i\frac{\theta}{2}Z}$, with
$-\frac{\theta}{2}$ the coefficient multiplying the Pauli string. If the
product contains an $X$ or $Y$ matrix, we apply a basis change in the
corresponding qubit, namely $X=HZH$ and $Y=R_{x}^{\dagger}ZR_{x}$, where $H$
is the Hadamard gate and $R_{x}$ the rotation $e^{-i\frac{\pi}{4}X}$. Figure 5
(middle) illustrates the procedure for the example implementation of
$e^{-i\frac{\theta}{2}X_{2}X_{3}Y_{4}X_{5}}$. If $e^{i\theta T_{rs}^{pq}}$
acts on non-adjacent qubits, we implement a change of basis through fermionic
SWAP (FSWAP) gates, so that only CNOTs applied to contiguous qubits are
needed. The FSWAP exchanges states while maintaining the correct parity,

$\text{FSWAP}=1+a_{i}^{\dagger}a_{j}+a_{j}^{\dagger}a_{i}-a_{i}^{\dagger}a_{i}-a_{j}^{\dagger}a_{j}.$
(14)

Using the staircase protocol, each parametrized layer $e^{i\theta
T_{rs}^{pq}}$ requires $16\,(L_{pqrs}-1)$ CNOT gates, where $L_{pqrs}$ is the
average length of the Pauli strings in the operator. $L_{pqrs}$ is bounded by
the number of qubits $N_{qb}$, implying that the maximum number of CNOTs per
ansatz layer is $16\,(N_{qb}-1)$ and that the depth per layer grows linearly
with the number of single-particle states in the valence space. If qubits are
linearly connected in hardware and non-adjacent qubit states are brought
together with FSWAPs, the depth per layer has a total linear overhead. The
precise overhead size depends on how qubits are arranged and connected to each
other. However, it is bounded by $4(N_{qb}-4)$.

Let us provide an example illustrating the simplicity of the ADAPT-VQE circuit
implementation. Obtaining the ground-state energy of simple nuclei only
demands a few operators. As shown in Results, ADAPT-VQE simulations for 18O
converge to an energy accuracy better than $10^{-6}$ with a five-layer ansatz,
reading

$|\psi_{{}^{18}{\rm
O}}\rangle=e^{i\theta_{4}T^{05}_{23}}e^{i\theta_{3}T^{05}_{9\,10}}e^{i\theta_{2}T^{05}_{14}}e^{i\theta_{1}T^{05}_{67}}e^{i\theta_{0}T^{05}_{8\,11}}X_{0}X_{5}|0\rangle^{\otimes
12}.$

Figure 6 shows the full circuit assuming one-dimensional connectivity between
qubits, and gives the parameter values. Our algorithm includes the multiqubit
operators $e^{i\theta T_{rs}^{pq}}$ involving CNOT gates acting on non-
adjacent qubits when these are laid out in a one-dimensional array. We
manipulate these operators to include only local two-qubit gates through a
series of FSWAPs.

### Measurement

Once the ADAPT-VQE ansatz $|\psi_{n}\rangle$ is prepared in the quantum
circuit at a given layer $n$, we measure the energy with the expectation value
$\langle\psi_{n}|H_{\rm{eff}}|\psi_{n}\rangle$. To this end, we build a series
of circuits that implement a change of basis to diagonalize separately each
term of the Hamiltonian. The number of terms in the shell-model Hamiltonian
scales with the number of qubits as $O(N_{qb}^{4})$, but we find a much milder
scaling of the circuit number with $N_{qb}$.

One-body (number) operators $n_{i}$ are diagonal and can be measured directly,

$\langle\psi_{n}|n_{i}|\psi_{n}\rangle=\frac{1}{2}\langle\psi_{n}|1-Z_{i}|\psi_{n}\rangle=p_{1}^{(i)},$
(15)

where $p_{1}^{(i)}$, the probability of measuring $``1"$ in qubit $i$, can be
extracted by measuring multiple times that qubit. Since all one-body operators
commute with each other, we can measure all of them simultaneously. The two-
body part of the Hamiltonian $h_{ijkl}$ can be divided into three kinds of
terms depending on whether indices ($i$, $j$, $k$, $l$) are two, three, or
four different integers. Local terms $h_{ijij}$ are the product of two number
operators $n_{i}$ and $n_{j}$ and they can be measured simultaneously,

$\langle\psi_{n}|h_{ijij}|\psi_{n}\rangle=-2\langle\psi_{n}|n_{i}n_{j}|\psi_{n}\rangle=-2p_{11}^{(ij)},$
(16)

with $p_{11}^{(ij)}$ the probability to measure $``1"$ in qubits $i$ and $j$.
The non-diagonal parts of $h_{ijik}$ and $h_{ijkl}$ swap two states in the
subspaces of qubits ($i$, $j$, $k$) and ($i$, $j$, $k$, $l$), respectively.
These operators can be disentangled through series of CNOT gates and reduced
to an $X$ gate acting on a single qubit. The Pauli matrix $X$ is then
diagonalized with a Hadamard gate, $X=HZH$. In turn, we diagonalize $h_{ijik}$
and $h_{ijkl}$ using $M_{jk}\equiv\,CX_{kj}\,H_{k}\,CX_{kj}$ and
$M_{ijkl}\equiv\,CX_{ij}CX_{ki}CX_{lk}H_{l}CX_{lk}CX_{ki}CX_{ij}$, where
$CX_{ij}$ represents a CNOT gate with control qubit $i$ and target qubit $j$.
The right block of Fig. 5 illustrates the corresponding circuit
implementation. After diagonalization, assuming contiguous indices, the
expectation values read

$\begin{split}\langle\psi_{n}|h_{ijik}|\psi_{n}\rangle=&p_{101}^{(ijk)}-p_{110}^{(ijk)},\end{split}$
(17)

and

$\begin{split}\langle\psi_{n}|h_{ijkl}|\psi_{n}\rangle=&p_{1100}^{(ijkl)}-p_{0011}^{(ijkl)},\end{split}$
(18)

with $p_{r_{1}\cdots r_{k}}^{(q_{1}\dots q_{k})}$ being the probabilities of
measuring results $r_{1}$ to $r_{k}$ in qubits $q_{1}$ to $q_{k}$ in the
statevector where the basis changes have been applied. We refer to the
Supplementary Information for a detailed derivation of Eq. (17) and Eq. (18).

The changes of basis needed for measurements add, for any nucleus, an overhead
of zero, two or six two-qubit gates depending on the Hamiltonian term
measured. This represents a small fraction of the circuit depth and a constant
scaling with the number of single-particle states in the valence space. We
discuss in the Supplementary Information details regarding to the number of
different measurement circuits required to measure the energy as well as the
gradients of Eq. (7).

### Error mitigation

Finally, expectation values of the Hamiltonian computed using the algorithm
described above are subject to statistical errors and quantum noise. The
former scale as the inverse of the number of shots,
$\sigma_{E}\propto\frac{1}{\sqrt{N_{s}}}$. In other words, given a target
error in the energy accuracy $\varepsilon_{\langle H\rangle}$, the number of
necessary shots scales as

$N_{s}\propto\frac{1}{\varepsilon_{\langle H\rangle}^{2}}.$ (19)

The specific factor may be estimated simulating the measurement protocol. A
straightforward and robust strategy to mitigate errors for ADAPT-VQE shell-
model simulations is to use symmetry considerations and discard measurements
that do not yield results consistent with the Fock basis of the simulated
nucleus. Since the JW mapping identifies Fock and computational states, this
amounts to excluding all states with different number of measured $``1$”s than
nucleons in the valence space. Likewise, one should also ignore states with
measured $``1$”s distributed in a set of qubits corresponding to a different
angular momentum or isospin than the simulated nucleus. This protocol should
be particularly effective in mitigating single bit-flip errors, which
effectively create or destroy nucleons, as well as multiple bit-flip errors
which do not preserve either nucleon number, angular momentum or isospin.
These simple but robust strategies may be key in future implementations of
this method on NISQ devices.

## Data availability

The data that support the findings of this study are available within the
paper and its Supplementary Information. Any additional information is
available from the corresponding authors upon request.

## References

  * Taniuchi  _et al._ (2019) R. Taniuchi _et al._ , Nature 569, 53 (2019), arXiv:arxiv:1912.05978 [nucl-ex] .
  * Butler  _et al._ (2019) P. A. Butler _et al._ , Nat. Commun. 10, 2473 (2019), arXiv:arxiv:2003.10147 [nucl-ex] .
  * Tsunoda  _et al._ (2020) N. Tsunoda, T. Otsuka, K. Takayanagi, N. Shimizu, T. Suzuki, Y. Utsuno, S. Yoshida, and H. Ueno, Nature 587, 66 (2020).
  * Schmidt  _et al._ (2020) A. Schmidt _et al._ (CLAS), Nature 578, 540 (2020), arXiv:arxiv:2004.11221 [nucl-ex] .
  * Mukha  _et al._ (2006) I. Mukha _et al._ , Nature 439, 298 (2006).
  * Hinke  _et al._ (2012) C. B. Hinke _et al._ , Nature 486, 341 (2012).
  * Walz  _et al._ (2015) C. Walz, H. Scheit, N. Pietralla, T. Aumann, R. Lefol, and V. Y. Ponomarev, Nature 526, 406 (2015).
  * Cowan  _et al._ (2021) J. J. Cowan, C. Sneden, J. E. Lawler, A. Aprahamian, M. Wiescher, K. Langanke, G. Martínez-Pinedo, and F.-K. Thielemann, Rev. Mod. Phys. 93, 15002 (2021), arXiv:arxiv:1901.01410 [astro-ph.HE] .
  * Aalbers  _et al._ (2023) J. Aalbers _et al._ , J. Phys. G 50, 013001 (2023), arXiv:arxiv:2203.02309 [physics.ins-det] .
  * Engel  _et al._ (2013) J. Engel, M. J. Ramsey-Musolf, and U. van Kolck, Prog. Part. Nucl. Phys. 71, 21 (2013), arXiv:arxiv:1303.2371 [nucl-th] .
  * Avignone  _et al._ (2008) F. T. Avignone, III, S. R. Elliott, and J. Engel, Rev. Mod. Phys. 80, 481 (2008), arXiv:arxiv:0708.1033 [nucl-ex] .
  * Mayer (1949) M. G. Mayer, Phys. Rev. 75, 1969 (1949).
  * Haxel  _et al._ (1949) O. Haxel, J. H. D. Jensen, and H. E. Suess, Phys. Rev. 75, 1766 (1949).
  * Brown and Wildenthal (1988) B. A. Brown and B. Wildenthal, Annu. Rev. Nucl. Part. Sci. 38, 29 (1988).
  * Caurier  _et al._ (2005) E. Caurier, G. Martinez-Pinedo, F. Nowacki, A. Poves, and A. P. Zuker, Rev. Mod. Phys. 77, 427 (2005), arXiv:arxiv:nucl-th/0402046 .
  * Otsuka  _et al._ (2020) T. Otsuka, A. Gade, O. Sorlin, T. Suzuki, and Y. Utsuno, Rev. Mod. Phys. 92, 015002 (2020), arxiv:1805.06501 .
  * Stroberg  _et al._ (2019) S. R. Stroberg, S. K. Bogner, H. Hergert, and J. D. Holt, Annu. Rev. Nucl. Part. Sci. 69, 307 (2019), arXiv:arxiv:1902.06154 [nucl-th] .
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Nature 574, 505 (2019).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018), arxiv:1801.00862 .
  * Peruzzo  _et al._ (2014) 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, 1 (2014), arxiv:1304.3061 .
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. 18, 023023 (2016), arxiv:1509.04279 .
  * Bharti  _et al._ (2022) K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, _et al._ , Rev. Mod. Phys. 94, 015004 (2022), arxiv:2101.08448 .
  * Cerezo  _et al._ (2021) M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Nat. Rev. Phys. 3, 625 (2021), arxiv:2012.09265 .
  * Tilly  _et al._ (2022) J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, _et al._ , Phys. Rep. 986, 1 (2022), arxiv:2111.05176 .
  * Anand  _et al._ (2022) A. Anand, P. Schleich, S. Alperin-Lea, P. W. Jensen, S. Sim, M. Díaz-Tinoco, J. S. Kottmann, M. Degroote, A. F. Izmaylov, and A. Aspuru-Guzik, Chem. Soc. Rev. 51, 1659 (2022), arxiv:2109.15176 .
  * McArdle  _et al._ (2020) S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Rev. Mod. Phys. 92, 015003 (2020), arxiv:1808.10402 .
  * Haidar  _et al._ (2022) M. Haidar, M. J. Rančić, T. Ayral, Y. Maday, and J.-P. Piquemal, “Open Source Variational Quantum Eigensolver Extension of the Quantum Learning Machine (QLM) for Quantum Chemistry,”  (2022), https://arxiv.org/abs/2206.08798, arxiv:2206.08798 .
  * Cade  _et al._ (2020) C. Cade, L. Mineh, A. Montanaro, and S. Stanisic, Phys. Rev. B 102, 235122 (2020), arxiv:1912.06007 .
  * Cervera-Lierta (2018) A. Cervera-Lierta, Quantum 2, 114 (2018), arxiv:1807.07112 .
  * Cervia  _et al._ (2021) M. J. Cervia, A. Balantekin, S. Coppersmith, C. W. Johnson, P. J. Love, C. Poole, K. Robbins, and M. Saffman, Phys. Rev. C 104, 024305 (2021), arxiv:2011.04097 .
  * Harsha  _et al._ (2018) G. Harsha, T. Shiozaki, and G. E. Scuseria, J. Chem. Phys. 148, 044107 (2018), arxiv:1711.00579 .
  * Faba  _et al._ (2022) J. Faba, V. Martín, and L. Robledo, Phys. Rev. A 105, 062449 (2022), arxiv:2203.09400 .
  * Wahlen-Strothman  _et al._ (2017) J. M. Wahlen-Strothman, T. M. Henderson, M. R. Hermes, M. Degroote, Y. Qiu, J. Zhao, J. Dukelsky, and G. E. Scuseria, J. Chem. Phys. 146, 054110 (2017), arxiv:1611.06273 .
  * Robin and Savage (2023) C. E. P. Robin and M. J. Savage, “Quantum Simulations in Effective Model Spaces (I): Hamiltonian Learning-VQE using Digital Quantum Computers and Application to the Lipkin-Meshkov-Glick Model,”  (2023), https://arxiv.org/abs/2301.05976, arxiv:2301.05976 .
  * Lacroix (2020) D. Lacroix, Phys. Rev. Lett. 125, 230502 (2020), arxiv:2006.06491 .
  * Ruiz Guzman and Lacroix (2022) E. A. Ruiz Guzman and D. Lacroix, Phys. Rev. C 105, 024324 (2022), arxiv:2111.13080 .
  * Qian  _et al._ (2022) W. Qian, R. Basili, S. Pal, G. Luecke, and J. P. Vary, Phys. Rev. Res. 4, 043193 (2022), arXiv:arxiv:2112.01927 [quant-ph] .
  * Grimsley  _et al._ (2019) H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, Nat. Commun. 10, 1 (2019), arxiv:1812.11173 .
  * Sapova and Fedorov (2022) M. D. Sapova and A. K. Fedorov, Commun. Phys. 5, 199 (2022), arxiv:2108.11167 .
  * Feniou  _et al._ (2023) C. Feniou, M. Hassan, D. Traoré, E. Giner, Y. Maday, and J.-P. Piquemal, “Overlap-ADAPT-VQE: Practical Quantum Chemistry on Quantum Computers via Overlap-Guided Compact Ansätze,”  (2023), https://arxiv.org/abs/2301.10196, arxiv:2301.10196 .
  * Dumitrescu  _et al._ (2018) E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean, and P. Lougovski, Phys. Rev. Lett. 120, 210501 (2018), 1801.03897 .
  * Lu  _et al._ (2019) H.-H. Lu, N. Klco, J. M. Lukens, T. D. Morris, A. Bansal, A. Ekström, G. Hagen, T. Papenbrock, A. M. Weiner, M. J. Savage, and P. Lougovski, Phys. Rev. A 100, 012320 (2019), arxiv:1810.03959 .
  * Stetcu  _et al._ (2022a) I. Stetcu, A. Baroni, and J. Carlson, Phys. Rev. C 105, 064308 (2022a), arxiv:2110.06098 .
  * Kiss  _et al._ (2022) O. Kiss, M. Grossi, P. Lougovski, F. Sanchez, S. Vallecorsa, and T. Papenbrock, Phys. Rev. C 106, 034325 (2022), arxiv:2205.00864 .
  * Shalit and Talmi (1963) A. Shalit and I. Talmi,  _Nuclear Shell Theory_ (Academic Press N. Y., 1963).
  * Talmi (1993) I. Talmi,  _Simple Models of Complex Nuclei: The Shell Model and Interacting Boson Model_ , Beitrage Zur Wirtschaftsinformatik (Harwood Academic Publishers, 1993).
  * Varshalovich  _et al._ (1988) D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii,  _Quantum theory of angular momentum_ (World Scientific, 1988).
  * Hjorth-Jensen  _et al._ (1995) M. Hjorth-Jensen, T. T. Kuo, and E. Osnes, Phys. Rep. 261, 125 (1995).
  * Epelbaum  _et al._ (2009) E. Epelbaum, H.-W. Hammer, and U.-G. Meissner, Rev. Mod. Phys. 81, 1773 (2009), arXiv:arxiv:0811.1338 [nucl-th] .
  * Poves and Zuker (1981) A. Poves and A. Zuker, Phys. Rep. 70, 235 (1981).
  * Cohen and Kurath (1965) S. Cohen and D. Kurath, Nuclear Physics 73, 1 (1965).
  * Brown and Richter (2006) B. A. Brown and W. A. Richter, Phys. Rev. C 74, 034315 (2006).
  * Poves  _et al._ (2001) A. Poves, J. Sánchez-Solano, E. Caurier, and F. Nowacki, Nuclear Physics A 694, 157 (2001).
  * Caurier and Nowacki (1999) E. Caurier and F. Nowacki, Acta Physica Polonica 30, 705 (1999).
  * Shimizu  _et al._ (2019) N. Shimizu, T. Mizusaki, Y. Utsuno, and Y. Tsunoda, Comput. Phys. Commun. 244, 372 (2019), arXiv:arxiv:1902.02064 [nucl-th] .
  * Brown and Rae (2014) B. Brown and W. Rae, Nuclear Data Sheets 120, 115 (2014).
  * Johnson  _et al._ (2018) C. W. Johnson, W. E. Ormand, K. S. McElvain, and H. Shan, “BIGSTICK: A flexible configuration-interaction shell-model code,”  (2018), https://arxiv.org/abs/1801.08432, arxiv:1801.08432 .
  * Seeley  _et al._ (2012) J. T. Seeley, M. J. Richard, and P. J. Love, J. Chem. Phys. 137, 224109 (2012), arXiv:1208.5986 .
  * Romero  _et al._ (2022) A. M. Romero, J. Engel, H. L. Tang, and S. E. Economou, Phys. Rev. C 105, 064317 (2022), arxiv:2203.01619 .
  * Ritz (1909) W. Ritz, Journal für die reine und angewandte Mathematik 135, 1 (1909).
  * Rayleigh (1870) J. Rayleigh, Phil. Trans 161, 16 (1870).
  * Tang  _et al._ (2021) H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, PRX Quantum 2, 020310 (2021), arxiv:1911.10205 .
  * Barkoutsos  _et al._ (2018) P. K. Barkoutsos, J. F. Gonthier, I. Sokolov, N. Moll, G. Salis, A. Fuhrer, M. Ganzhorn, D. J. Egger, M. Troyer, A. Mezzacapo, _et al._ , Phys. Rev. A 98, 022322 (2018), arxiv:1805.04340 .
  * Childs  _et al._ (2021) A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Phys. Rev. X 11, 011020 (2021), arxiv:1912.08854 .
  * Pellow-Jarman  _et al._ (2021) A. Pellow-Jarman, I. Sinayskiy, A. Pillay, and F. Petruccione, Quantum Inf. Process. 20, 1 (2021), arxiv:2106.08682 .
  * Stetcu  _et al._ (2022b) I. Stetcu, A. Baroni, and J. Carlson, “Projection algorithm for state preparation on quantum computers,”  (2022b), https://arxiv.org/abs/2211.10545, arxiv:2211.10545 .
  * Gigena and Rossignoli (2015) N. Gigena and R. Rossignoli, Phys. Rev. A 92, 042326 (2015), arxiv:1509.05970 .
  * Robin  _et al._ (2021) C. Robin, M. J. Savage, and N. Pillet, Phys. Rev. C 103, 034325 (2021), arxiv:2007.09157 .
  * Johnson and Gorton (2023) C. W. Johnson and O. C. Gorton, Journal of Physics G: Nuclear and Particle Physics 50, 045110 (2023).
  * Bulgac  _et al._ (2023) A. Bulgac, M. Kafker, and I. Abdurrahman, Phys. Rev. C 107, 044318 (2023).
  * Pazy (2023) E. Pazy, Phys. Rev. C 107, 054308 (2023).
  * Bulgac (2022) A. Bulgac, “Entanglement entropy, single-particle occupation probabilities, and short-range correlations,”  (2022), https://arxiv.org/abs/2203.12079, arxiv:2203.12079 .
  * Lee  _et al._ (2023) S. Lee, J. Lee, H. Zhai, Y. Tong, A. M. Dalzell, A. Kumar, P. Helms, J. Gray, Z.-H. Cui, W. Liu, _et al._ , Nature Communications 14, 1952 (2023).
  * Di Matteo  _et al._ (2021) O. Di Matteo, A. McCoy, P. Gysbers, T. Miyagi, R. Woloshyn, and P. Navrátil, Phys. Rev. A 103, 042405 (2021), arxiv:2008.05012 .
  * Siwach and Arumugam (2021) P. Siwach and P. Arumugam, Phys. Rev. C 104, 034301 (2021).
  * Faba  _et al._ (2021a) J. Faba, V. Martín, and L. Robledo, Phys. Rev. A 104, 032428 (2021a), arxiv:2106.15993 .
  * Faba  _et al._ (2021b) J. Faba, V. Martín, and L. Robledo, Phys. Rev. A 103, 032426 (2021b), arxiv:2012.01292 .
  * Kirby  _et al._ (2023) W. Kirby, M. Motta, and A. Mezzacapo, Quantum 7, 1018 (2023).
  * Efthymiou  _et al._ (2021) S. Efthymiou, S. Ramos-Calderer, C. Bravo-Prieto, A. Pérez-Salinas, D. García-Martín, A. Garcia-Saez, J. I. Latorre, and S. Carrazza, Quantum Sci. Technol. 7, 015018 (2021), arxiv:2009.01845 .
  * Lam  _et al._ (2015) S. K. Lam, A. Pitrou, and S. Seibert, in  _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_, LLVM ’15 (Association for Computing Machinery, New York, NY, USA, 2015).
  * Pérez-Obiol  _et al._ (2022) A. Pérez-Obiol, A. Pérez-Salinas, S. Sánchez-Ramírez, B. G. M. Araújo, and A. Garcia-Saez, Phys. Rev. A 106, 052408 (2022), arxiv:2204.03013 .
  * Okuta  _et al._ (2017) R. Okuta, Y. Unno, D. Nishino, S. Hido, and C. Loomis, in  _Proceedings of Workshop on ML Systems in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS)_ (2017).
  * Jordan and Wigner (1993) P. Jordan and E. P. Wigner, in  _The Collected Works of Eugene Paul Wigner_ (Springer, 1993) pp. 109–129.
  * McClean  _et al._ (2020) J. R. McClean, N. C. Rubin, K. J. Sung, I. D. Kivlichan, X. Bonet-Monroig, Y. Cao, C. Dai, E. S. Fried, C. Gidney, B. Gimby, _et al._ , Quantum Sci. Technol. 5, 034014 (2020), arxiv:1710.07629 .
  * Sawaya  _et al._ (2020) N. P. D. Sawaya, T. Menke, T. H. Kyaw, S. Johri, A. Aspuru-Guzik, and G. G. Guerreschi, npj Quantum Inf. 6, 49 (2020), arxiv:1909.12847 .
  * Carbone  _et al._ (2013) A. Carbone, A. Cipollone, C. Barbieri, A. Rios, and A. Polls, Phys. Rev. C 88, 054326 (2013).
  * Hebeler (2021) K. Hebeler, Physics Reports 890, 1 (2021), arXiv:arxiv:2002.09548 [nucl-th] .

## Acknowledgements

A.M.R. thanks J. Engel for the support and fruitful discussions in the
conception of this project. A. P-O. and A. G-S. thank the QUANTIC group at BSC
for insightful comments and discussions along the realization of this work. A.
G-S. received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement No 951911 (AI4Media). This work is
financially supported by the Ministry of Economic Affairs and Digital
Transformation of the Spanish Government through the QUANTUM ENIA project call
- Quantum Spain project, by the European Union through the Recovery,
Transformation and Resilience Plan - NextGenerationEU within the framework of
the Digital Spain 2026 Agenda, by grants PID2020-118758GB-I00 and
PID2020-114626GB-I00 funded by MCIN/AEI/10.13039/501100011033; by the ”Ramón y
Cajal” grants RYC-2017-22781 and RYC2018-026072 funded by MCIN/AEI
/10.13039/501100011033 and FSE “El FSE invierte en tu futuro”; and by the
“Unit of Excellence María de Maeztu 2020-2023” award to the Institute of
Cosmos Sciences, Grant CEX2019-000918-M funded by
MCIN/AEI/10.13039/501100011033.

## Author contributions statement

A.P.O. and A.M.R. designed and performed research; A.P.O. analyzed data; and
A.P.O., A.M.R., J.M., A.R., A.G.S., and B.J.D. wrote the paper.

## Additional information

Supplementary information The online version contains supplementary material.

Competing interests The authors declare no competing interests.

Correspondence and requests for materials should be addressed to AMR or APO.

## Supplementary information

## Circuit design strategy

### Number of different measurement circuits

Here we discuss the number of different measurement circuits that are
necessary to compute expectation values of the energy as well as of the
products of operators required in the gradient calculations. Local terms
$n_{i}$ and $h_{ijij}$ can be measured simultaneously. We analyze and optimize
the number of different circuits needed to measure the expectation value of
the non-local part of $H_{\rm{eff}}$, $h_{ijki}$ and $h_{ijkl}$, for the $p$,
$sd$ and $pf$ shell valence spaces.

All terms $h_{ijki}=-n_{i}(a_{j}^{\dagger}a_{k}+a_{k}^{\dagger}a_{j})$ with
the same hopping (same indices j, k) and different local terms $n_{i}$ can be
measured simultaneously since they commute,
$[h_{ijik},h_{i^{\prime}ji^{\prime}k}]=0$. The local part of $h_{ijik}$
conserves the third components of the angular momentum and isospin, $m$ and
$t_{z}$, implying that the complementary hopping term involves only indices in
the same vertical axis in the panel (b) diagram of Fig. 1 in the main text.
For example, considering the $sd$ shell with only neutrons, this amounts to a
total of eight terms:

$\\{(j,k)\\}=\\{(1,8),(2,6),(6,9),(2,9),(3,7),(7,10),(3,10),(4,11)\\}.$ (20)

The number of different circuits needed to measure all $h_{ijki}$ terms is
then equivalent to the number of different $m$\- and $t_{z}$-conserving
single-excitation operators in the shell. This scales, in the worst case, as
$O(N_{qb}^{2})$, representing a relatively small number of circuits. Each term
is diagonalized with the circuit $M_{jk}=CX_{kj}H_{k}CX_{kj}$, which for
continguous indices, $k=j+1$, results in the operator $|101\rangle\langle
101|-|110\rangle\langle 110|$, where the indices $(i,j,k)$ have been omitted.
Therefore, $\langle h_{ijik}\rangle=p_{101}^{(ijk)}-p_{110}^{(ijk)}$, with
$p_{101}^{(ijk)}$ and $p_{110}^{(ijk)}$ the probabilities of measuring $101$
and $110$ in qubits $(i,j,k)$ after the change of basis.

The double-hopping terms $h_{ijkl}$ that involve different sets of orbitals
($i$, $j$, $k$, $l$) also commute and can be measured with the same circuit.
Given a group of self-commuting terms, products of $Z$s of one or more terms
$h_{ijkl}$ appearing in the JW mapping may overlap with the indices of another
term $h_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ in the group. A product of
an even number of overlapping $Zs$, for example
$P_{\textrm{even}}=Z_{i^{\prime}}Z_{j^{\prime}}$, commutes with
$M_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ and the same circuit $M_{ijkl}$
can be used for both. If there is a product of an odd number of overlapping
$Z$s, $P_{\textrm{odd}}$, then
$[P_{odd},M_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}]\neq 0$ and all the
different $h_{ijkl}$ operators need to be diagonalized simultaneously. Some
terms that share two indices also commute, but for simplicity we do not group
them into the same measurement.

###  .1 Simultaneous diagonalization of double-hopping terms with different
indices

Measuring the expected value of the Hamiltonian requires then a simultaneous
diagonalization of each term $h_{ijkl}$ with different values for the indices
($i$, $j$, $k$, $l$). These operators consist of the product
$h_{ijkl}=P_{ij}^{kl}O_{ijkl}$, where $P_{ij}^{kl}$ is a diagonal Pauli string
containing only $Z$s and $O_{ijkl}$ is the non-diagonal part,

$\begin{split}O_{ijkl}\equiv\,&\left(\sigma^{-}_{i}\sigma^{-}_{j}\sigma^{+}_{k}\sigma^{+}_{l}+\sigma^{-}_{k}\sigma^{-}_{l}\sigma^{+}_{i}\sigma^{+}_{j}\right)\\\
=&\frac{1}{8}(X_{i}X_{j}X_{k}X_{l}-X_{i}X_{j}Y_{k}Y_{l}+X_{i}Y_{j}X_{k}Y_{l}\\\
&+X_{i}Y_{j}Y_{k}X_{l}+Y_{i}Y_{j}Y_{k}Y_{l}-Y_{i}Y_{j}X_{k}X_{l}\\\
&+Y_{i}X_{j}Y_{k}X_{l}+Y_{i}X_{j}X_{k}Y_{l}),\\\ =&|0011\rangle\langle
1100|+|1100\rangle\langle 0011|,\end{split}$ (21)

where in the last line the indices ($i,j,k,l$) have been omitted, see Table 1
in the main text. To diagonalize a single term $h_{ijkl}$ we use the change of
basis $M_{ijkl}\equiv\,CX_{ij}CX_{ki}CX_{lk}H_{l}CX_{lk}CX_{ki}CX_{ij}$, such
that

$\begin{split}M_{ijkl}^{\dagger}O_{ijkl}M_{ijkl}\equiv\,&D_{ijkl}=|1100\rangle\langle
1100|-|0011\rangle\langle 0011|.\end{split}$ (22)

For contiguous indices, $j=i+1$, $l=k+1$, then $P_{ij}^{kl}=1$, and we have
$\langle h_{ijkl}\rangle=p_{1100}^{(ijkl)}-p_{0011}^{(ijkl)}$, dependent on
the probabilities of measuring 1100 and 0011 in qubits ($i$, $j$, $k$, $l$)
after applying the change of basis, as stated in Eq. (18) in the main text. In
the general case, $j>i+1$, $l>k+1$, and $P_{ij}^{kl}\neq 1$, the expected
value needs to account for the product of $Z$ matrices. For example,
considering $\langle Z_{q}\rangle=p_{0}^{(q)}-p_{1}^{(q)}$ and $\langle
Z_{q}Z_{r}\rangle=p_{00}^{(qr)}-p_{01}^{(qr)}-p_{10}^{(qr)}+p_{11}^{(qr)}$,

$\begin{split}\langle
Z_{q}O_{ijkl}\rangle=&[p_{01100}^{(qijkl)}-p_{00011}^{(qijkl)}]-[p_{11100}^{(qijkl)}-p_{10011}^{(qijkl)}]\\\
\langle
Z_{q}Z_{r}O_{ijkl}\rangle=&[p_{001100}^{(qrijkl)}-p_{000011}^{(qrijkl)}]-[p_{011100}^{(qrijkl)}-p_{010011}^{(qrijkl)}]-[p_{101100}^{(qrijkl)}-p_{100011}^{(qrijkl)}]+[p_{111100}^{(qrijkl)}-p_{110011}^{(qrijkl)}].\end{split}$
(23)

In the case where two terms $h_{ijkl}$,
$h_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ are simultaneously
diagonalized, the indices from the product of $Z$ matrices in each term might
overlap. If there is an even number of overlapping $Z$ matrices, $P_{ij}^{kl}$
commutes with $M_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ and the same
circuit $M_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ to diagonalize
$h_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ can be used, since
$P_{ij}^{kl}$ can be factored out. The same holds for $M_{ijkl}$. For example,
if there are two overlapping $Z$s,

$\begin{split}&M_{ijkl}^{\dagger}M_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}^{\dagger}\left(Z_{i^{\prime}}Z_{j^{\prime}}O_{ijkl}\right)\left(Z_{i}Z_{j}O_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}\right)M_{ijkl}M_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}\\\
&=\left(M_{ijkl}^{\dagger}O_{ijkl}M_{ijkl}\right)\left(M_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}^{\dagger}O_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}M_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}\right)Z_{i}Z_{j}Z_{i^{\prime}}Z_{j^{\prime}}\\\
&=D_{ijkl}D_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}Z_{i}Z_{j}Z_{i^{\prime}}Z_{j^{\prime}},\end{split}$
(24)

with $D_{ijkl}$ and $D_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ the
corresponding diagonal operators. If $P_{ij}^{kl}$ contains a product of three
$Z$s overlapping with ($i^{\prime}$, $j^{\prime}$, $k^{\prime}$,
$l^{\prime}$), then two can be factored out so that the problem is reduced to
simultaneously diagonalizing operators $Z_{l^{\prime}}O_{ijkl}$ and
$Z_{l}O_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$.

In practice, we only need to build new circuits that diagonalize a $2$-qubit
subspace, instead of the full $8$-qubit space. The non-diagonal part
$O_{ijkl}$ exchanges the states $|0011\rangle$ and $|1100\rangle$, effectively
operating in this two-state subspace through an $X$ gate. The circuit in the
right dashed box of Fig. 2 in the main text can be interpreted as a three-step
protocol. First, a change of basis through a set of CNOT gates such that $X$
operates only in the last qubit; second, a Hadamard gate acting on that qubit
to diagonalize $X$, $HXH=Z$; and third, the inverse sequence of CNOTs to
switch back to the original basis. If one term has an overlapping $Z$, then
instead of the Hadamard gate acting separately on each $4$-qubit circuit, we
need to diagonalize the corresponding $2$-qubit space. For example, if we want
to measure $Z_{i^{\prime}}O_{ijkl}$ and
$Z_{i}O_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ with the same circuit, we
need to diagonalize $X_{l}Z_{l^{\prime}}$ and $Z_{l}X_{l^{\prime}}$, and embed
the corresponding circuit, $CX_{ll^{\prime}}H_{l^{\prime}}CX_{ll^{\prime}}$,
within the change of basis, see Fig. 7.

###  .2 Circuits to diagonalize products of Hamiltonian and pool operators

In order to measure gradients using Eq. (7) in the main text, we need to
compute expected values of $h_{ijkl}T_{pq}^{rs}$. Similarly to $O_{ijkl}$,
this operator effectively swaps two states in the computational basis,

$\begin{split}h_{ijkl}T_{pq}^{rs}&=i|11001100\rangle\langle
00110011|-i|00110011\rangle\langle 11001100|\\\ &+i|00111100\rangle\langle
11000011|-i|11000011\rangle\langle 00111100|,\end{split}$ (25)

where we have assumed $P_{ij}^{kl}=P_{pq}^{rs}=1$. This operator can be
disentangled through a series of CNOT gates up to the 2-qubit operator
$X_{i}Y_{p}$, which is then diagonalized with the basis change
$CX_{ip}{R_{x}}_{i}CX_{ip}$. Figure. 8 illustrates the full circuit to
diagonalize $h_{ijkl}T_{pq}^{rs}$.

Figure 7: Quantum circuit to implement the change of basis to diagonalize
$Z_{l^{\prime}}O_{ijkl}Z_{l}O_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$ for
double-hopping terms. Figure 8: Quantum circuit $M_{ijkl}^{pqrs}$ to
diagonalize $h_{ijkl}T_{pq}^{rs}$ when all eight indices are different. The
corresponding expectation value,
$\langle\psi_{n}|h_{ijkl}T_{pq}^{rs}|\psi_{n}\rangle=-p_{00100010}+p_{00101010}-p_{10100010}+p_{10101010}$,
depends on $p_{m}$, the probabilities of measuring $m$ in the corresponding
qubits ($i$, $j$, $k$, $l$, $p$, $q$, $r$, $s$) in the statevector
$M_{ijkl}^{pqrs}|\psi_{n}\rangle$.

## Discussion on the complete simulation set

We choose the Cohen-Kurath interaction Cohen and Kurath (1965) in the $p$
shell, USDB Brown and Richter (2006) in the $sd$ shell and KB3G in the $pf$
shell Poves  _et al._ (2001). Explicit three-nucleon interactions are
typically neglected because their leading effects can be written as an
effective two-body term Carbone  _et al._ (2013); Hebeler (2021).

Figure 9 shows the dependence on the number of ansatz layers of the energy
error $\varepsilon_{\rm E}$ (top panels), infidelities $I$ (second-row
panels), number of CNOTs $N_{\textrm{CNOT}}$ (third-row panels) and number of
cost-function calls used by the classical optimizer $N_{\textrm{fc}}$ (bottom
panels) for all nuclei considered in this work. The iterative evolution shown
by Fig. 9 presents similar features to Figs. 4 and 5 of the main text, where
results are shown only for selected nuclei.

The first column of Fig. 9 indicates that all nuclei in the $p$ shell are
relatively straightforward to implement. They all converge quickly, reaching a
relative ground-state energy error $\varepsilon_{\rm E}<10^{-3}$ with only a
dozen layers. Only 8Be and 10Be, with 2 valence protons and 2 and 4 valence
neutrons, respectively, require circuit architectures with $\approx$ 50 layers
in order to capture their open-shell correlations, converging to a precision
below $\varepsilon_{\rm E}=10^{-6}$. We show $N_{\textrm{CNOT}}$ and
$N_{\textrm{fc}}$ only up to this point, since this is the accuracy threshold
of the classical minimizer. For all cases, the number of CNOT gates increases
smoothly with numbers between 65 and 85 gates per layer. Thus, the
implementation of $p$-shell nuclei in quantum circuits is promising in terms
of both width (number of qubits, $N_{qb}=12$ in this case) and depth (number
of total CNOTs).

Using a single Slater determinant as a reference state is usually enough for
the adaptive iterative procedure to reach the ground-state energy and
wavefunction exponentially by increasing the number of parameters. In some
cases, for particularly correlated systems, the initial state may be closer in
structure to an excited state than the ground state, and one may land into the
local minimum corresponding to the excited state. The only such situation we
encountered is 6Li, where a simple change of reference state was sufficient to
converge into the ground state. We also note that 6Be is represented in the
figure, but it converges in only $2$ layers.

The second column of Fig. 9 shows results for oxygen isotopes (with no valence
protons) and neon (two valence protons) in the $sd$ shell, studied with
circuits of $N_{qb}=12$ and $24$ qubits, respectively. We observe a stark
difference in the simulation of both isotopic chains: the adaptive procedure
—starting from a single Slater determinant reference state— needs
significantly more layers to capture the many-body correlations present in
open-shell neon isotopes. This is due to the relatively large many-body basis
dimension of these neon isotopes $\dim_{\rm{mb}}\approx 10^{4}-10^{5}$ (see
the right panel of Fig. 1 in the main text). Nevertheless, the number of CNOT
gates scales at most polynomially with the number of layers, with between $90$
and $100$ gates per layer for oxygen and between $110$ and $150$ for neon
isotopes. This relatively mild non-exponential scaling is promising toward the
implementation of ADAPT-VQE in NISQ devices. The bottom panel shows that the
number of calls to the cost function used by the classical optimizer at a
given iteration is similar for $sd$\- and $p$-shell nuclei. This suggests that
there is no bottleneck in resources associated to the classical optimizer.

Finally, the third column of Fig. 9 presents the results for calcium isotopes
(with no valence protons) in the $pf$ shell, using circuits with
$N_{\textrm{qb}}=20$ qubits. The first isotope, 42Ca, convergences extremely
quickly, within $10$ layers. In contrast, calcium isotopes with more than $2$
valence neutrons result in a slow convergence, similar to the one for neon
isotopes. Again, these calcium isotopes have $\dim_{\rm{mb}}\approx
10^{4}-10^{5}$, and the algorithm needs more updates of the wavefunction to
capture the strong correlations in their ground states. We find the slowest
convergence for 44Ca, a midshell isotope between the closed-shell 40Ca and
48Ca. Likewise, the infidelity of 44Ca seems to stall around $I\approx 3\times
10^{-2}$ and even the number of CNOT gates per layer grows beyond the range
found for the rest of isotopes. This suggests that a different choice of
reference state, involving more many-body basis states, may be required for a
faster convergence and, as a result, a reduction in quantum resources. In
contrast, we find again that the number of cost-function calls for all calcium
isotopes follows a similar trend to the $p$\- and $sd$-shell nuclei.

Figure 9: Evolution of the relative error for the ground-state energy,
$\varepsilon_{\rm E}$, (top row), infidelity $I$ (second row), number of CNOT
gates in the ansatz circuit $N_{\rm CNOT}$ (third row) and number of cost-
function calls $N_{\rm fc}$ in the classical optimizer (bottom row) as a
function of the number of ansatz layers for simulations of all $p$-shell
(first column), $sd$-shell (second column) and $pf$-shell (third column)
nuclei considered in this work. The bands in the number of CNOT gates panels
are meant to guide the eye and correspond to lower (upper) limits of CNOT
gates per layer of 65 (85) in $p$-shell nuclei, 90 (100) in oxygen isotopes
and 110 (150) in both neon and calcium isotopes. The number of CNOT gates
increases polynomially even in least favorable cases of convergence of 44Ca
and 24Ne. The relative energy error and infidelities follow analogous trends
during the iterative process. This indicates that the algorithm captures the
correlations in the nuclear wavefunctions. The number of calls to the cost-
function for the classical optimization presents a similar trend for all
nuclei, mildly increasing on average with the number of layers.