Ab initio calculation of the potential bubble nucleus Si
Abstract
 Background

The possibility that an unconventional depletion (referred to as “bubble”) occurs in the center of the charge density distribution of certain nuclei due to a purely quantum mechanical effect has attracted theoretical and experimental attention in recent years. Based on a meanfield rationale, a correlation between the occurrence of such a semibubble and an anomalously weak splitting between low angularmomentum spinorbit partners has been further conjectured. Energy density functional and valencespace shell model calculations have been performed to identify and characterize the best candidates, among which Si appears as a particularly interesting case. While the experimental determination of the charge density distribution of the unstable Si is currently out of reach, (d,p) experiments on this nucleus have been performed recently to test the correlation between the presence of a bubble and an anomalously weak splitting in the spectrum of Si as compared to S.
 Purpose

We study the potential bubble structure of Si on the basis of the stateoftheart ab initio selfconsistent Green’s function manybody method.
 Methods

We perform the first ab initio calculations of Si and S. In addition to binding energies, the first observable of interest are the charge density distribution and the charge root mean square radius for which experimental data exist in S. The second observable of interest is the lowlying spectroscopy of Si and S obtained from (d,p) experiments along with the spectroscopy of Al and P obtained from knockout experiments. The interpretation in terms of the evolution of the underlying shell structure is also provided. The study is repeated using several chiral effective field theory Hamiltonians as a way to test the robustness of the results with respect to input internucleon interactions. The convergence of the results with respect to the truncation of the manybody expansion, i.e. with respect to the manybody correlations included in the calculation, is studied in detail. We eventually compare our predictions to stateoftheart multireference energy density functional and shell model calculations.
 Results

The prediction regarding the (non)existence of the bubble structure in Si varies significantly with the nuclear Hamiltonian used. However, demanding that the experimental charge density distribution and the root mean square radius of S are well reproduced, along with Si and S binding energies, only leaves the NNLO Hamiltonian as a serious candidate to perform this prediction. In this context, a bubble structure, whose fingerprint should be visible in an electron scattering experiment of Si, is predicted. Furthermore, a clear correlation is established between the occurrence of the bubble structure and the weakening of the splitting in the spectrum of Si as compared to S.
 Conclusions

The occurrence of a bubble structure in the charge distribution of Si is convincingly established on the basis of stateoftheart ab initio calculations. This prediction will have to be revisited in the future when improved chiral nuclear Hamiltonians are constructed. On the experimental side, present results act as a strong motivation to measure the charge density distribution of Si in future electron scattering experiments on unstable nuclei. In the meantime, it is of interest to perform oneneutron removal on Si and S in order to further test our theoretical spectral strength distributions over a wide energy range.
I Introduction
The possibility that an unconventional depletion, referred to as “bubble”, occurs in the center of the pointnucleon and/or charge density distributions of a nucleus has attracted theoretical as well as experimental attention in recent years. Earlier, socalled (semi)bubble structures had been invoked mainly in connection with hypothetical (super) hyperheavy nuclei characterized by a very large charge () Dechargé et al. (2003). Indeed, singlereference (SR) energy density functional (EDF) calculations predicted that the groundstate configuration of these nuclei may display a depletion in the center of their density distribution Dechargé et al. (2003); Bender et al. (1999) as a result of a collective quantum mechanical effect sustained by the compromise between the large repulsive Coulomb interaction and the strong force that binds them Bender and Heenen (2013).
The case of present interest relates to nuclei characterized by more conventional masses and possibly displaying a semi bubble^{1}^{1}1As for the terminology, we use indistinguishably ”bubble”, ”semi bubble” or ”central depletion” in the following although strictly speaking ”bubble” should be kept for speculative hyperheavy nuclei possibly displaying a null density in their center, which is not the case for the nucleus of present interest. in their center. The rationale here solely relates to the quantum mechanical effect that finds its source in the sequence of occupied and unoccupied singleparticle states near the Fermi energy in an independentparticle or a meanfield picture. While () orbitals display a radial distribution that is peaked at the center of the nucleus, orbitals with nonzero angular momenta () are suppressed in the nuclear interior such that they do not contribute to the central density. As a result, any vacancy of orbitals embedded among larger orbitals near the Fermi level is expected to produce a depletion of the central density. These hypothetical nuclei are of interest as they must be modelled via meanfield potentials that differ from those associated with Fermitype density distributions that fit the vast majority of nuclei. In turn, a nonzero density derivative in the nuclear interior has been conjectured to cause a sharp increase of ”nonnatural” sign of the effective onebody spinorbit potential, eventually inducing a reduction of the splitting between spinorbit partners characterized by low angular momenta ToddRutel et al. (2004); Burgunder (2012).
Going beyond this meanfield scenario, a small energy difference between the unoccupied shell and the last occupied/next unoccupied shells can favor collective correlations and thus lower or even wash out the depletion at the center of the potential bubble nucleus. Therefore the search for the best bubble candidates must be oriented towards nuclei that can be reasonably modelled by an orbital well separated from nearby singleparticle states such that correlations are weak. In turn, this feature underlines the necessity to employ theoretical methods that explicitly incorporate longrange correlations modifying the density on a length scale of about fm, which is the typical expected spatial extent of the depletion at the center of bubble nuclei as discussed below.
In recent years, SR Richter and Brown (2003); ToddRutel et al. (2004); Khan et al. (2008); Grasso et al. (2009) and multireference (MR) Yao et al. (2012, 2013); Wu et al. (2014) EDF calculations along with shell model calculations Grasso et al. (2009) have been performed for O, Si, Ar, Hg and Hg. Indeed, these nuclei appeared as favorable candidates based on the naive filling of singleparticle shells their number of protons and/or neutrons correspond to. Among those, Si (, ) stands out as the most viable case as its depletion factor defined as
(1) 
is predicted to be the highest among all candidates in SREDF calculations. In Eq. 1, and denote central and maximum (pointnucleon or charge) density values, respectively. For , the naive filling of proton shells leaves the singleparticle state as the first unoccupied level above the Fermi energy. Furthermore, the magic character of Si translates into a first excitation energy ( MeV) and a transition probability Ibbotson et al. (1998) that are similar to the doublymagic Ca nucleus^{2}^{2}2See Ref. Sorlin and Porquet (2008) and references therein for the systematic of and in the isotonic chain.. The low electric monopole transition strength Rotaru et al. (2012) completes the picture of a doublymagic system Baumann et al. (1989). These features leaves the hope that the naive rationale based on an independentparticle model only needs to be slightly perturbed by the inclusion of longrange correlations^{3}^{3}3The perfect counter example is Si (, ) for which the even more promising picture provided by the naive filling of spherical shells is eventually fully invalidated, resulting in a charge density distribution that does not display any depletion in its center Miessen (1982)..
A bubble structure mainly driven by protons, as in Si, can be probed directly by measuring the charge density distribution via electron scattering. However, it is presently not possible to perform electron scattering on unstable nuclei as light as Si with sufficient luminosity. Such an experiment may become feasible in the next decade at ELISeFAIR Simon (2007) or after an upgrade of the SCRIT facility at RIKEN Suda et al. (2009).
Because the presence of the central depletion is believed to correlate with specific quantum mechanical properties and to feedback on other observables, one may think of alternative ways to probe its presence indirectly, e.g. via direct reactions. In the present case of interest, we specifically wish to test the correlation between the presence of the bubble and the evolution of the spinorbit splitting when going from S to Si. The establishment of this correlation is performed in the eye of the capacity of ab initio calculations to reproduce the lowlying spectroscopy of nuclei obtained via the addition of a neutron Thorn et al. (1984); Eckle et al. (1989); Burgunder et al. (2014) or the removal of a proton Khan et al. (1985); Mutschler et al. (2016a, b) on S and Si.
While potential bubble nuclei such as Si have already been investigated quite thoroughly within the frame of EDF and shellmodel manybody methods, the goal of the present work is to provide the first study based on ab initio manybody calculations. As mentioned above, our aim is to perform a coherent analysis of both density distributions and oneneutron addition and removal spectral strength distributions. Ideally, one would like to further correlate these observables with spectroscopic information in Si itself as was done in Refs. Yao et al. (2012, 2013). However, the manybody scheme employed does not allow to do it yet. This will hopefully become possible in a not too distant future. Also, one of the objectives of the present study is to characterize the sensitivity of the results to the input Hamiltonian and to outline the role of threenucleon forces.
The paper is organized as follows. Section II focuses on the computational scheme employed, paying particular attention to the convergence of the observables of interest with respect to the basis used to represent the Schrödinger equation and to the manybody truncation implemented to solve it. Section III analyses in detail the characteristic of pointproton and charge density distributions of S and Si. The impact of manybody correlations and the sensitivity of the results to the utilized Hamiltonian are discussed. Results from our ab initio calculations are further compared to those obtained from stateoftheart MREDF and SM calculations. Section IV concentrates first on the reproduction of the spectroscopy of neighboring S, P, Si and Al. In particular we correlate the evolution of the spinorbit splitting when going from S to Si and the presence/absence of a bubble structure in the ground state of Si/S. The interpretation in terms of the evolution of the underlying onenucleon shell structure is also provided.
Ii Computational scheme
ii.1 Manybody methods and Hamiltonians
The present analysis is based on selfconsistent Green function (SCGF) calculations Dickhoff and Barbieri (2004); Barbieri and HjorthJensen (2009); Somà et al. (2011) of Si and S. This is done employing the following scheme

