Nuclear shell-model simulation in digital quantum computers

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.

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 [1][2][3], strong short-range correlations between pairs of nucleons [4], or decay modes driven by the strong [5], weak [6] or electromagnetic [7] forces.Furthermore, nuclear decays are crucial to understand the origin of heavy elements in the universe [8], and experiments using nuclei aim to answer fundamental physics questions such as which is the nature of dark matter [9], why matter dominates over antimatter in the universe [10], or whether neutrinos are their own antiparticles [11].
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 [12,13].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 [14][15][16][17], 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 [18].In the current noisy intermediate-scale quantum (NISQ) device era [19], variational quantum eigensolvers (VQE) [20,21] are among the most successful algorithms [22] exploiting the benefits of quantum computing to deal with complex many-body problems in physics [23,24] and chemistry [25][26][27].
In general, a VQE implementation requires a series of  well-defined stages [24], 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 [41][42][43][44].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.

Nuclear shell model
The nuclear shell model [14][15][16][17] 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 [45].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 [46].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 [47].
The nuclear Hamiltonian in a given configuration space can be written as where ε i is the energy of the single-particle state i and vijkl = v ijkl −v ijlk are antisymmetrized two-body matrix elements.a i and a † i are fermionic annihilation and creation operators associated to each single-particle state, i.The matrix elements vijkl can be obtained [17,48] from an effective field theory of the underlying theory of the nuclear force, quantum chromodynamics [49].Here, instead, we use standard phenomenological Hamiltonians, with components adjusted to better reproduce key properties of selected nuclei [50].We choose the Cohen-Kurath interaction in the p shell [51], USDB in the sd shell [52] and KB3G in the pf shell [53].
A suitable many-body basis, also referred to as Fock space, for shell-model calculations is provided by the so-called M −scheme [46], 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, and nuclear wavefunctions and their corresponding energies are eigenvectors and eigenvalues of the Hamiltonian matrix in the basis of Slater determinants.The c α coefficients are obtained through diagonalization employing state-of-the-art nuclear shell-model codes [54][55][56][57] 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 singleparticle basis of a nuclear shell consisting of several orbitals nl j is 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 where N CI (Z 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 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 SD [15].However, the scaling in either dim mb or N 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 [14][15][16][17].

Variational algorithm
Here, we implement the nuclear shell model in a quantum computer following a standard Jordan-Wigner (JW) mapping [43,44,58,59].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 memorystorage 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 [60,61] to calculate the ground-state of a Hamiltonian starting from an initial ansatz.Our algorithm of choice is ADAPT-VQE [27,38,40,59,62], which iteratively builds a wavefunction of the form where |ref⟩ is an initial (reference) state of the quantum system, k is the iteration (or layer) index, A k are particlehole excitation operators, and θ = {θ i , i = 1, . . ., n} are a set of variational parameters.We stress that the adapted wavefunction in Eq. ( 5) is free of Trotter-Suzuki approximation errors [63,64].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 [25,44].The minimization of the energy of this wavefunction with respect to the parameters θ, can be performed classically [65] 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, The new operator A k is selected according to the largest energy gradient computed as Thus, at every layer, the wavefunction adapts to the new information acquired in the previous optimization.The set of parameters θ 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 [38,40].
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 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 [43,44], 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 [38,62].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 [59].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 resourceintensive 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 I 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 × N tot × N f c .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 I.Finally, N f c is the number of function calls from the classical optimizer, which we analyze in the Supplementary Information.

Simulations
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 en- , where E SM is the corresponding classical shell-model diagonalization result.clei 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 CNOT , as a quantitative indicator of circuit resources.
Figure 2 shows the evolution of ε E (top panel) and N 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 ε 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 II.As an example, finding the ground-state energy of 22 O with an error of few percent, requires about 20 ansatz layers and ≈ 2000 CNOT gates.We provide more details for all the nuclei studied in this work in the Supplementary Information.Figure 2  compared to previous estimates of circuit depth based on UCC-VQE on the p shell and on two oxygen isotopes on the sd shell [43,44].For 8 Be, Stetcu et al. require 112 variational parameters to reach ε E ≈ 1% even after including triple and quadruple excitation operators [43].Our implementation of ADAPT-VQE, with two-body excitation operators only, requires 48 parameters to reach ε E = 10 −7 .In 22 O, the UCC-VQE ansatz leads to ε E ≈ 3% with 35 parameters [43], whereas Fig. 3 indicates that ADAPT-VQE reaches a similar level of accuracy with about 20 layers.For 6 Li, we find that 9 layers suffice to get a converged result up to 10 −7 , in contrast to the observations of Ref. [44], where an alternative ADAPT-VQE implementation reaches only ε E ≈ 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) [44].Moreover, UCC-VQE shellmodel implementations have so far relied on Hartree-Fock reference states, which may not be optimal starting points for VQEs [59,66].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 chem-  istry [27].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.ADAPT-VQE predicts the ground-state energy of the nucleus, but one also has access to the nuclear wavefunction |ψ(θ)⟩, although reconstructing it from quantum hardware may require costly quantum tomography.One can quantify the quality with respect to a given benchmark wavefunction, |ψ b ⟩, by employing the infidelity I = 1 − |⟨ψ b |ψ(θ)⟩| 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, , bound between 0 and 1, to evaluate the importance of quantum correlations in the ansatz [67][68][69][70][71][72].
These two indicators provide quantitative complementary information on the quality of the wavefunction and the variational process.
Focusing on the test case example of 20 O, 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, ε S (1) = 1 i (dotted line).These two quantities follow closely ε 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 ε 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 [43,68,69], we expect nuclear-structure features to correlate with single-particle states entanglement properties.The panel shows the quantum simulated singleorbital 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, 0d 5/2 , reaching almost the maximal value.These are followed by the 1s 1/2 and the 0d 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 (≈ 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 II do not demonstrate an exponential quantum advantage [73] 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 ≃ 50 the number of CNOT gates scales roughly as the number of Slater determinants.
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 [58] (instead of a JW transformation) suggests a ≈ 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 [74,75] 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 [32,76,77], 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 meanfield level [43,44] or mixing many-body configurations, may improve the overall performance [36,59] 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 [78].

METHODS
We simulate circuits for several p-, sd-and pf -shell nuclei using the statevector simulator qibo [79], together with the qibojit package, which harnesses multi-core parallelization based on JIT (just-in-time) compilation and the numba compiler [80].qibo has been found to be specially efficient when compared to other simulators for similar fermionic quantum-circuit simulations [81].At each layer, we execute the quantum circuit to extract a statevector |ψ n ⟩ 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 [82] to accelerate computations.
Next, we describe the five different stages [24] of our VQE circuit design strategy.

Mapping
We consider the JW mapping [58,83], which transforms nucleonic creation and annihilation operators as where σ ± j = 1 2 (X j ± 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 III 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 hpqpr = −nphqr and T pr pq = npTqr, with q < r.We note that hpqpq = −2npnq and T pq pq = 0.
(self-adjoint) terms appearing in the nuclear shell-model Hamiltonian H eff in Eq. ( 1).We use an auxiliary operator Table III also indicates the JW transformation for the pool operators T pq rs , and for single-excitation operators which appear when indices are repeated in either h pqrs or T pq rs .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 pq rs 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 pr pq , as indicated in Table III.These consist of two Pauli strings of length L (1) pqr = r − q + 1 and two other strings of length L (2)

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 Examples of main circuit blocks, separated by dashed boxes, in ADAPT-VQE for the simulation of 6 Be.Left: preparation of the reference state defined in Eq. ( 11) and Eq. ( 12).Middle: implementation of e −i θ 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 M0123 needed to diagonalize h0123.The subcircuit in qubits q2 and q3 containing two CNOTs and a Hadamard gate H corresponds to the basis change M23.
the qubits corresponding to the occupied orbitals using X gates.Considering for example the case of 6 Be, an isotope in the p shell (panel (b) of Fig. 1) and for our interaction of choice, the lowest-energy Slater determinant is where |vac⟩ is the vacuum state with no particles in the valence space.After a JW mapping, the state is translated into the computational basis as 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 CI and Z CI , finding the lowest energy Slater determinant requires at most N SD operations.This task can be performed relatively quickly in a classical computer, and is a one-off preprocessing overhead that we do not incorporate in the circuit resources discussed below.

Variational optimization
The variational ansatz is parametrized as in Eq. ( 5), with pool operators A k = T pq rs given in Table III after the JW transformation.We convert the pool operators T pq rs to Pauli strings using the OpenFermion package [84], and for the circuits for the unitaries e iθT pq rs 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 pq rs can be exponentiated separately and there is no need for a Trotter-Suzuki approximation.This results in the expression e iθT pq rs =e −iθ with θ ′ = θ/8 and P pq rs given in Eq. ( 10).The exponential of a single Pauli string is particularly easy to implement with the staircase algorithm [85].If the Pauli string contains only Z matrices, the circuit contains two cascades of CNOTs and a Z rotation, R z (θ) ≡ e −i θ 2 Z , with − θ 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 ZR x , where H is the Hadamard gate and R x the rotation e −i π 4 X .Figure 5 (middle) illustrates the procedure for the example implementation of e −i θ 2 X2X3Y4X5 .If e iθT pq rs 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, Using the staircase protocol, each parametrized layer e iθT pq rs 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 18 O converge to an energy accuracy better than 10 −6 with a five-layer ansatz, reading 23 e iθ3T 05 9 10 e iθ2T 05 14 e iθ1T 05 67 e iθ0T 05 8 11 X 0 X 5 |0⟩ ⊗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θT pq rs involving CNOT gates acting on non-adjacent Circuit to prepare the 18 O 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 pq rs (θ) ≡ e iθT pq rs and θ0 = −0.157263,θ1 = −0.437238,θ2 = 0.604663, θ3 = 0.214431, θ4 = −0.785469.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 |ψ n ⟩ is prepared in the quantum circuit at a given layer n, we measure the energy with the expectation value ⟨ψ n |H eff |ψ n ⟩.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 4 qb ), 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, where p 1 , 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 twobody 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, 11 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 ≡ CX kj H k CX kj and M ijkl ≡ 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 and with p (q1...q k ) r1•••r 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 twoqubit gates depending on the Hamiltonian term measured.This represents a small fraction of the circuit depth and a constant scaling with the number of singleparticle 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, σ E ∝ 1 √ Ns .In other words, given a target error in the energy accuracy ε ⟨H⟩ , the number of necessary shots scales as 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.
measured with the same circuit.Given a group of selfcommuting terms, products of Zs of one or more terms h ijkl appearing in the JW mapping may overlap with the indices of another term h i ′ j ′ k ′ l ′ in the group.A product of an even number of overlapping Zs, for example P even = Z i ′ Z j ′ , commutes with M i ′ j ′ k ′ l ′ and the same circuit M ijkl can be used for both.If there is a product of an odd number of overlapping Zs, P odd , then [P odd , M i ′ j ′ k ′ l ′ ] ̸ = 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.

A. 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 kl ij O ijkl , where P kl ij is a diagonal Pauli string containing only Zs and O ijkl is the non-diagonal part, 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 For contiguous indices, j = i + 1, l = k + 1, then P kl ij = 1, and we have 0011 , 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 kl ij ̸ = 1, the expected value needs to account for the product of Z matrices.For example, considering ⟨Z q ⟩ = p (q) 0 −p (q) 1 and ⟨Z q Z r ⟩ = p In the case where two terms h ijkl , h i ′ j ′ k ′ l ′ 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 kl ij commutes with M i ′ j ′ k ′ l ′ and the same circuit M i ′ j ′ k ′ l ′ to diagonalize h i ′ j ′ k ′ l ′ can be used, since P kl ij can be factored out.The same holds for M ijkl .For example, if there are two overlapping Zs, with D ijkl and D i ′ j ′ k ′ l ′ the corresponding diagonal operators.If P kl ij contains a product of three Zs overlapping with (i ′ , j ′ , k ′ , l ′ ), then two can be factored out so that the problem is reduced to simultaneously diagonalizing operators 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⟩ and |1100⟩, 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 ′ O ijkl and Z i O i ′ j ′ k ′ l ′ with the same circuit, we need to diagonalize X l Z l ′ and Z l X l ′ , and embed the corresponding circuit, CX ll ′ H l ′ CX ll ′ , within the change of basis, see Fig. 7.

B. 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 rs pq .Similarly to O ijkl , this operator effectively swaps two states in the computational basis, where we have assumed P kl ij = P rs pq = 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 xi CX ip .Figure .8 illustrates the full circuit to diagonalize h ijkl T rs pq . qi . Quantum circuit M pqrs ijkl to diagonalize h ijkl T rs pq when all eight indices are different.The corresponding expectation value, ⟨ψn|h ijkl T rs pq |ψn⟩ = −p00100010 + p00101010 − p10100010 + p10101010, depends on pm, the probabilities of measuring m in the corresponding qubits (i, j, k, l, p, q, r, s) in the statevector M pqrs ijkl |ψn⟩.

DISCUSSION ON THE COMPLETE SIMULATION SET
We choose the Cohen-Kurath interaction [51] in the p shell, USDB [52] in the sd shell and KB3G in the pf shell [53].Explicit three-nucleon interactions are typically neglected because their leading effects can be written as an effective two-body term [86,87].
Figure 9 shows the dependence on the number of ansatz layers of the energy error ε E (top panels), infidelities I (second-row panels), number of CNOTs N CNOT (third-row panels) and number of cost-function calls used by the classical optimizer N fc (bottom panels) for all nuclei considered in this work.The iterative evolution shown by Fig. 9 presents similar features to Figs.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 groundstate energy error ε E < 10 −3 with only a dozen layers.Only 8 Be and 10 Be, with 2 valence protons and 2 and 4 valence neutrons, respectively, require circuit architectures with ≈ 50 layers in order to capture their open-shell correlations, converging to a precision below ε E = 10 −6 .We show N CNOT and N 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 6 Li, where a simple change of reference state was sufficient to converge into the ground state.We also note that 6 Be 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 mb ≈ 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 qb = 20 qubits.The first isotope, 42 Ca, 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 mb ≈ 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 44 Ca, a midshell isotope between the closed-shell 40 Ca and 48 Ca.Likewise, the infidelity of 44 Ca seems to stall around I ≈ 3 × 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.9. Evolution of the relative error for the ground-state energy, εE, (top row), infidelity I (second row), number of CNOT gates in the ansatz circuit NCNOT (third row) and number of cost-function calls N 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 44 Ca and 24 Ne.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.

11 FIG. 1 .
FIG. 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 mb , in the M -basis as a function of the number of active neutrons in the configuration space, NCI.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.

2 FIG. 3 .
FIG. 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, εE, the infidelity, I, and the average relative error of single-orbital entropies, ε S (1) for20 O as a function of the number of ansatz layers (top panel).Evolution of S (1) i 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(1) k is 1, very close to the value of the 0d 5/2 orbitals.
⟨Z q Z r O ijkl ⟩ =[p 4 and 5 of the main text, where results are shown only for selected nuclei.

TABLE I .
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 h and N 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 h + N hh + 1, accounting also for the single circuit needed to measure ⟨ni⟩ and ⟨h 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.
Table II lists the number of ADAPT-VQE layers needed in an ansatz state to achieve a given value of ε E for a series of nu-ijij ⟩.
and Table II demonstrate that ADAPT-VQE converges exponentially as the number of layers, or equivalently CNOT gates, is increased.Our results are either commensurate or competitive shell N qb NSD nucleus N layers εE bound NC (bound) TABLE II.Ansatz and circuit depth for a given energy bound.Number of ansatz layers (N layers ) and relative-error (ε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) NSD.The last column reports the average number of CNOT gates per layer NC together with its upper bound, 16(N qb − 2) (see Methods).For nuclei with N layers > 100, the average only accounts for the first 100 layers.