Towards heavy-mass ab initio nuclear structure: Open-shell Ca, Ni and Sn isotopes from Bogoliubov coupled-cluster theory

Recent developments in nuclear many-body theory enabled the description of open-shell medium-mass nuclei from first principles by exploiting the spontaneous breaking of symmetries within correlation expansion methods. Once combined with systematically improvable inter-nucleon interactions consistently derived from chiral effective field theory, modern ab initio nuclear structure calculations provide a powerful framework to deliver first-principle predictions accompanied with theoretical uncertainties. In this Letter, controlled ab initio Bogoliubov coupled cluster (BCC) calculations are performed for the first time, targeting the ground-state of all calcium, nickel and tin isotopes up to mass A ~ 180. While showing good agreement with available experimental data, the shell structure evolution in neutron-rich isotopes and the location of the neutron drip-lines are predicted. The BCC approach constitutes a key development towards reliable first-principles simulations of heavy-mass open-shell nuclei.

Following such a strategy, heavy closed-shell nuclei (see Refs. [29][30][31]) have become recently accessible, all the way to 208 Pb [32].At the same time, nuclei away from shell closures still provide a formal and computational challenge.One way to overcome this limitation relies on valence-space (VS) techniques that build VS interactions from a variety of expansion methods for closed-shell nuclei [7,20,33,34] before proceeding to an exact diagonalization within that VS via modern shell-model codes [35,36].In particular, the VS formulation based on the IMSRG has become the method of choice to describe medium-mass open-shell nuclei up to the lightest tin isotopes [29].Still, VS approaches eventually suffer from the exponential scaling of the valencespace diagonalization, which makes it extremely difficult to push them to yet heavier open-shell nuclei.Thus, there have been extensive efforts to complement VS formulations with expansion methods building upon symmetry-broken reference states grasping strong static correlations associated with nuclear superfluidity and/or quadrupolar deformation [11,25,37,38].In particular, breaking  (1) symmetry associated with particle-number conservation leads to employing so-called Bogoliubov reference states.While such a paradigm was initially implemented within SCGF theory [37,39,40], it was later employed to design Bogoliubov MBPT (BMBPT) [11,12,41,42], Bogoliubov CC (BCC) [43] and Bogoliubov IMSRG [44] [10,13,18,[45][46][47][48][49].In this context, the non-perturbative BCC theory is presently employed for the first time to address the heaviest open-shell nuclei ever computed ab initio.

Many-body framework
Particle-number-breaking many-body frameworks are formulated starting from a Bogoliubov reference state |Φ⟩ generated from the physical vacuum |0⟩ through where  denotes a complex normalization constant.The state |Φ⟩ is a vacuum for the set of quasi-particle operators {  ,  †  } generated from standard particle operators {  ,  †  } via a unitary Bogoliubov transformation maintaining standard fermionic anti-commutation rules [50].
The variationally optimal Bogoliubov transformation is determined by solving Hartree-Fock-Bogoliubov (HFB) mean-field equations.In this work, rotational invariance is enforced to make |Φ⟩ an eigenstate of the total angular momentum  2 and its projection   with eigenvalues  = 0 and  = 0, respectively.
where the set of tensors   1 ... 2 constitute the unknown cluster amplitudes.The central quantity of the many-body formalism is the (non-Hermitian) similarity-transformed grand potential Ω ≡ (  − Ω  ) c , where the lower index 'c' stipulates the connected character of the operator [51].Once cluster amplitudes have been determined by minimizing the BCC residual  ≡ ⟨Φ| Ω|Φ⟩, where  projects on the manifold of single, double, . . .quasi-particle excitations of ⟨Φ|, the ground-state grand potential is computed as Ω 0 ≡ ⟨Φ| Ω|Φ⟩.
While BCC provides an exact parametrization of |Ψ⟩,  must be truncated to make the problem numerically tractable.In this Letter, the BCC with singles and doubles (BCCSD) truncation is employed, i.e.,  BCCSD ≡  1 +  2 .The numerical implementation of BCCSD employs the quasi-linear form of the amplitude equation [43].Rotational invariance is enforced on the amplitude equations thanks to angularmomentum-coupling techniques using the AMC program [52].The iterative procedure to determine single and double amplitudes is initialized using first-order BMBPT wave functions and convergence is accelerated by employing the direct inversion of the iterative subspace (DIIS) approach [53].The BCC iterations are carried out until all single and double residuals fall below ‖‖ ≤ 10 −3 such that the ground-state energy and chemical potentials are stable at sub-keV level.
Since |Φ⟩ breaks U(1) symmetry, so does the approximate BCCSD wave function.Moreover, many-body correlations tend to induce a shift of the average particle number with respect to |Φ⟩ such that, e.g., ⟨Ψ BCCSD ||Ψ BCCSD ⟩ ≠ N 0 .The BCCSD wave function is thus further constrained to carry the correct average neutron (proton) number, 1 i.e.
⟨Φ| Ñ( Z)|Φ⟩ = N 0 (Z 0 ), by readjusting the chemical potentials during the BCC iterations leading to the redefinition [54] 1 Since  and  commute with  , the associated values in the BCC groundstate can, just as the energy, be evaluated using a projective approach.

Fig. 1.
Convergence of the ground-state energy as a function of  (3)   max using  max = 12 with an oscillator frequency of ℏΩ = 16 MeV.In all cases the "EM1.8/2.0"interaction is employed. (3)

Interaction and model-space details
Numerical simulations of even-even Ca, Ni and Sn isotopes employ the "EM1.8/2.0"nuclear Hamiltonian [55] derived within chiral effective field theory [56][57][58][59].While the three-nucleon (3N) interaction operator is fully accounted for at the HFB level, BCC equations are solved using the particle-number-conserving normal-ordered two-body (PNO2B) approximation [60].The nuclear Hamiltonian is represented using a spherical harmonic oscillator (HO) single-particle basis truncated according to the principal quantum number  max ≡ (2 + ) max of the last included state, where  denotes the radial quantum number and  the orbital angular momentum.To keep the number of 3N matrix elements tractable, an additional cut is performed by restricting three-body basis states to  1 +  2 +  3 ≤  (3)   max .In all calculations the model space is truncated according to  max = 12 and  (3)  max = 24.A fully quantitative prediction of properties of neutron-rich systems eventually require a proper treatment of continuum effects (see e.g.Refs.[61,62]) that are not well resolved in the truncated harmonic oscillator single-particle basis employed in this work.

Many-body uncertainties
Gauging the quality of ab initio results necessitates an assessment of both Hamiltonian and many-body uncertainties.While the former, due to working at finite truncation order in the chiral power counting, is not addressed in this work, related discussions can be found in, e.g., Refs.[63][64][65].The latter is schematically given by and is reported as a band in the figures below.The finite basis-size ( max ) error  FBS depends on the mass regime and varies among the studied isotopic chains: calcium (1%), nickel (1.5%) and tin (2%) at the employed  max = 12 value.Based on previous assessments in midmass closed-shell nuclei [66], the error  PNO2B induced by omitting the residual normal-ordered 3N interaction is taken to be at the 2% level.

Three-body basis truncation error
The error  3NB due to further truncating (via  (3)   max ) the 3N basis has been the limiting factor to extend ab initio calculations to nuclei beyond Ni isotopes.Evaluating this error in details and eventually keeping it small enough up to neutron-rich tin isotopes is presently made possible by extending the novel three-body format of Ref. [31] to the context of  (1)-breaking expansion methods such as BMBPT and BCC.
For each isotopic chain a representative (neutron-rich) nucleus was chosen: 60 Ca, 88 Ni and 140 Sn.The convergence depends on the nuclear mass.While 60 Ca is reasonably well converged at  (3)  max = 16, there are still sizable contributions from discarded three-body matrix elements in heavier systems.In 140 Sn, working with  (3)  max = 16 induces an error of the order of several tens of MeV, demonstrating the need for a significantly larger 3N basis when targeting nuclei with  ≳ 80 [31].For sufficiently high values of  (3)  max , Δ 3NB follows an empirically observed geometric progression Δ 3NB () ∼   .Hence the residual uncertainty can be estimated from the value obtained for the largest available  (3)  max = 24 times ∕(1 − ), which for presently investigated systems amounts to a few tens of keV on total binding energies.

Many-body truncation error
The error  CC relates to the truncation of the BCC expansion, i.e. to the omission of high-rank cluster operators.At the BCCSD level, the error is essentially carried by missing triple excitations.In this respect, available VS-IMSRG(2) results [22] obtained in open-shell calcium isotopes with the same Hamiltonian constitute a reference given that static correlations are fully accounted for via the diagonalization within the A similar error model is also built for two-neutron separation energies ( 2 ).Indeed, the difference between BCCSD and VS-IMSRG  2 is shown to be well approximated by taking 0.1% of the BCCSD correlation energy in open-shell Ca isotopes (see Fig. 2).

Particle-number violation error
The error  PNR is associated with the unphysical residual breaking of particle-number symmetry carried by the approximate BCCSD solution.This attractive contribution, presently estimated for each studied nucleus via effective field theory [67], is characterized in Sec. 6.

Calcium isotopic chain
Results of ab initio calculations along the calcium isotopic chain ( = 20) are reported in Fig. 2. For the employed Hamiltonian 60-70% of the total binding energy is captured at the HFB level while the remaining is due to many-body correlations.Due to the softness of the "EM1.8/2.0"Hamiltonian, BMBPT(2) results are in good agreement with BCCSD.Eventually, BCCSD energies typically underestimate known experimental masses ( ≤ 58) by about 5 − 13 MeV (1.7 − 2.9%).
As testified by available VS-IMSRG(2) results [22], missing triple excitations would bring BCC calculations consistently closer to the data.Two-neutron separation energies Calculations are carried out using  max = 12 and  (3) max = 24 with an oscillator frequency of ℏΩ = 16 MeV.In all cases the "EM1.8/2.0"Hamiltonian is employed.
Experimental values are shown as black bars [70].Bayesian uncertainties for  2n values from VS-IMSRG are taken from Ref. [22].are also displayed in Fig. 2. Available experimental values [68,69] are reproduced within the 1-2 MeV many-body uncertainty.The drops at 52,54 Ca indicate shell closures at neutron number  = 32, 34.Furthermore, the location of the neutron dripline ( 2 (, ) < 0) is predicted at  = 60 from both BMBPT(2) and BCCSD in correspondence with the closing of the 1 5∕2 single-particle shell.These finding are consistent with VS-IMSRG(2) results.
Nuclear stability is commonly assessed from the excitation energy of the first 2 + 1 state and its electromagnetic transition to the ground state.It can also be quantified from nuclear masses by computing the so-called two-neutron shell gap

Nickel isotopic chain
Fig. 3 shows that predicted ground-state energies are consistent with experimental values up to the most neutron-rich observed isotope ( 78 Ni).The systematic underbinding is again due to currently discarded triple excitations.Neutron-rich nickel isotopes ( = 50 and beyond) provide a prime example for the limitations of valence-space techniques.In 78 Ni the use of a 0ℏΩ space entails FCI dimensions of dim(  ) ∼ 2 ⋅ 10 10 that cannot be handled by the most advanced shell-model code without severely truncating the configuration basis [71][72][73].Hence the use of no-core approaches such as BCC provide a scalable solution in cases where the VS dimension becomes intractable and/or truncations of the shell model basis induce sizeable uncertainties.
The predicted dripline differs by two neutrons for BMBPT(2) ( 88 Ni) and BCCSD ( 86 Ni): the slightly positive BMBPT(2) value  2 (88, 28) = 35 keV is contrasted by the BCC prediction  2 (88, 28) = −388 keV yielding 86 Ni as the last bound nickel isotope.While convergence with respect to model-space parameters is guaranteed, the relaxation of the many-body truncation and a more adequate treatment of the particle continuum may slightly impact the predicted dripline location.BMBPT(2) and BCCSD two-neutron shell gaps further predict magicity at  = 28, 40, 50.While this is confirmed experimentally for  = 28 and  = 50 [72], it is not the case at  = 40.A more definite statement requires the inclusion of  3 contributions and a systematic variation of the employed chiral Hamiltonian.

Tin isotopic chain
Finally, the most challenging tin ( = 50) isotopic chain is targeted all the way to  = 180.While light tin isotopes were studied via VS-IMSRG calculations based on  (3)  max = 16 [29], heavier isotopes are presently addressed via BMBPT and BCC calculations with  (3)  max = 24.
Based on the error model built in Ca isotopes, Fig. 4 demonstrates that BMBPT(2) and BCCSD binding energies are consistent with experiment up to  = 138 [70], although with a significant uncertainty due to missing correlations.
In Fig. 5,  2 from BCCSD calculations are displayed for three different values of the harmonic oscillator frequency around the empirical minimum of ℏΩ = 12 MeV.While  FBS was globally estimated to be around 2% along the tin isotopic chain, the use of three different ℏΩ values testifies that such a basis truncation error is enhanced in the most neutron-rich isotopes and of the same order as  CC dominating the uncertainty band. 2 While BCCSD  2 are in good agreement with experiment up to 132 Sn, they are significantly too low from 134 Sn till 138 Sn. 3 The difference to the data relates to an overestimation of the  = 82 magic shell gap and is inconsistent with our many-body error estimate.Consequently, it is probably attributable to the Hamiltonian uncertainty. The

Particle-number breaking and restoration
While the HFB reference state breaks particle-number conservation in open-shell systems, Fig. 6 displaying HFB, BMBPT(2) and BCCSD neutron-number variances Δ 2 ≡ ⟨( − ⟨⟩) 2 ⟩ demonstrates that the "EM1.8/2.0"Hamiltonian induces only weak pairing correlations at the mean-field level.Indeed the HFB variance is very close to the minimal boundary corresponding to the zero-pairing limit [76] along Ca, Ni and Sn isotopic chains except, typically, for the most neutron-rich open-shell isotopes.Because |Ψ⟩ is an eigenstate of  , one expects that the symmetry violation is reduced by improving on the BMBPT/BCC truncation.While the neutron-number variance essentially increases in BMBPT(2) calculations, it typically only displays a slight reduction at the BCCSD level compared to HFB.This shows both that the variance reduction is not monotonic in BMBPT [42,77] and that the symmetry restoration cannot effectively be achieved via the summation of low-rank individual excitations.This calls for an explicit particle-number-projection in open-shell nuclei.While this is a standard procedure at the mean-field level [50], extensions to correlation expansions have been explored only recently [28,[78][79][80][81].In absence of an explicit projection, an effective field theory (EFT) is presently used to estimate, at each truncation order, the missing energy contribution  PNR ∼  −1 Δ 2 , where  is a characteristic low-energy scale [67].The bottom panel of Fig. 6 illustrates that this (non-size-extensive) contribution decreases from about 1 MeV in the Ca-Ni region to about 0.5 MeV in Sn isotopes with significant local variation associated with the filling of successive neutron shells.

Conclusions and perspectives
In this Letter, large-scale ab initio non-perturbative Bogoliubov coupled cluster (BCC) calculations were performed for the first time to systematically describe ground-state energies along Ca, Ni and Sn isotopic chains from first principles.These constitute the heaviest openshell nuclei ( ∼ 180) ever computed within an ab initio framework.Such an achievement was made possible by extending to  (1)-breaking expansion methods such as BCC and the optimized storage scheme of chiral three-body forces recently employed in the computation of doubly closed-shell 132 Sn [31] and 208 Pb [32].
Various research directions open up for the near future.First the BCC truncation will be relaxed to incorporate the effect of triple excitations [23,25,82] and reduce the dominating many-body uncertainty.To further reduce uncertainties in heavy neutron-rich isotopes, calculations will be pushed to  max > 12.To cope with the increased dimensionality, this will be combined with basis-optimization techniques accelerating the convergence of many-body observables [25,83,84].Moreover, the formulation of an equation-of-motion [51] BCC extension will provide access to the spectroscopy of even-even and even-odd open-shell nuclei beyond  = 100.Eventually, the implementation of a particle-number restoration [79,80] is envisioned to explicitly account for  PNR as well as to compute rotational excitations and pair-transfer probabilities.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
valence shell.Consequently, a model is presently built for  CC by scaling the difference between BCCSD and VS-IMSRG(2) binding energies in (open-shell) Ca isotopes (see Fig.2) with the BCCSD correlation energy.The former is systematically attractive and indeed shown to amount to about 8% of the latter along open-shell Ca isotopes.Eventually,  CC dominates the many-body error budget (Eq.(4)).With this model at hand,  CC can be propagated to Ni and Sn isotopes based on our sole BCCSD results.Of course, the goal in the future is to avoid the use of this rough error model by explicitly including (perturbative) triple corrections in BCC calculations.
increase of which indicates an extra stability associated in a mean-field picture with a closed-shell nucleus displaying a large Fermi gap.The bottom panel of Fig.2reveals pronounced kinks at  = 20, 28 and less pronounced kinks at  = 32, 34 in agreement with experimental data.Eventually, BMBPT(2) and BCCSD predict a significant shell closure at  = 40 that is however less pronounced in VS-IMSRG(2) calculations.This feature remains to be elucidated experimentally.

Fig. 4 .
Fig. 4. Ground-state energy for 100−180 Sn.Results from HFB, BMBPT(2) and BCCSD are compared to experimental data.Model-space parameters are identical to the ones employed in Figs. 2 and 3.
flat evolution of BCCSD  2 beyond 140 Sn leads, within estimated many-body uncertainties, to a large uncertainty on the location of the predicted neutron drip-line:  ∈ [140 − 162].In spite of its large span, the predicted interval is inconsistent with the energy density functional (EDF) prediction,  ∈ [172 − 176], from Ref. [74] (partially) accounting for systematic and statistical uncertainties of the EDF method.The location of the BCCSD prediction at significantly smaller neutron numbers clearly originates from the too large drop at  = 82.Consequently, the predicted interval should be enlarged to become consistent with the EDF prediction once the Hamiltonian uncertainty is included.Eventually, the uncertainty should be drastically reduced by including triple excitations and going to larger  max on the one hand and by reducing the (presently omitted) Hamiltonian uncertainty on the other hand.

Fig. 5 .
Fig. 5. BCCSD two-neutron separation energies for 100−180 Sn using three different values for the harmonic oscillator frequency around the empirical minimum of ℏΩ = 12 MeV are compared to experimental data.The gray box locates the drip-line prediction from energy density functional calculations (including statistical uncertainty from model parameters) of Ref. [75].

Fig.
Fig.6.Neutron-number variance (top) and leading-order EFT estimate of the energy gain from neutron-number restoration (bottom) for calcium, nickel and tin isotopic chains at HFB, BMBPT(2) and BCCSD truncation level.The low- expansion methods.Alternatively, static correlations can be accounted for through multi-reference techniques https://doi.org/10.1016/j.physletb.2024.138571Received 15 January 2024; Received in revised form 7 March 2024; Accepted 7 March 2024 that are built upon a non-product reference state [43]e particle-number symmetry is spontaneously broken by |Φ⟩, the Hamiltonian must be replaced by the (zero-temperature) grand potential Ω ≡  −    −   .The neutron (proton) chemical potential   (  ) is tailored to constrain |Φ⟩ to carry the physical number of neutrons N 0 (protons Z 0 ) on average.Building on such a reference state, Bogoliubov coupled-cluster theory[43]parametrizes the correlated ground-state wave function as |Ψ⟩ ≡   |Φ⟩ where the cluster operator  ≡  1 +  2 + ... sums components inducing 2 quasi-particle excitations on top of |Φ⟩