By default, a combination of twonucleon (2N) and threenucleon (3N) interactions obtained within the frame of chiral effective field theory (EFT) at nexttonextto leading order (NLO) and denoted as NNLO Ekström et al. (2015) is used. For comparison, we occasionally employ a different set of NLO interactions Ekström et al. (2013), denoted as NNLO, as well as the combination of NLO 2N and NLO 3N chiral forces (NN+3N400) of Refs. Entem and Machleidt (2003); Navrátil (2007); Roth et al. (2012) evolved down to a lowmomentum scale by means of a similarity renormalization group (SRG) transformation Bogner et al. (2010). The comparison of the results based on the three sets of EFT interactions is of interest given that NNLO was specifically designed Ekström et al. (2015) to address the impossibility of any existing set to convincingly describe binding energies and nuclear radii (saturation) of midmass nuclei Somà et al. (2014a); Hergert et al. (2014); Lapoux et al. (2016) (infinite nuclear matter Carbone et al. (2014); Hagen et al. (2014)) at the same time.

Calculations expand one, two and threebody operators on a spherical harmonic oscillator basis containing up to 14 harmonic oscillator (HO) shells [ max = 13]. All matrix elements of one and twobody operators are used whereas those of 3N interactions are limited to configurations characterized by ++=16.

As nuclei under study display a closed (sub)shell character in a naive independentparticle approximation, both Dyson Dickhoff and Barbieri (2004); Barbieri (2009); Cipollone et al. (2013) and Gorkov Somà et al. (2011, 2014b, 2014a) SCGF calculations can be safely performed. The latter framework is employed to test whether the explicit inclusion of pairing correlations via the breaking of U(1) global gauge symmetry associated with particlenumber conservation modifies, e.g. improves, the theoretical predictions or not. The manybody truncation scheme is based on the socalled norder adiabatic diagrammatic construction (ADC(n)) Schirmer et al. (1983). While Gorkov SCGF calculations can currently be performed at ADC(1), i.e. HartreeFock(Bogoliubov) (HF(B)), and secondorder (i.e. ADC(2)) levels, Dyson SCGF further accesses ADC(3) calculations. While ADC(2) calculations already resum the bulk of dynamical correlations, ADC(3) provides wellconverged bulk properties of the Abody system of interest along with a refined description of spectroscopic properties of A1 nuclei.
ii.2 Convergence of groundstate energies
The groundstate energy of Si computed at the ADC(2) level is displayed in Fig. 1 as a function of the oscillator frequency and for increasing values of . As the number of major shells increases, the energies become more independent of while their minima shift progressively towards lower values to eventually reach MeV for . The behavior is similar for S.
In principle, the choice of a specific harmonic oscillator frequency is arbitrary in the limit of very large as manybody quantities must become independent of it. At workable values of , observables are only approximately independent of and the optimal value of the latter may differ from one observable to the other and from one eigenstate to the other. For example, lower values of than the one minimizing the groundstate energy might better approach the infinitebasis value for longrange, e.g. meansquare radii, operators Shin et al. (2016). As discussed next, however, in the present case there is little impact of the specific value of on density distributions, which constitute the focus of the present paper. Consequently, and given the lack of a well defined extrapolation procedure for density distributions, the value MeV corresponding to the minimum of the energy for is considered in the following sections.
ADC(1)  ADC(2)  ADC(3)  Experiment  

Si  84.481  274.626  282.938  283.427 
S  90.007  296.060  305.767  308.714 
Groundstate energies computed at various orders in the manybody truncation scheme are compared to experimental data in Tab. 1. At the ADC(2) level, theoretical results are within of experimental data, which is consistent with missing ADC(3) correlations and the intrinsic uncertainty of the input Hamiltonian Lapoux et al. (2016); Ekström et al. (2015). Going to ADC(3) indeed brings about 810 MeV additional binding, which represents about of the correlation energy generated at the ADC(2) level. Extrapolating the pattern of reduction in the correlation energy added at each ADC(n) order, the ADC(3) results can be safely believed to be about 12 MeV (i.e. less than ) away from the fully converged values. With the presently used NNLO Hamiltonian, this happens to be of the order of the difference to experimental data.
ii.3 Convergence of groundstate radii
Before addressing pointnucleon and charge density distributions, let us focus on the integrated information constituted by pointnucleon and charge rootmeansquare (rms) radii. In Fig. 2, ADC(2) calculations of the charge rms radius^{4}^{4}4In the present work charge radii are computed from pointproton radii by accounting for the finite charge radii of both protons and neutrons in addition to the DarwinFoldy correction, see Ref. Cipollone et al. (2015) for details. of Si are displayed for different values of and . As increases, the dependence on becomes weaker, totalling to about for for MeV. Table 2 reports charge rms radii of Si and S computed within different manybody truncation schemes. The convergence pattern is similar for the two nuclei, with tiny differences between ADC(2) and ADC(3) results. This indicates that rms radii are essentially converged already at the ADC(2) level.
It is currently a challenge for ab initio calculations to describe both the binding energy and the size of mediummass nuclei at the same time Lapoux et al. (2016). This situation lead recently to the construction of the (unconventional) NNLO EFT Hamiltonian Ekström et al. (2015) that is presently used and that indeed improves the situation significantly Garcia Ruiz et al. (2016); Lapoux et al. (2016). The computed value fm in S is very close to the experimental measurement. Comparatively, the rms charge radius computed from the NN+3N400 Hamiltonian processed through a SRG transformation is significantly too small, e.g. it is predicted to be fm at the ADC(2) level for fm.
Experimental charge radii are unavailable for the unstable Si nucleus. While charge radii for stable isotopes can be measured by means of electron scattering, collinear laser spectroscopy experiments Campbell et al. (2016) currently constitute the most appropriate way to access charge radii of unstable nuclei with lifetimes as low as a few milliseconds. However, Si elements are highly reactive and require a high evaporation temperature, thus are extremely difficult to produce and extract via ISOL techniques. Inflight facilities, e.g. NSCL at Michigan State University, are able to provide highintensity beams of Si isotopes. Future developments of highresolution laser spectroscopy experiments should enable a measure of the rms charge radius of Si Garcia Ruiz (2016).
ADC(1)  ADC(2)  ADC(3)  Experiment  

Si  3.270  3.189  3.187   
S  3.395  3.291  3.285  3.2985 0.0024 
Si  3.085  3.258  3.188  3.187 

