Shell-model study of $^{58}$Ni using quantum computing algorithm

This study presents a simulated quantum computing approach for the investigation into the shell-model energy levels of $^{58}$Ni through the application of the variational eigensolver (VQE) method in combination with a problem-specific ansatz. The primary objective is to achieve a fully accurate low-lying energy spectrum of $^{58}$Ni. The chosen isotope, $^{58}$Ni is particularly interesting in nuclear physics through its role in astrophysical reactions while also being a simple but not-trivial nucleus for shell-model study, it being two particles outside a closed shell. Our ansatz, along with the VQE method are shown to be able to reproduce exact energy values for the ground state and first and second excited states. We compare a classical shell model code, the values obtained by diagonalization of the Hamiltonian after qubit mapping, and a noiseless simulated ansatz+VQE simulation. The exact agreement between classical and qubit-mapped diagonalisation shows the correctness of our method, and the high accuracy of the simulation means that the ansatz is suitable to allow a full reconstruction of the full nuclear wave function.


I. INTRODUCTION
Atomic nuclei are systems of interacting spin-1/2 fermions.As such, their simulation is a strong candidate for a possibly revolutionary treatment using quantum computation, as has been explored in nuclei [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] and other manyfermion systems [17][18][19].A standard microscopic (treating nucleons as the degree of freedom) method to interpret observed states of nuclei is the configuration interaction nuclear shell model [20][21][22] in which nucleons interact in a defined single-particle basis to produce correlated eigenstates of a model Hamiltonian.Finding the wavefunctions and consequent properties of these eigenstates is the job of a shell model practitioner; a job which has already been explored on real and simulated quantum computers in the case of the lightest nuclei [23,24] with various aspects of the specific algorithms needed for quantum computation examined [25,26].
If quantum computation is to form a future tool for large-scale shell model calculations on near-term hardware, much work will be needed to understand appropriate algorithms to find the eigenstates and to deal with larger problem sizes.Here we take a step in the latter direction by exploring a heavier system than previously explored, looking at a nucleus ( 58 Ni) with valence particles in the f p-shell.In our exploration, a simulated variational quantum eigensolver (VQE) algorithm is used to compare aspects of the method, in comparison with an exact classical shell model calculation.The choice of 58 Ni is primarily motivated as being the simplest two-valence-particle system in the f p shell, and one with which we can compare our quantum simulation calculations with the classical shell model.This isotope is also of particular interest in astrophysics in the s-process in AGB stars [27][28][29], as a well-studied nuclear for use as a nuclear data benchmark [30], and as a nuclear material [31].
The primary objective of our research is to achieve a high level of precision in determining the low-lying energy levels of 58 Ni, by constructing a problem-based ansatz able to reproduce exact results based on the shell-model interaction JUN45 [32].We give a survey of the shell model and quantum computation methods in section II, before presenting results in section III and some concluding remarks in section IV.