S  3.184  3.285  3.240  3.285 
For completeness, pointproton, pointneutron, matter and charge radii computed at the ADC(3) level are reported in Tab. 3. Recent works have demonstrated that matter radii can be reliably extracted from elastic proton scattering data Lapoux and Alamanos (2015); Lapoux et al. (2016). Given that such data are available for several sulfur isotopes including S Maréchal et al. (1999), it would be interesting to compare a similar evaluation of to our present results.
ii.4 Convergence of pointnucleon densities
Expanding the manybody Schroedinger equation on a HO basis, the spurious centreofmass (COM) contribution to the density distribution can be removed exactly as long as the COM part of the manybody state factorizes Navrátil (2004), which is authorized by truncating the basis in term of a fixed number of Abody HO excitations. It is not presently ensured given that the basis truncation is performed directly at the level of onebody, twobody and theebody Hilbert spaces. It happens that such a truncation procedure leads to an effective factorization of the COM part of the wave function as demonstrated in coupledcluster calculations Hagen et al. (2009). Although it remains to be explicitly demonstrated for SCGF calculations, the proximity of both methods and of their results gives confidence that a potential contamination of the density distribution is, at most, small in this nuclei. More on this point in Sec. III.4.
The pointproton density distribution of Si is displayed in Fig. 3 for different values of at MeV and for different harmonic oscillator frequencies at . The overall profile shows a very weak dependence on the model space parameters. In particular, the dependence of the central density on is minor such that considerations about the potential bubble structure are little affected. As expected from the use of a harmonic oscillator basis, the asymptotic of the density distribution is altered by the change of , which however does not impact the analysis presented below.
Given that we are primarily interested in features associated with potential bubble structures, we compute the pointproton parameter (Eq. 1) at in order to estimate the uncertainty associated with the choice of . For MeV we find , where the error is the standard deviation of the dataset. This shows that the uncertainty associated with the choice of on the depletion factor is negligible compared to other sources of error (see below).
Iii Bubble structure
iii.1 Pointnucleon density distributions
The onebody density matrix associated with the ground state of Si or S is computed in the HO basis of the onebody Hilbert space through
(2) 
where denotes the (normal part of the) onebody Green’s function in the energy representation Somà et al. (2011) and where the integral is performed along a closed contour located in the upper half plane of the complex plane. The pointproton density distribution reads as
(3) 
where denotes HO singeparticle wavefunctions and where the sum is obviously restricted to proton singleparticle states. Similar expressions hold for pointneutron and matter density distributions.
Theoretical pointproton and point neutron density distributions of Si and S computed at the ADC(3) level are displayed in Fig. 4. Most clearly, the pointproton density distribution of Si displays a pronounced depletion at the center while the one of S presents a maximum. This results in a noticeable (null) pointproton depletion factor () in Si (S). On the other hand, pointneutron density distributions are very similar, i.e. the creation of the proton bubble in Si associated with the removal of two protons from S affects the spatial distribution of neutrons only weakly by pulling them slightly away from the very center to the maximum of the proton density around fm.
iii.2 Theoretical analysis
Pointnucleon density distributions can be analyzed, internally to the theoretical scheme^{5}^{5}5See Sec. IV.4 for a brief discussion on the nonobservable character of quantities that are internal to the theory, i.e. that have no counterpart in the empirical world Duguet et al. (2015), such as the presently introduced singleparticle state occupations Furnstahl and Hammer (2002)., by expressing them in the socalled natural basis , i.e. in the orthonormal basis of that diagonalizes the onebody density matrix . Natural orbitals are thus obtained in the HO basis by solving the eigenvalue equation
(4) 
In the natural basis, the pointproton density distribution reduces to a sum of positive singleparticle contributions. Taking into account the spherical symmetry characterizing the ground state of the nuclei of interest, the pointproton density distribution eventually reads as
(5) 
where denotes the radial part of the natural singleparticle wavefunctions and where the partialwave contributions to the density have also been introduced.
Figure 5 displays the partialwave decomposition of pointproton density distributions of Si and S. In Si, the very interior of the density is entirely built from the partialwave and is depleted compared to its maximum at about fm that is dictated by the and partialwaves. The wave contributes at larger radii and dominates at the surface. One also observes small contributions from the and waves.
Ab initio calculations describe the complete dynamics of the interacting nucleons in large model spaces. Consequently, each partial wave builds from several orbitals corresponding to different principal quantum numbers . Figure 6 decomposes the occupation of each partial wave into individual proton natural orbital occupation. As expected, the partialwave is dominated by the protons occupying the states. Still, the () protons occupying the () states do contribute for () of the density at . This surprisingly large contribution originates from the fact that the () wavefunction is () times larger than the wavefunction at (see Fig. 7).
As visible in Fig. 5, the change in pointproton density from Si to S entirely originates from the partial wave. The proton density almost doubles at when adding two protons to Si, leading to the disappearance of the bubble structure in S. As expected from a naive independent particle picture, the rise of the central density is driven by the states whose occupation increases from in Si to protons in S. Interestingly though, this increase in occupation is somewhat compensated by a lowering of the associated wave function at , as is visible from Fig. 7. Consequently, the rise of the central density of S due to the increase of is only 55 of what it would have been with a frozen natural orbital wavefunction. Correspondingly, the () states decrease (increase) the central density by an amount that corresponds to 20 (8) of the rise generated by the states. These subleading contributions mainly originate from the fact that the two added protons lead to a lowering (rising) of the () natural wavefunctions at . Overall, even though the disappearance of the proton bubble when going from Si to S does mainly reflect the increased occupation of the states, one observes that in quantitative terms the net effect results from the combination of several intricate features in our ab initio calculations.
As testified by Fig. 8, the fact that pointneutron density distributions of Si and S are (nearly) identical reflects the (essentially) equal occupations of neutron natural orbitals and their (essentially) unchanged wave functions.
iii.3 Impact of correlations
In panel (a) of Fig. 9, pointproton density distributions of Si and S are compared at various levels of manybody truncations, i.e. as obtained from Dyson ADC(1), ADC(2) and ADC(3) calculations. Moving from ADC(1) (i.e. HF) to ADC(2), the amplitude of the central depletion diminishes in Si. This erosion of the bubble structure is the consequence of explicit dynamical correlations added at the ADC(2) level, knowing that correlations added at the ADC(3) level do not further change the picture. Eventually, this erosion results in a decrease of the pointproton F factor from to when going from HF to ADC(3) calculations (see Tab. 4). The impact of correlations on the pointproton density distribution of S is less pronounced.
Si  ADC(1)  ADC(2)  ADC(3) 

0.49  0.34  0.34 
The impact of manybody correlations can be analyzed on the basis of the natural orbital decomposition of the density. Given that the onebody density matrix reflects the correlations included in the calculation of , it is clear that not only the natural occupations but also the natural orbital wavefunctions change with the manybody truncation scheme employed.
First, dynamical correlations partially promote protons from (, , , ) states into (, , , , ) as is visible from Fig. 11. Second, longrange correlations lead to a contraction of the proton natural wavefunctions, the effect being the most drastic for the orbitals that are originally unoccupied at the ADC(1) level, e.g. for the and wavefunctions on Fig. 10. This is particularly true for the wavefunction that becomes strongly peaked at . In spite of partially promoting protons into orbitals that are located more towards the outside, longrange correlations induce a reduction of the charge rms radius of the nucleus associated with the contraction of the natural orbital wave functions. With the presently used NNLO Hamiltonian, this further improves the agreement of the predicted charge rms radius of S with experimental data as visible in Tab. 2.
The net result of manybody correlations on the pointproton density of Si is visible in Fig. 12, where the variation due to each partial wave is displayed. One observes that the nonzero occupation of the states and the contraction of their wave function increases the central density by about 18 of its maximum value, leading to the erosion of the bubble structure mentioned above. The fact that the bubble structure can be significantly suppressed by the inclusion of longrange correlations was already pointed out on the basis of MREDF Yao et al. (2012, 2013) and shellmodel (SM) Grasso et al. (2009) calculations.
On panel (b) of Fig. 9, pointproton density distributions of Si and S obtained from both Dyson and Gorkov SCGF calculations at the ADC(2) level are compared. Pairing correlations are very weak in these two closed subshell nuclei such that their explicit account does not impact the results in any significant way. This feature, presently obtained on the basis of realistic 2N and 3N internucleon interaction, mirrors the situation at play in SREDF calculations Grasso et al. (2009).
iii.4 Charge density distribution
Generically speaking, the electromagnetic charge density (and current) operator is expressed as an expansion in manybody operators acting on nucleonic degrees of freedom. This operator not only accounts for the point distribution of protons but also for their own charge distribution, along with the one of neutrons, and for charge (and current) distributions associated with the light charged mesons they exchange. To first approximation, the nuclear charge density can be obtained through the folding of the nuclear pointproton density distribution with the charge density distribution of the proton. In doing so, one omits neutrons’ contribution^{6}^{6}6We remind however that neutron’s charge density contribution to charge radii is presently taken into account Cipollone et al. (2015). as well as relativistic spinorbit corrections, both typically relatively small Chandra and Sauer (1976). We thus compute the charge density distribution as
(6) 
where the sets come from having parameterized the charge density distribution of the proton as a linear superposition of three Gaussians and have been adjusted to reproduce the proton charge form factor from electron scattering data Chandra and Sauer (1976). The proton rms radius that results from this parameterization is , consistent with the value used to compute the rms charge radius Mohr et al. (2012). Let us note that eventual smaller values of Pohl et al. (2010) would lead to an increase of the depletion factor (see also Tab. 8 and corresponding discussion).
Furthermore, one needs to correct for spurious center of mass and include DarwinFoldy relativistic correction. Assuming^{7}^{7}7While this has been proven for coupledcluster Hagen et al. (2009) and inmedium similarity renormalization group Hergert et al. (2016) calculations, a similar study remains to be done for SCGF calculations. Given the proximity of these manybody methods, one is however confident that the center of mass factorization does indeed occur in the same way in SCGF calculations as is assumed here. that the center of mass wavefunction factorizes in the groundstate of a harmonic oscillator Hamiltonian characterized by the frequency , the inclusion of spurious centerofmass and DarwinFoldy relativistic corrections can be performed at the price of proceeding to the replacement Negele (1970); Chandra and Sauer (1976)
(7) 
in Eq. 6, where is the nucleon mass, hence ^{8}^{8}8We use here, as everywhere throughout the paper, . fm, and . Employing Bethe’s formula Negele (1970), the latter term can be approximated with . We note that, for O, such an approximation is consistent with the value of found in Ref. Hagen et al. (2009) and is thus safe to use in present calculations of Si and S.
Theoretical charge density distributions of Si and S computed at the ADC(3) level are compared to their pointproton counterpart in Fig. 13 and to the experimental charge density of S Rychel et al. (1983). The excellent agreement between theoretical and experimental charge density distributions of S gives confidence in the SCGF prediction obtained with NNLO for Si. While the folding operated in Eq. 6 weakly reduces the peak at the center of the density distribution of S, it significantly smears out the depletion in the pointproton density distribution of Si. This effect could be expected given that the folding takes place over a distance set by the proton charge radius that is consistent with the size of the proton bubble. The fact that the bubble structure could be strongly suppressed in the observable charge density of Si was already pointed out on the basis of SR and MREDF calculations Yao et al. (2012, 2013). This reflects in the value of the F factor of Si that goes down from to when going from the pointproton to the charge density distribution (see Tab. 8).
iii.5 Form factor
Accessing the charge density distribution of Si would require to scatter electrons on radioactive ions. This would lead to measuring the electromagnetic charge form factor, which relates to the nuclear charge density distribution through
(8) 
where is the transferred momentum, itself related to the incident momentum and the scattering angle via .
In Ref. Khan et al. (2008), HartreeFock densities based on two Skyrme EDF parameterizations were used to demonstrate that the diffraction pattern of a semibubble nucleus differs significantly from the one of the same nucleus without a bubble. Similarly, Fig. 14 displays the angular dependence of the form factor obtained at the ADC(2) and ADC(3) levels for 300 MeV electron scattering on Si and S. From 60 to 95 degrees, the angular distribution is located at higher magnitude in Si than in S. This is accompanied by a shift of about 20 degrees between both second minima such that the two angular distributions are out of phase at about 110 degrees. Furthermore, the comparison between ADC(2) and ADC(3) results demonstrates that this pattern is solid against the further addition of manybody correlations whenever the bulk of them have been included. All in all, the depletion in the charge density of Si is predicted to be large enough for the fingerprint of the bubble structure to be potentially visible in a future electron scattering experiment.
iii.6 Choice of the Hamiltonian
We wish to probe the sensitivity of the results to the choice of the input nuclear Hamiltonian. Ideally, one would like to propagate systematic uncertainties associated with the EFTbased Hamiltonian to manybody observables of interest^{9}^{9}9This should be performed at various chiral orders, i.e. going from LO to NLO, to NNLO etc, in order to check the consistency of the extrapolated uncertainties.. Although efforts in that direction are currently being made Carlsson et al. (2016), it is not possible to implement this program convincingly yet. Consequently, we presently limit ourselves to repeating the calculation for the four representative EFTbased Hamiltonians introduced in Sec. II.1.
S  NN+3N400(1.88)  NNLO  NNLO 

2.864  3.033  3.291 
The charge density distribution of Si (S) computed at the ADC(2) level from the four different Hamiltonians is displayed in the left (right) panel of Fig. 15. Starting with S for which experimental data exist, one can appreciate the significant spread of the results. While the charge density obtained with NNLO is in close agreement with data, those obtained from NNLO, NN+3N400(1.88) and NN+3N400(2.24) are much too large in the bulk and do not extend far enough. The shortcoming of the last three Hamiltonians relates to their incapacity to account for both nuclear binding and nuclear size at the same time Somà et al. (2014a); Hergert et al. (2014); Lapoux et al. (2016). As visible in Tab. 5, the pattern in the charge density does correlate with the charge rms radius. As a matter of fact, fixing the latter by using charge radii of C and O in the fit of NNLO, the overall charge density is very well reproduced. Empirically speaking, and omitting the questions raised by this way of fixing the problem, this result strongly favors the prediction made with NNLO in the present study.
The variability of the charge density of Si follows the same pattern as in S and thus reflects the capacity of the Hamiltonian to account for nuclear sizes. The choice of the Hamiltonian slightly impacts the radial extension of the bubble but not its depth on an absolute scale that appears to be a robust prediction. However, as the density in the bulk overall shrinks as one goes from NN+3N400, to NNLO and to NNLO, the F factor grows accordingly as reported in Tab. 6. Given the empirical superiority of NNLO, the corresponding value ( at both ADC(2) and ADC(3) levels) provides the most probable ab initio characterization of the charge bubble in Si as of today.
Last but not least, one notices that varying the momentum scale of the SRG transformation applied to the NN+3N400 Hamiltonian from fm to fm has negligible impact on the charge density of Si and S. However, it is to be noticed that comparing results obtained from the NN+3N400 Hamiltonian to which a SRG transformation was applied to those obtained from the unevolved NNLO Hamiltonian raises the question of the SRG transformation of the charge density operator in the former case. Although such an evolution is unlikely to impact the bubble structure significantly, it is a caveat to be mentioned here and investigated in the future.
Si  NN+3N400(1.88)  NNLO  NNLO 

0.08  0.11  0.15 
iii.7 Impact of 3N interactions
To single out the effect of the 3N interaction, charge density distributions of Si and S computed at the ADC(2) level with and without it are displayed in Fig. 16 for both NNLO and NN+3N400(1.88) Hamiltonians.
For NNLO, the inclusion of the 3N interaction has a severe impact on the charge density distribution and on the rms charge radius of S, which is eventually in good agreement with data as already discussed above. Generically, the charge density is strongly reduced in the bulk and pushed towards the outside. The maximum of the density in Si moves down and away from the center, i.e. from about fm to about fm, enlarging radially the bubble structure. At the same time, the 3N force increases significantly the depletion in the center. Given that the charge density overall shrinks in the bulk, this leads to a pronounced enlargement of the F factor from with 2N interaction only to when the 3N interaction is incorporated.
At first sight, the effect of the 3N interaction, which first enters at NNLO in Weinberg’s power counting Bedaque and van Kolck (2002); Nogga et al. (2005); Epelbaum et al. (2009), seems anomalously large compared to the 2N only result. To investigate this point thoroughly, one would need to perform the calculation at various orders in the EFT power counting, not just remove the NNLO 3N contribution from the NNLO Hamiltonian as done here. If a strong impact of the NNLO contributions to the Hamiltonian were to be confirmed, it would additionally question the present use of a leadingorder approximation to the chargedensity operator^{10}^{10}10Independently of that, it is anyway necessary in principle to use consistent Hamiltonian and charge operators, e.g. see Ref. Krebs et al. (2016).. While being beyond the scope of the present paper, these issues need to be investigated in the future.
Si  