II. THEORETICAL FRAMEWORK
In this study, we give details of our methods to investigate the properties of the 58 Ni nucleus using quantum computing simulation techniques, along with a classical shell model comparison.We take the model space comprising orbitals 1p 3/2 , 0f 5/2 , and 1p 1/2 above the 56 Ni core.This section describes the Hamiltonian employed along with its adaptation for quantum computing simulation.We present the tailored ansatz as appropriate to reach the ground state of the system.Furthermore, we discuss the application of the Variational Quantum Eigensolver (VQE) algorithm and associated optimizers for energy minimization.The implementation of excitation operators, crucial for constructing quantum circuits, is detailed.
The shell model Hamiltonian can be expressed in terms of second quantization, where creation and annihilation operators act on states representing single-particle wavefunctions.The shell model Hamiltonian H consists of the single-particle energy term and the two-body interaction term: Here, â † i and âi are the creation and annihilation operators, ϵ i is the single-particle energy, V ijkl represents the twobody matrix elements (TBMEs) describing interactions between nucleons in states |i⟩, |j⟩, |k⟩, and |l⟩.The nucleon state |i⟩ is characterized by quantum numbers: radial, orbital angular momentum, total spin, and projections of spin and isospin as |n, l, j, m α , m tα ⟩.We use the JUN45 [32] interaction, which consists of the 1p 3/2 , 0f 5/2 , 1p 1/2 and 0g 9/2 orbital in the usual spectroscopic notation.In this work, we have neglected the 0g 9/2 orbital, choosing only f p model space for a simpler approach comprising 12 single-particle states, which remain appropriate for the lowest-energy states.Since we consider the case of 58 Ni, consisting of a core plus two valence neutrons, we only consider the neutron single-particle states.We represent these orbitals with 12 qubits, one for each orbital, as detailed in Table I.
Since matrix elements of shell model interactions are generally given in the J-scheme in which single-particle states are coupled to good total angular momentum J, while a quantum computation implementation relies on a mapping of the uncoupled single-particle states (the "M -scheme") to each qubit, a transformation from the J-scheme to the M -scheme is applied.The transformation is performed using [22] The normalization constants N ab (JT ) and N cd (JT ) are defined as follows: Here, we follow the notation of Suhonen [22], in which a,b,c,d are the single-particle orbitals |a⟩ = |n α , l α , j α ⟩, |b⟩ = |n β , l β , j β ⟩, |c⟩ = |n γ , l γ , j γ ⟩, |d⟩ = |n δ , l δ , j δ ⟩ carrying no projection (m, m t ) labels.The Greek letters α, β, γ, δ are the nucleon states as a complete set of quantum numbers containing the isospin parts |α⟩ = |n α , l α , j α , m α , m tα ⟩ etc.In (2), the four quantities inside round brackets are the Clebsch-Gordan coefficient of angular momenta and isospins, ⟨ab; JT |V |cd; JT ⟩ is the J-scheme TBMEs.The quantum numbers represented by Greek letter indices are tabulated in Table I, and the m-scheme TBMEs (v αβγδ ) resulting from equation ( 2) are given in Appendix A, as transformed from the J-scheme values tabulated in [32].TABLE I: Qubit mapping for the valence space of 58 Ni.Each qubit corresponds to a specific combination of j the total angular momentum, m α projection on the z axis, and m tα the third component of the isospin, is fixed at 1/2 for all qubits to represent the neutron.
The Jordan-Wigner transformation [33] maps fermionic creation and annihilation operators onto Pauli spin matrices.Since the states of the qubits of quantum computers can be written as spinors, with Pauli spin matrices forming a complete set with which to define Hermitian operators acting on the qubits, this JW mapping is often used to translate general many-fermion Hamiltonians to a qubit representation.We follow this method, while noting that other mapping schemes exist [34].
The Jordan-Wigner transformation is defined as follows.Consider a system of N fermionic modes (e.g., nucleons in different energy levels in a nucleus).Each mode can be occupied (state |1⟩) or unoccupied (state |0⟩).The occupation number operators are denoted as nj = a † j a j , where a † j and a j are the creation and annihilation operators for fermions in mode j, respectively.The Jordan-Wigner transformation maps these fermionic operators to spin-1 2 (qubit) operators: where σ + j = 1 2 (σ x j + iσ y j ) and σ − j = 1 2 (σ x j − iσ y j ) are the spin raising and lowering operators for the j-th qubit and σ x j and σ y j are the Pauli X and Y operator, respectively, and σ z j is the Pauli Z operator for the j-th qubit.The product j−1 k=1 σ z k ensures that the fermionic anti-commutation relations are preserved.The Jordan-Wigner transformation and similar methods, like the Bravyi-Kitaev transformation [35], allow a successful mapping of fermionic Hamiltonians onto qubit degrees of freedom.Each possible transformation has its strengths and weaknesses in terms of resulting Hamiltonian complexity, circuit depth and locality/non-locality of representation.We choose the JW method as a standard against which future comparisons can be made.