NNLO  2N  0.01 
NNLO  2N+3N  0.15 
NN+3N400(1.88)  2N  0.08 
NN+3N400(1.88)  2N+3N  0.08 
For NN+3N400(1.88), the impact of the 3N interaction is much weaker. In S, the 3N interaction only slightly pushes the charge density away from the center, which is not sufficient to bring the charge radius in agreement with data. The charge density of Si is essentially untouched and so is the bubble structure. Independently of the inclusion of the 3N force, the F factor remains moderate and equal to .
The present analysis demonstrates that, as expected on general grounds, the effect of the 3N interaction depends of the Hamiltonian, i.e. it is a function of the partner 2N interaction and of the overall quality of the Hamiltonian it is part of. When considering the overall (2N+3N) Hamiltonian, only the latter remains and the variability of the results is more transparent and systematic, as discussed in Sec. III.6.
iii.8 Comparison to SM and MREDF calculations
We wish to compare results of ab initio SCGF calculations based on the NNLO Hamiltonian to those obtained from alternative theoretical methods. As EDF calculations of Refs. Yao et al. (2012, 2013) demonstrated, correlations beyond any HartreeFock or static SREDF method significantly affect potential bubble structures. Consequently, we focus on MREDF calculations that do incorporate longrange correlations explicitly. Table. 8 compares proton and charge depletion factors obtained in Si from present ab initio SCGF calculations to those obtained from stateoftheart MREDF Yao et al. (2012, 2013) and SM calculations. It is to be noticed that the experimental charge density of S is equally well reproduced by ab initio SCGF and MREDF calculations. The reproduction is fair in the SM calculation, knowing that the charge density is computed by weighting onebody wavefunctions derived from a well tailored WoodsSaxon potential with the singleparticle occupations obtained from the SM calculation.
Multireference EDF calculations of Refs. Yao et al. (2012) and Yao et al. (2013) are based on nonrelativistic and relativistic energy functionals, respectively. Both sets of calculations include longrange correlations associated with the restoration of and symmetries, i.e. with the restoration of good nucleon numbers and angular momentum, as well as with the (large amplitude) fluctuation of the intrinsic axial quadrupole deformation. The calculations are benchmarked against the known spectroscopy of Si (i.e. , , and ). While the agreement is satisfactory in both cases, the calculations based on a the relativistic pointcoupling PCPK1 parameterization Zhao et al. (2010), augmented with a contact interaction to treat pairing correlations within the BCS approximation, are particularly performing and display an appropriate amount of intrinsic shape mixing in the lowlying states of interest. Shellmodel calculations of Ref. Grasso et al. (2009) employ the USD interaction Brown and Wildenthal (1988) in the sd valence space. As such, longrange correlations are included via the full diagonalization of the effective interaction in the valence space. A good reproduction of and states in Si is however likely to necessitate 2p2h intruders outside the sd shell such that the benchmarking against the known spectroscopy of Si is out of the reach of sd SM calculations Brown (2016). Similarly, the excitation energy of and states in Si is not yet available in ab initio SCGF calculations for benchmarking. However, SCGF calculations can be tested against spectra of neighboring nuclei accessed via, e.g., onenucleon addition and removal experiments. This is postponed to the next section.
Si  SCGF  SCGF*  MREDF Yao et al. (2012)  MREDF Yao et al. (2013)  SM Grasso et al. (2009) 

0.34  0.34  0.21  0.22  0.41  
0.15  0.19*  0.09  0.11  0.28 
As seen in Tab. 8, charge depletion factors of Si obtained from both MREDF calculations are essentially identical and predicted to be around . The depletion factor is predicted to be larger in ab initio calculations (*) and significantly larger in SM calculations (), with the caveat that the charge density is obtained through a rather had hoc procedure in the latter case. Interestingly, the difference between SCGF and MREDF calculations essentially originates from the underlying spherical meanfields. While present ADC(1) calculation predicts *, the charge depletion factor is in the SREDF calculation based on a spherical HartreeFock reference state Yao et al. (2012). Adding longrange correlations, the suppression of the bubble structure is identical in both sets of calculations in spite of the fact that the manybody schemes employed to do so are very different. This may reflect the need for EDF parameterizations to be fitted at the MR level, i.e. once longrange correlations are included. All in all, predictions from present ab initio calculations and from SM calculations leave more hope to observe a bubble structure in the charge density of Si than presentday MREDF calculations.
Iv Spectroscopy
iv.1 Onenucleon addition and removal spectra
The spectral strength distribution displays onenucleon separation energies
(9) 
against spectroscopic factors
(10) 
for all final states of the systems reached by adding/removing one nucleon to/from the Abody groundstate of interest. Spectroscopic factors are computed from onenucleon addition and removal spectroscopic probability matrices defined through
(11a)  
(11b) 
Selfconsistent Green’s function calculations of Si and S ground states automatically access the information on neighboring systems associated with Eqs. 911. For reference, oneneutron and oneproton addition and removal spectral strength distributions associated with the ground state of Si (S) and calculated at the ADC(3) level on the basis of the NNLO Hamiltonian are displayed over a wide energy range in Fig. 17 (Fig. 18). Bars above (below) the dashed line denote states in the nucleus with one nucleon more (less). The length of the bars characterizes the fragmentation of the strength in the present manybody calculation.
iv.2 Comparison to experimental data
To better typify present theoretical predictions and compare them to available experimental data, oneneutron additional energies to the lowestlying states of Si and S () with dominant strength are shown in Fig. 19 against experimental data obtained from (d,p) reactions on Si Burgunder et al. (2014) and S Thorn et al. (1984); Eckle et al. (1989). One notices that the theoretical ordering of the , and states is correct in Si and S and the distance between the states is in fair agreement with data. Even though the theoretical spectra are slightly too diluted, they are incomparably better than with the NN+3N400 Hamiltonian for which they come out much more spread out (see Ref. Papuga et al. (2014) for a typical example in K isotopes). Finally, the three lowestlying separation energies are also well reproduced on an absolute scale, which relates to the fact that the error on the separation energy between the ground states of Si (S) and Si (S) is only keV ( keV). While the 3N interaction is key to obtain the correct density of states, manybody correlations are essential to position the spectrum on an absolute scale, e.g. going from ADC(2) to ADC(3) lowers and () by slightly more (less) than MeV in the direction of experimental data in Si. In view of this, we can speculate that missing ADC(4) correlations might shift onenucleon separation energies of dominant peaks on a level of keV. Overall, the largest disagreement with data relates to the state in Si that is wrongly predicted unbound by keV, i.e. keV above the experimental state. In addition to the possible improvement brought about by a better EFT Hamiltonian and by the inclusion of ADC(4) correlations, this low angularmomentum state near the continuum threshold is probably significantly impacted by the use of a truncated HO basis to expand the manybody Schrödinger equation. The residual sensitivity of this state to the value of () testifies that it is less well converged than the other lowlying states and should probably lie few hundred keV lower.
Similarly, one can analyse oneproton removal energies to the lowestlying states of Al () and P () as obtained from knockout experiments on Si Mutschler et al. (2016b) and S Khan et al. (1985); Mutschler et al. (2016a). The comparison is shown in Fig. 20. In P, the sequence of , and is qualitatively consistent with data, which is a key feature with regards to the occurrence of the proton bubble. Quantitatively though, the separation energy to the ground state is MeV too small whereas excitation energies of the and states are MeV and MeV too large, respectively. The theoretical spectrum displays additional lowlying and states with small spectroscopic strength. However, these small fragments are unlikely to be fully converged with respect to the manybody truncations even at the ADC(3) level.
In Al, no state except for the groundstate has firm spinparity assignment. The groundstate spinparity is reproduced in the calculation although the corresponding separation energy is MeV too small. The experimental spectrum appears much denser, i.e. more fragmented, between and MeV excitation energy than in P. It is not what is obtained theoretically where only two states arise in this energy window. While at least two tentative states are seen experimentally below MeV with small cross sections, the first calculated state is located at about MeV. This situation most probably reflects that Al is on the edge of the socalled island of inversion. Capturing the strong associated quadrupolequadrupole correlations is probably beyond the reach of our manybody truncation scheme. As a matter of fact, improving our description of the lowlying spectroscopy of Al and/or calculating its groundstate quadrupole moment Heylen et al. (2016) constitutes a challenging task for theory. As for ab initio SCGF method, one possible way out consists in allowing for spherical symmetry to be spontaneously broken in the calculation. This would however require for good angular momentum to be eventually restored, which is yet to be formulated within the frame of the SCGF method.
All in all, the spectroscopy of nuclei obtained from Si and S by the addition (removal) of a neutron (proton) is in fair agreement with data at the present stage of development of ab initio calculations in midmass nuclei. Still, the significant disagreement with data in Al is at variance with the other cases and deserves more attention in the future. At this point in time, the comparison to such refined spectroscopic observables give further confidence into the predictions made about the occurrence of a bubble structure in Si.
iv.3 Spinorbit splitting and bubble structure
S  Si  SSi  