C. VQE and the wavefunction ansatz
The Variational Quantum Eigensolver (VQE) [36][37][38] is a hybrid quantum-classical algorithm proposed as a promising approach for finding ground states in many-body quantum systems.It has been applied in different guises in nuclear physics problems [3,6,10,12,16,39] thanks to its ability to be implemented, in some cases, on current quantum computing hardware.The basis of the VQE algorithm is to have a parameterized wave function ansatz on a quantum computer, in which the parameters are free variables; the expectation value of the Hamiltonian in its transformed Pauli matrix form is evaluated on the quantum computer, then a classical algorithm adjusts the parameters of the wave function ansatz to seek a minimum of the Hamiltonian expectation value.This is done through a repeated cycle of expectation value measurements on the quantum computer and directed parameter adjustments on the side of the classical algorithm.An appropriate ansatz (a guess for the wave function) must be chosen.In nuclear physics, this could be inspired by traditional nuclear models or techniques, like the Unitary Coupled Cluster (UCC) method [40,41], or by problem-specific considerations.
In our study, we choose a problem-inspired ansatz in which Pauli X gates initially create the correct number of particles in a state with quantum numbers appropriate to the desired target state.For the ground state of 58 Ni, characterized by J = 0 (and hence M = 0), the initial state takes the form of a product state, such as |0, 1⟩ ≡ |000000000011⟩ -i.e. the states |n = 1, l = 1, j = 3/2, m = −3/2⟩ and |n = 1, l = 1, j = 3/2, m = 3/2⟩ are occupied.This state is straightforward to prepare as each qubit can be addressed individually and the default initial |0⟩ state of a qubit can be flipped to |1⟩ by application of an X gate.
Having created a simple single configuration reference state for a given target state, parameterized operators are applied to explore the other configurations expected to contribute to the given M -valued wave function.We use parameterized Givens rotation gates [42] to access selected configurations in 58 Ni with correct M values.These gates are particle number preserving, keeping us in the correct Fock space, and are suitably parameterized with enough freedom to find the lowest energy configuration.The excitations in the quantum simulations have been made using the single and double Givens excitation gates.
Single excitations, shown in Figure 1, involve the application of creation and annihilation operators to excite a single particle from one orbital to another, as represented by a † i a j .Since the nucleon carries an intrinsic angular momentum (spin), the single excitation process is applied so that it conserves m α by choosing the individual spins of the nucleons involved in the excitation.
Double excitations (Figure 2), on the other hand, simultaneously promote two neutrons from their respective orbitals to other orbitals, characterized by expressions like a † i a † j a k a l .The set of double excitation processes is chosen so that the total angular momentum projection along the z-axis remains constant by accounting for the combined spins of the two nucleons undergoing excitation.
The prepared ansatz for the ground state, first excited state, and second excited state is represented in Figure 3.Note that while all qubits are active in the ground state wave function, only three are active in the J = 4 state as only they can give a total projection M = 4, while nine are active in the J = 2 state.
In the model space for 58 Ni, only a limited number of combinations are feasible for the second excited state (4 + ).In this context, we have formulated the ansatz employing solely a single excitation, as illustrated in Figure 3b.The system is initialized in the state |1, 5⟩, explicitly applying the single excitation {1, 7} among the qubits.