SCGF  2.18  1.16  1.02 (47) 
(d,p)  1.99  0.91  1.08 (54) 
As alluded to in the introduction, it has been speculated that the presence of a semibubble in the center of a light nucleus could feedback on the splitting between low angularmomentum spinorbit partners in the systems. Consequently, the question presently arises to which extent the appearance of a significant bubble () when removing two protons from S to Si impacts the size of the splitting between the lowlying and states when going from S to Si. As reported in Tab. 9, this splitting is predicted to decrease from MeV to MeV, i.e. to be reduced by MeV (47). This sudden reduction by about 50 is unique over the nuclear chart and is in good agreement with data both on an absolute and a relative scale.
To better establish whether there is an actual correlation between the appearance of the bubble in Si and the reduction of the splitting when going from S to Si, the analysis is first repeated with NN+3N400(1.88) for which the bubble structure is predicted to be less pronounced, i.e. instead of for NNLO at the ADC(2) level. The splitting between the and states only reduces by 27^{11}^{11}11The spectra being significantly more spread out with NN+3N400(1.88), it is not meaningful to compare the reduction of the splitting in absolute values. Still, let us mention that the splitting decreases by about keV out of the MeV value it takes in S, i.e. it actually decreases less in absolute value than with NNLO is spite of being MeV larger (too large) in the first place. when going from S to Si, which is indeed in proportion of of the reduction seen with NNLO. To establish the correlation on firmer a ground, the reduction of the splitting is plotted in the upper panel of Fig. 21 against the charge depletion factor in Si for all presently available SCGF calculations^{12}^{12}12The theoretical data set contains calculations performed with and without 3N interactions at the ADC(1), ADC(2) and ADC(3) levels for the four Hamiltonians under present consideration.. This systematic gives quantitative credit to the existence of a correlation between the appearance of a bubble structure and the reduction of spinorbit splittings of lowlying/low angularmomentum states.
Last but not least, the lower panel of Fig. 21 highlights the fact that the depletion factor of Si is significantly correlated with the difference of charge rms radius between S and Si. Furthermore, the sensitivity to the bubble leads to a variability of the order of fm between the case where there would be no bubble () in Si and our current prediction (), which is several times larger than the typical 1 uncertainty for the measurement of absolute charge radii whenever doable Fricke and Heilig (2004). The latter acts as a strong motivation to measure the absolute charge radius of Si or its isotopic shift with respect to Si whose charge radii are known Fricke and Heilig (2004).
iv.4 Effective singleparticle energies
The original speculation that the presence of a bubble could lead to a reduction of the splitting between spinorbit partners characterized by moderate angular momenta ToddRutel et al. (2004) was based on a meanfield picture in which the effective onebody spinorbit potential is essentially proportional to the derivative of the pointnucleon density distribution. We have seen above that this correlation does indeed manifest when looking at manybody energies in spite of the fact that no onebody spinorbit potential proportional to the derivative of the density is explicitly at play in the ab initio resolution of the manybody Schrödinger equation. It is thus of interest to reverse engineer and study to which extent a consistent picture does indeed emerge in the onebody shell structure that underlines the correlated manybody system. It is indeed useful to interpret the behavior of observable onenucleon separation energies in simple terms. One must however be clear that effective singleparticle energies are not observable, i.e. the picture provided by this interpretation is only valid within the theoretical scheme used and must not be extrapolated a priori to other theoretical schemes^{13}^{13}13A theoretical scheme is fully defined by the given of interacting degrees of freedom and of the Hamiltonian that drives their dynamics, without any further freedom to perform unitary transformations of that Hamiltonian. Indeed, the value of effective singleparticle energies can be changed within the theory through a unitary transformation without changing observable onenucleon separation energies Duguet et al. (2015). This demonstrates that both quantities possess different ontological status within the frame of manybody quantum mechanics and that a quantitative link between them can only be made within the theory by fixing the freedom allowed by unitary transformations. Duguet et al. (2015).
Having the complete set of onenucleon addition/removal energies and of spectroscopic probability matrices at hand, the centroid Hamiltonian can be constructed to give access, via its diagonalization, to the onenucleon shell structure associated with effective single particle energies (ESPEs) Baranger (1970). In the corresponding eigenbasis, ESPEs write as
(12) 
The set of ESPEs computed for Si and S appear in Figs. 17 and 18. By comparing ESPEs to the associated spectral strength distribution, one appreciates the fragmentation of the latter. Correspondingly, any given ESPE is not in onetoone correspondence with the dominant fragment of appropriate spin and parity but is a centroid of all the fragments with the same spin and parity. This is particularly striking in S where the ordering of neutron ESPEs above the Fermi level differs from the ordering of the main peaks in the spectral strength distribution of S. Indeed, the is above the whereas the main fragment is below the main fragment. The orderings of both spectra become consistent when removing two protons given that the strength is less fragmented in Si than in S.
S  Si  SSi  