D. Classical optimizers
Classical optimizers play a key role in hybrid quantum-classical algorithms, such as the Variational Quantum Eigensolver (VQE).The role of the classical optimizer is to adjust the parameters in the wave function ansatz -in our case the angles in each Givens rotation -to optimize the objective function -in our case minimization of the energy expectation value.Here is an overview of the different classical optimizers used in this work, which we take the built-in version from the IBM qiskit package [44].q q q q q q q q q q q 10 q 11 FIG. 3: The quantum circuit (a) ground state, (b) second-excited state, and (c) first-excited state used for the VQE simulations for 58 Ni.As mentioned in the main text, the circuits used are composed of a set of initial state preparation gates and single and double-excitation gates.Numbers on gates indicate qubits between which gates are acting.
1. COBYLA The COBYLA (Constrained Optimization BY Linear Approximations) [45] optimizer is a derivative-free method suited for nonlinear optimization problems (which are often present in nuclear physics models), especially where gradients are unavailable or costly to compute.It employs linear approximations to construct a model of the objective function, optimizing within a dynamically adjusted trust region.COBYLA is effective for handling moderate-sized problems and constraints, making it a practical choice in various applications, but may be less efficient for large-scale optimization tasks or when gradient-based methods are applicable.

SLSQP
The SLSQP (Sequential Least Squares Programming) [46] optimizer is a gradient-based method effective for solving optimization problems with both linear and nonlinear constraints.It operates by iteratively solving a series of quadratic subproblems, using gradient information to refine solutions at each step.SLSQP is particularly valuable in scenarios where precision in parameter control and cohesion to constraints are crucial, as in nuclear many-body problems.While efficient for problems with available and accurate gradient information, its performance can be limited in cases where such data is difficult to compute or unreliable.

SPSA
The SPSA (Simultaneous Perturbation Stochastic Approximation) [47] optimizer efficiently tackles high-dimensional optimization by using stochastic methods to approximate gradients, perturbing all parameters simultaneously.This reduces computational demands, especially in problems with many variables.SPSA excels in complex, noisy scenarios, although it might need more iterations and careful tuning to converge.Its suitability for complex nuclear many-body problems comes from its ability to manage numerous variables and work effectively without exact gradient information.

Gradient Descent
Gradient Descent is a basic gradient-based optimization method that minimizes functions by iteratively adjusting parameters against the gradient of the objective functions.While effective for problems with well-defined gradients, its performance may suffer with local minima, non-convexities, or noisy gradients.Commonly used in machine learning, its applicability in nuclear many-body problems is limited due to difficulties in obtaining precise gradients and navigating complex energy landscapes.Despite these challenges, it can be effective if accurate gradients are available.

III. RESULTS
This section presents the results obtained from the previously prepared parameterized quantum circuit in Sec.II C. We employ the four optimization algorithms given in the previous section: COBYLA, SLSQP, SPSA, and Gradient Descent.We investigate the performance and convergence behaviour of the different optimization algorithms in VQE simulations for computing the expectation value of 58 Ni.The optimization algorithms each navigate the parameter space adeptly, at least in the noiseless simulations we present here, avoiding barren plateaus [48] at the level of noiseless simulation.Table II shows a detailed breakdown of the properties associated with each circuit.Through the use of our problem-specific ansatz, the number of gates is kept low compared with more general ansatzes [25].We have not imposed nearest-neighbour qubit interaction which is a requirement for some quantum computing hardware, and for the J = 2 and J = 4 states, doing this would increase the gate count, though this increase could be minimized with careful re-mapping of single-particle states to qubits.Details of the performance of each optimizer are given in the following section.

A. Ground state
Ground state calculations are performed using the ground state ansatz given and motivated previously.In Figure 4, we compare the different optimizer results discussed in Sec.II D for the ground state.As illustrated in the figure, ( (  ( our prepared ansatz converges with all four optimizers.Notably, COBYLA converges fastest within 80 iterations, achieving an accuracy of 10 −3 MeV shown in Table III.The SLSQP also shows a similar convergence rate.This efficient convergence may be attributed to the problem-based ansatz preparation, which mitigates the presence of secondary minima.However, the other two optimizers are less advantageous for the prepared ansatz, requiring over 2000 iterations to converge to the desired state with a less accurate expectation value.The comparison between the VQE wavefunction and the classical shell model wavefunction is shown in Figure 5.The wavefunction is reproduced nicely with the VQE calculation.As 58 Ni is above magic nuclei 56 Ni (Z = 28, N =28), it is expected to exhibit single-particle nature as seen in the wavefunction for the ground state.Since the maximal occupied Qubit states represent the same orbital (in this case 1p 3/2 with > 80% probability), this shows the single-particle nature of that wavefunction.

B. First excited state
In Figure 6, we compare the different optimizer results discussed in Sec.II D for the first excited state.Similar to the ground state, our prepared ansatz converges with all four optimizers.Again, COBYLA exhibits the fastest convergence, achieving convergence in 105 iterations with an accuracy of 10 −2 MeV, as detailed in Table III.The remaining three optimizers converge close to the exact value but fall short of achieving precise accuracy.Similar to the ground state, SPSA and Gradient Descent take longer (more iterations) to converge.The comparison between the VQE wavefunction and the Shell model wavefunction is shown in Figure 7.The VQE calculation with COBYLA wellreproduces the reference shell model wavefunction, while the other optimizers are reasonable but miss e.g. a correct distribution across the (3/2, −1/2, 5/2, 5/2) and (3/2, 1/2, 5/2, −1/2) configurations.In the m-scheme wavefunction, this first excited J = 2 state is dominated by the (1p 3/2 ) 2 and (1p 3/2 1p 1/2 ) components, with probabilities of 63.49% and 12.54%, respectively.These characteristics are mirrored in the VQE wavefunction.

C. Second excited state
In Figure 8, we compare the optimizer results discussed in Sec.II D for the second excited state.The prepared ansatz converges with all four optimizers to the exact reference value as detailed in Table III.Notably, COBYLA and SLSQP converge in approximately 10 iterations.The perturbation parameter in SPSA has been tuned to 0.01 to get better convergence since it was showing oscillatory behaviour at the default calibration [49].The comparison between the VQE wavefunction and the Shell model wavefunction is shown in Figure 9.The VQE calculation closely reproduces the dominant reference wavefunction, which is 1p 3/2 0f 5/2 .The nuclear shell model provides a method to calculate properties of atomic nuclei in which nucleons interact in a given model space of single-particle states.By mapping single particle states to qubits, and using a suitable method to prepare wavefunctions -in our case, tailored ansatz and VQE algorithm -we have been able to extend the simulated quantum computation of nuclei up to the 58 Ni isotope, finding exact results for the lowest 0 + , 2 + , and 4 + states.
Prospects for the method include a straightforward extension to other nuclei in the f p shell, which requires no new qubits but more complicated circuits; possible simplification of the circuits through compilation or approximation techniques; and extension to heavier nuclei where the increasing single particle model space involves a linear increase in qubit number, automatically accompanied by an exponential increase in the size of the Hilbert space which hinders classical shell model calculations.Such simulated quantum computer calculations as we have shown will form the check of real quantum computation for future calculations.
TABLE IV: The m-scheme TBMEs v αβγδ in MeV (Equation 2) of JUN45 interaction [32].There is one to one mapping between the m-scheme basis (α) and the qubits (Table I).The core is taken as Ni, and the single-particle energies are considered -9.8280, -8.7087 and -7.8388 MeV for the p 3/2 , f 5/2 and p 1/2 orbital, respectively.The TBMEs should be multiplied by a mass-scaling factor (A/58) −0.3 in J-scheme, for any nuclei of mass A.Here T = 1 and M T = m tα + m tβ , as we have calculated only for the same type of nuclei (nn and pp).

TABLE II :
Properties of ansatz used for the ground state, first and second excited state with the number of parameters, CNOT gates(2-qubit), RY gates(1-qubit), H gate(1-qubit), X gate(1-qubit) and circuit depth.

TABLE III :
Summary of the results for the ground state (g.s.), first excited state (1 st e.s.), and second excited state (2 nd e.s.) alongside the shell model value for 58 Ni in pf model space.The exact result, obtained with exact diagonalization of the Hamiltonian and the E U CC energy is obtained by minimizing the Hamiltonian using ansatz.