SCGF  1.50  1.07  0.43 (29) 
SM      0.38 (25) 
In connection with the study performed in Sec. IV.1 on the splitting between lowlying and states reached by onenucleon addition processes, we now focus on the evolution of the neutron 1p1p ESPE spinorbit splitting when going from S to Si. As qualitatively visible in Figs. 17 and 18, and as quantitatively reported in Tab. 10, the ESPE spinorbit splitting reduces from MeV to MeV, i.e. it is lowered by . This feature, as well as the correlation this reduction of the ESPE splitting entertains with the charge depletion factor, is also illustrated in Fig. 21.
The reduction of the 1p1p ESPE spinorbit splitting obtained from of a full sdpf SM calculation Burgunder et al. (2014) is also reported in Tab. 10. Within that theoretical scheme, which reproduces by construction the experimental energies of the main fragments in S to Si once full correlations are included, the ESPE spinorbit splitting reduces by , which is close to the reduction obtained within the SCGF calculation based on the NNLO Hamiltonian. While ESPEs do not have to agree, even if both theoretical calculations equally reproduce observable manybody energies, the reduction of the ESPE spinorbit splitting happen to be similar in both schemes.
The analysis of the reduction of the ESPE splitting as a function of the depletion factor allows to better understand how the reduction of the observable manybody splitting operates in the present ab initio calculation. The splitting being a coherent combination of the ESPE contribution and of manybody correlations Duguet et al. (2015), one sees that the former contribution is responsible for about of the total reduction while the latter generates the other in the full calculation with NNLO. The contribution from manybody correlations relates to the fact that the strength is fragmented in S but not in Si. This is testified by the presence of the second state just above the state in the left panel of Fig. 18. As a result, the lowestlying state is mechanically pushed away from the state, thus adding to the splitting beyond the sole migration of the ESPE centroids. As the bubble disappears and the charge radius decreases in the lowest two panels of Fig. 21, both the relative migration of 1p and 1p ESPEs and the fragmentation of the strength (not shown here) diminish, thus leading to a less pronounced reduction of the splitting when going from S to Si.
V Conclusions
Semibubble or bubble structures have been invoked in connection with hypothetical superheavy or hyperheavy nuclei as a result of a collective quantum mechanical effect, further sustained by the compromise between the large repulsive Coulomb interaction and the strong force that binds them. In lighter systems, semibubble structures have been conjectured on the basis of the sole quantum mechanical effect, which finds its source in the sequence of occupied and unoccupied singleparticle states near the Fermi energy in an independentparticle or a meanfield picture.
In connection with the latter category, we have studied the potential bubble Si nucleus on the basis of stateoftheart ab initio quantum manybody methods. The occurrence of a depletion in the center of the charge density distribution and its correlation with an anomalously weak splitting between low angularmomentum spinorbit partners was originally postulated on the basis of a meanfield rationale. It was thus of interest to investigate this speculation on the basis of fully correlated solution of the manybody Schrödinger equation.
We have performed ab initio selfconsistent Green’s function manybody calculations of Si and S. The focus was put on binding energies, charge density distributions and charge root mean square radii of Si and S as well as on the lowlying spectroscopy of nuclei obtained via the addition (removal) of a neutron (proton). Predictions regarding the (non)existence of the bubble structure in Si have been shown to vary significantly with the input chiral effective field theory Hamiltonian. However, demanding that the experimental charge density distribution and the root mean square radius of S, as well as Si and S binding energies, are well reproduced has only left the NNLO Hamiltonian as a reliable candidate to perform these predictions. Additionally, the role of the threenucleon interaction in the generation of the bubble structure has been scrutinized.
Employing the NNLO Hamiltonian, a bubble structure has been convincingly predicted in Si. In particular, it has been demonstrated that the depletion in the center of the charge density distribution should be large enough to leave a visible signal in the form factor to be measured in a future electron scattering experiment. Furthermore, a clear correlation has been established between the occurrence of the bubble structure and the weakening of the splitting in the spectrum of Si. Consequently, the latter does offer an indirect hint for the existence of the bubble in Si Mutschler et al. (2016b). The interpretation in terms of the underlying onenucleon shell structure has also been provided.
Present results will have to be revisited with future generations of EFT(based) interactions. Indeed, the use of binding energies and radii of nuclei containing about nucleons to adjust the lowenergy constants of two and threebody interactions entering the NNLO Hamiltonian leaves many questions open regarding epistemology aspects of ab initio calculations and their associated predictive power.
Present results provide a strong motivation to measure the charge density distribution of Si in future electron scattering experiments on unstable nuclei. It was also shown that the measurement of its charge radius would already bring key insight. Furthermore, it is of interest to perform oneneutron removal experiments on Si and S in order to further test our theoretical spectral strength distributions over a wide energy range Mutschler (2015). Last but not least, other bubble nuclei candidates, possibly in excited states, should be investigated theoretically and experimentally in the future.
Acknowledgements
The authors warmly thank B. Bally, M. Bender, A. Brown, T. Cocolios, G. Hagen, V. Lapoux, M. Martini, W. Nazarewicz, G. Neyens, R. F. Garcia Ruiz, O. Sorlin and P. van Duppen for useful discussions and comments. Calculations were performed by using HPC resources from GENCITGCC (Contract No. 2016057392) and at the DiRAC Complexity system at the University of Leicester (BIS National Einfrastructure capital grant No. ST/K000373/1 and STFC grant No. ST/K0003259/1). This work was supported by the United Kingdom Science and Technology Facilities Council (STFC) under Grant No. ST/L005816/1 and in part by the NSERC Grant No. SAPIN201600033. TRIUMF receives federal funding via a contribution agreement with the National Research Council of Canada.
References
 Dechargé et al. (2003) J. Dechargé, J.F. Berger, M. Girod, and K. Dietrich, Nucl. Phys. A 716, 55 (2003).
 Bender et al. (1999) M. Bender, K. Rutz, P.G. Reinhard, J. A. Maruhn, and W. Greiner, Phys. Rev. C 60, 034304 (1999).
 Bender and Heenen (2013) M. Bender and P.H. Heenen, J. Phys. Conf. Series 420, 012002 (2013).
 ToddRutel et al. (2004) B. G. ToddRutel, J. Piekarewicz, and P. D. Cottle, Phys. Rev. C 69, 021301 (2004).
 Burgunder (2012) G. Burgunder, ‘‘Etude de l’intéraction nucléaire spinorbite par réactions de transfert S(d,p) S et Si(d,p) Si,” (2012), Ph.D. Thesis ; https://tel.archivesouvertes.fr/tel00695010.
 Richter and Brown (2003) W. A. Richter and B. A. Brown, Phys. Rev. C 67, 034317 (2003).
 Khan et al. (2008) E. Khan, M. Grasso, J. Margueron, and N. Van Giai, Nucl. Phys. A 800, 37 (2008).
 Grasso et al. (2009) M. Grasso, L. Gaudefroy, E. Khan, T. Niksic, J. Piekarewicz, O. Sorlin, N. V. Giai, and D. Vretenar, Phys. Rev. C 79, 034318 (2009).
 Yao et al. (2012) J.M. Yao, S. Baroni, M. Bender, and P.H. Heenen, Phys. Rev. C 86, 014310 (2012).
 Yao et al. (2013) J. Yao, H. Mei, and Z. Li, Physics Letters B 723, 459 (2013).
 Wu et al. (2014) X. Y. Wu, J. M. Yao, and Z. P. Li, Phys. Rev. C 89, 017304 (2014).
 Ibbotson et al. (1998) R. W. Ibbotson, T. Glasmacher, B. A. Brown, L. Chen, M. J. Chromik, P. D. Cottle, M. Fauerbach, K. W. Kemper, D. J. Morrissey, H. Scheit, and M. Thoennessen, Phys. Rev. Lett. 80, 2081 (1998).
 Sorlin and Porquet (2008) O. Sorlin and M.G. Porquet, Prog. Part. Nucl. Phys. 61, 602 (2008).
 Rotaru et al. (2012) F. Rotaru, F. Negoita, S. Grévy, J. Mrazek, S. Lukyanov, F. Nowacki, A. Poves, O. Sorlin, C. Borcea, R. Borcea, A. Buta, L. Cáceres, S. Calinescu, R. Chevrier, Z. Dombrádi, J. M. Daugas, D. Lebhertz, Y. Penionzhkevich, C. Petrone, D. Sohler, M. Stanoiu, and J. C. Thomas, Phys. Rev. Lett. 109, 092503 (2012).
 Baumann et al. (1989) P. Baumann, A. Huck, G. Klotz, A. Knipper, G. Walter, G. Marguier, H. L. Ravn, C. RichardSerre, A. Poves, and J. Retamosa, Phys. Lett. B 228, 458 (1989).
 Miessen (1982) H. Miessen, “Bestimmung der Magnetisierungsdichte von Si und P durch Elektronenstreuung,” (1982), Ph.D. Thesis; University of Mainz, Germany.
 Simon (2007) H. Simon, Nucl. Phys. A 787, 102 (2007).
 Suda et al. (2009) T. Suda et al., Phys. Rev. Lett. 102, 102501 (2009).
 Thorn et al. (1984) C. E. Thorn, J. W. Olness, E. K. Warburton, and S. Raman, Phys. Rev. C 30, 1442 (1984).
 Eckle et al. (1989) G. Eckle, H. Kader, H. Clement, F. J. Eckle, F. Merz, R. Hertenberger, H. J. Maier, P. Schiemenz, and G. Graw, Nucl. Phys. A 491, 205 (1989).
 Burgunder et al. (2014) G. Burgunder, O. Sorlin, F. Nowacki, S. Giron, F. Hammache, M. Moukaddam, N. de Séréville, D. Beaumel, L. Càceres, E. Clément, G. Duchêne, J. P. Ebran, B. FernandezDominguez, F. Flavigny, S. Franchoo, J. Gibelin, A. Gillibert, S. Grévy, J. Guillot, A. Lepailleur, I. Matea, A. Matta, L. Nalpas, A. Obertelli, T. Otsuka, J. Pancin, A. Poves, R. Raabe, J. A. Scarpaci, I. Stefan, C. Stodel, T. Suzuki, and J. C. Thomas, Phys. Rev. Lett. 112, 042502 (2014).
 Khan et al. (1985) S. Khan, T. Kihm, K. T. Knopfle, G. Mairle, V. Bechtold, and L. Freidrich, Phys. Lett. 156B, 155 (1985).
 Mutschler et al. (2016a) A. Mutschler, O. Sorlin, A. Lemasson, D. Bazin, C. Borcea, A. Gade, H. Iwasaki, E. Khan, A. Lepailleur, F. Recchia, T. Roger, F. Rotaru, M. Stanoiu, S. R. Stroberg, J. A. Tostevin, M. Vandebrouck, D. Weisshaar, and K. Wimmer, Phys. Rev. C 93, 034333 (2016a).
 Mutschler et al. (2016b) A. Mutschler, A. Lemasson, O. Sorlin, D. Bazin, C. Borcea, Z. Dombrádi, J.P. Ebran, A. Gade, H. Iwasaki, E. Khan, A. Lepailleur, F. Recchia, T. Roger, F. Rotaru, D. Sohler, M. Stanoiu, S. R. Stroberg, J. A. Tostevin, M. Vandebrouck, D. Weisshaar, and K. Wimmer, Nature Phys. , in press (2016b).
 Dickhoff and Barbieri (2004) W. H. Dickhoff and C. Barbieri, Prog. Part. Nucl. Phys. 52, 377 (2004).
 Barbieri and HjorthJensen (2009) C. Barbieri and M. HjorthJensen, Phys. Rev. C 79, 064313 (2009).
 Somà et al. (2011) V. Somà, T. Duguet, and C. Barbieri, Phys. Rev. C 84, 064317 (2011).
 Ekström et al. (2015) A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen, T. Papenbrock, B. D. Carlsson, C. Forssén, M. HjorthJensen, P. Navrátil, and W. Nazarewicz, Phys. Rev. C 91, 051301 (2015).
 Ekström et al. (2013) A. Ekström, G. Baardsen, C. Forssén, G. Hagen, M. HjorthJensen, G. R. Jansen, R. Machleidt, W. Nazarewicz, T. Papenbrock, J. Sarich, and S. M. Wild, Phys. Rev. Lett. 110, 192502 (2013).
 Entem and Machleidt (2003) D. R. Entem and R. Machleidt, Phys. Rev. C 68, 041001 (2003).
 Navrátil (2007) P. Navrátil, FewBody Systems 41, 117 (2007).
 Roth et al. (2012) R. Roth, S. Binder, K. Vobig, A. Calci, J. Langhammer, and P. Navrátil, Phys. Rev. Lett. 109, 052501 (2012).
 Bogner et al. (2010) S. K. Bogner, R. J. Furnstahl, and A. Schwenk, Prog. Part. Nucl. Phys. 65, 94 (2010).
 Somà et al. (2014a) V. Somà, A. Cipollone, C. Barbieri, P. Navrátil, and T. Duguet, Phys. Rev. C 89, 061301(R) (2014a).
 Hergert et al. (2014) H. Hergert, S. Bogner, T. Morris, S. Binder, A. Calci, et al., Phys. Rev. C 90, 041302 (2014).
 Lapoux et al. (2016) V. Lapoux, V. Somà, C. Barbieri, H. Hergert, J. D. Holt, and S. R. Stroberg, Phys. Rev. Lett. 117, 052501 (2016).
 Carbone et al. (2014) A. Carbone, A. Rios, and A. Polls, Phys. Rev. C 90, 054322 (2014).
 Hagen et al. (2014) G. Hagen, T. Papenbrock, A. Ekström, K. A. Wendt, G. Baardsen, S. Gandolfi, M. HjorthJensen, and C. J. Horowitz, Phys. Rev. C 89, 014319 (2014).
 Barbieri (2009) C. Barbieri, Phys. Rev. Lett. 103, 202502 (2009).
 Cipollone et al. (2013) A. Cipollone, C. Barbieri, and P. Navrátil, Phys. Rev. Lett. 111, 062501 (2013).
 Somà et al. (2014b) V. Somà, C. Barbieri, and T. Duguet, Phys. Rev. C 89, 024323 (2014b).
 Schirmer et al. (1983) J. Schirmer, L. S. Cederbaum, and O. Walter, Phys. Rev. A 28, 1237 (1983).
 Shin et al. (2016) I. J. Shin, Y. Kim, P. Maris, J. P. Vary, C. Forssén, J. Rotureau, and N. Michel, (2016), arXiv:1605.02819 [nuclth] .
 Audi et al. (2012) G. Audi, M. Wang, A. H. Wapstra, F. G. Kondev, M. MacCormick, X. Xu, and B. Pfeiffer, Chinese Phys. C 36, 1287 (2012).
 Cipollone et al. (2015) A. Cipollone, C. Barbieri, and P. Navrátil, Phys. Rev. C 92, 014306 (2015).
 Garcia Ruiz et al. (2016) R. F. Garcia Ruiz et al., Nature Phys. 12, 594 (2016).
 Campbell et al. (2016) P. Campbell, I. Moore, and M. Pearson, Progress in Particle and Nuclear Physics 86, 127 (2016).
 Garcia Ruiz (2016) R. F. Garcia Ruiz, (2016), private communication.
 Angeli and Marinova (2013) I. Angeli and K. Marinova, Atomic Data and Nuclear Data Tables 99, 69 (2013).
 Lapoux and Alamanos (2015) V. Lapoux and N. Alamanos, Eur. Phys. J. A 51, 91 (2015).
 Maréchal et al. (1999) F. Maréchal, T. Suomijärvi, Y. Blumenfeld, A. Azhari, E. Bauge, D. Bazin, J. A. Brown, P. D. Cottle, J. P. Delaroche, M. Fauerbach, M. Girod, T. Glasmacher, S. E. Hirzebruch, J. K. Jewell, J. H. Kelley, K. W. Kemper, P. F. Mantica, D. J. Morrissey, L. A. Riley, J. A. Scarpaci, H. Scheit, and M. Steiner, Phys. Rev. C 60, 034615 (1999).
 Navrátil (2004) P. Navrátil, Phys. Rev. C 70, 014317 (2004).
 Hagen et al. (2009) G. Hagen, T. Papenbrock, and D. J. Dean, Phys. Rev. Lett. 103, 062503 (2009).
 Duguet et al. (2015) T. Duguet, H. Hergert, J. D. Holt, and V. Somà, Phys. Rev. C 92, 034313 (2015).
 Furnstahl and Hammer (2002) R. J. Furnstahl and H. W. Hammer, Phys. Lett. B 531, 203 (2002).
 Chandra and Sauer (1976) H. Chandra and G. Sauer, Phys. Rev. C 13, 245 (1976).
 Mohr et al. (2012) P. J. Mohr, B. N. Taylor, and D. B. Newell, Rev. Mod. Phys. 84, 1527 (2012).
 Pohl et al. (2010) R. Pohl et al., Nature 466, 213 (2010).
 Hergert et al. (2016) H. Hergert, S. K. Bogner, T. D. Morris, A. Schwenk, and K. Tsukiyama, Phys. Rep. 621, 165 (2016).
 Negele (1970) J. W. Negele, Phys. Rev. C 1, 1260 (1970).
 Rychel et al. (1983) D. Rychel, H. J. Emrich, H. Miska, R. Gyufko, and C. A. Wiedner, Phys. Lett. 130B, 5 (1983).
 Carlsson et al. (2016) B. D. Carlsson, A. Ekstrom, C. Forssen, D. F. Stromberg, G. R. Jansen, O. Lilja, M. Lindby, B. A. Mattsson, and K. A. Wendt, Phys. Rev. X 6, 011019 (2016).
 Bedaque and van Kolck (2002) P. F. Bedaque and U. van Kolck, Ann. Rev. Nucl. Part. Sci. 52, 339 (2002).
 Nogga et al. (2005) A. Nogga, R. G. E. Timmermans, and U. v. Kolck, Phys. Rev. C 72, 054006 (2005).
 Epelbaum et al. (2009) E. Epelbaum, H.W. Hammer, and U.G. Meissner, Rev. Mod. Phys. 81, 1773 (2009).
 Krebs et al. (2016) H. Krebs, E. Epelbaum, and U. G. MeiÃner, (2016), arXiv:1610.03569 [nuclth] .
 Zhao et al. (2010) P. W. Zhao, Z. P. Li, J. M. Yao, and J. Meng, Phys. Rev. C 82, 054319 (2010).
 Brown and Wildenthal (1988) B. A. Brown and B. H. Wildenthal, Annu. Rev. Nucl. Part. Sci. 38, 29 (1988).
 Brown (2016) B. A. Brown, (2016), private communication.
 Papuga et al. (2014) J. Papuga, M. L. Bissell, K. Kreim, C. Barbieri, K. Blaum, M. D. Rydt, T. Duguet, R. F. G. Ruiz, H. Heylen, M. Kowalska, R. Neugart, G. Neyens, W. Nortershauser, M. M. Rajabali, R. Sánchez, N. Smirnova, V. Somà, and D. T. Yordanov, Phys. Rev. C 90, 034321 (2014).
 Heylen et al. (2016) H. Heylen, M. D. Rydt, G. Neyens, M. L. Bissell, L. Caceres, J. M. Daugas, Y. Ichikawa, Y. Ishibashi, O. Kamalou, T. J. Mertzimekis, P. Morel, J. Papuga, A. Poves, M. M. Rajabali, C. Stodel, J. C. Thomas, H. Ueno, Y. Utsuno, N. Yoshida, and A. Yoshimi, Phys. Rev. C 94, 034312 (2016).
 Fricke and Heilig (2004) G. Fricke and K. Heilig, Nuclear Charge Radii (SpringerVerlag, 2004).
 Baranger (1970) M. Baranger, Nucl. Phys. A 149, 225 (1970).
 Mutschler (2015) A. Mutschler, “Le noyaubulle de Si : Un outil expérimental pour étudier l’interaction spinorbite?” (2015), Ph.D. Thesis; https://tel.archivesouvertes.fr/tel01206188.