Skyrme functional with tensor terms from ab initio calculations
of neutronproton drops
Abstract
A new Skyrme functional devised to account well for standard nuclear properties as well as for spin and spinisospin properties is presented. The main novelty of this work relies on the introduction of tensor terms guided by ab initio relativistic BruecknerHartreeFock calculations of neutronproton drops. The inclusion of tensor term does not decrease the accuracy in describing bulk properties of nuclei, experimental data of some selected spherical nuclei such as binding energies, charge radii, and spinorbit splittings can be well fitted. The new functional is applied to the investigation of various collective excitations such as the Giant Monopole Resonance (GMR), the Isovector Giant Dipole Resonance (IVGDR), the GamowTeller Resonance (GTR), and the SpinDipole Resonance (SDR). The overall description with the new functional is satisfactory and the tensor terms are shown to be important particularly for the improvement of the SpinDipole Resonance results. Predictions for the neutron skin thickness based on the nonenergy weighted sum rule of the SpinDipole Resonance are also given.
UTF8gbsn
I Introduction
Nuclear energy density functionals such as Skyrme, Gogny, and relativistic ones have achieved great success since a few decades Vautherin and Brink (1972); Dechargé and Gogny (1980); Bender and Heenen (2003); Meng (2016). Actually, by fitting to a reduced set of experimental data (mainly binding energies and charge radii), the obtained functionals can give very good description of groundstate as well as excited state properties (such as the excitation energies of Giant Resonances) along the whole nuclear chart. On the other hand, it should also be noted that there are still many open questions regarding the existing functionals. For example, certain properties are difficult to be tightly constrained from experimental data: one example is the symmetry energy Baldo and Burgio (2016); RocaMaza et al. (2018); RocaMaza and Paar (2018), that is at the same time extremely relevant to describe exotic neutronrich nuclei or neutron stars. Another key point is that there is not a unique ansatz to write down the structure of a functional; neither one knows clearly which is the path for the systematic improvement of existing functionals. In the context of the present paper, it is relevant to stress that the role of tensor terms within nuclear DFT is not completely elucidated despite significant efforts and partial achievements Sagawa and Colò (2014).
It is well known that the tensor force is an important component of the bare nucleonnucleon interaction Brown et al. (2010), but its relevance for energy density functionals is another question. In earlier implementations of DFT, based on effective Hamiltonians, the tensor force was either deemed not important or often ignored for the sake of simplicity. Nowadays, the situation has changed in particular for what concerns the Skyrme and Gogny functionals (the discussion about relativistic functionals is outside our scope here but the reader can consult Refs. Sagawa and Colò (2014); Meng (2016)). Skyrme functionals have the special property that the contribution of the tensor force to the functional goes together with that from the exchange terms of the central force, if we start from an effective interaction (or pseudopotential). These contributions are called terms of the functional (see Sec. II, where we define precisely the socalled spinorbit density ). Since these terms vanish in spinorbit saturated nuclei, and many nuclei used in the fits belong to this category, one can easily understand why not only the tensor terms but the whole terms were ignored in most of the early Skyrme parametrizations (see, for example, the discussion in Ref. Lesinski et al. (2007)).
The revival of the discussion on tensor terms occurred in the beginning of the 21st century. With the development of new radioactive ion beam facilities, more and more exotic neutronrich nuclei have been studied and new phenomena regarding the evolution of the spinorbit splittings with neutron excess have been highlighted Schiffer et al. (2004). Soon afterwards, Otsuka et al. suggested the importance of an effective neutronproton tensor force to explain this evolution Otsuka et al. (2005). A number of experimental data could be reproduced, to a certain extent better by including tensor terms, within the Gogny Otsuka et al. (2006), Skyrme Brown et al. (2006); Colò et al. (2007); Brink and Stancu (2007); Lesinski et al. (2007), and relativistic frameworks Long et al. (2008) (cf. also Refs. Jiang et al. (2015); Li et al. (2016); Karakatsanis et al. (2017); Brink and Stancu (2018); LópezQuelle et al. (2018); Zong and Sun (2018); Wang et al. (2018)).
In general, however, it is difficult to single out observables that can uniquely pinpoint the effects of the tensor force. For example, if one looks at singleparticle data, it is well known that particlevibration coupling (PVC) is also important for singleparticle properties Colò et al. (1994); Niu et al. (2012); Litvinova and Afanasjev (2011); Afanasjev and Litvinova (2015); as a consequence, trying to determine the strength of the tensor terms by looking at experimental singleparticle properties without including PVC effects in the theory, may not be well justified. This may explain the difficulties found in the systematic study carried out in Ref. Lesinski et al. (2007), where 36 Skyrme functionals were produced, characterized by different strengths of tensor terms, and properties like singleparticle levels, spinorbit splittings, mass residuals and anomalies of radii were compared with experimental data. It has been found very hard to satisfy all the experimental constraints, so that the Skyrme ansatz itself has been blamed for this deficiency.
In such a situation, ab initio calculations should become a benchmark. Based on the use of realistic interactions which can describe twonucleon scattering phase shifts and deuteron properties, ab initio calculations have no free parameters and can provide in some cases a reliable connection between fewbody and manybody physics Drut et al. (2010). The obtained results can then be used as pseudodata, in order to complement the experimental data and calibrate nuclear density functionals. For example, pseudodata obtained within the BruecknerHartreeFock (BHF) theory Day (1967) do not contain effects like particlevibration coupling, and can be used directly as a benchmark for DFT calculations.
In the case of electronic systems, the derivation of energy density functionals from ab initio calculations has been developed for a long time and has achieved many successes Perdew and Kurth (2003). Because of the more complicated nature of the nuclear interaction as compared with the Coulomb interaction, ab initio calculations in nuclear physics are much more difficult Dickhoff and Barbieri (2004); Lee (2009); Liu et al. (2012); Barrett et al. (2013); Hagen et al. (2014); Carlson et al. (2015); Hergert et al. (2016); Shen et al. (2016, 2017) and fall behind that in atomic physics. In fact, the study of the nucleonnucleon interaction itself, especially its derivation from the QCDinspired chiral Effective Field Theory (EFT) vs. the traditional mesonexchange force, or the precise assessment of threebody forces, is a very active and important field of research where the community is still making progress now Epelbaum et al. (2009); Machleidt and Entem (2011). Thus, despite the difficulties, one can say that the nuclear ab initio calculations can provide valuable information for the derivation of nuclear energy density functionals Drut et al. (2010); Liang et al. (2018).
In this context, recently, studies on the neutron drops are receiving more and more attention Pudliner et al. (1996); Gandolfi et al. (2011); Maris et al. (2013); Potter et al. (2014); Tews et al. (2016); Zhao and Gandolfi (2016); Shen et al. (2018a, b); Bonnard et al. (2018). The neutron drop is an ideal system composed of a finite number of neutrons confined in an external field. Similarly as nuclear matter, which is also an ideal and relatively easy system due to translational invariance that yet can provide us with valuable information RocaMaza and Paar (2018), neutron drops are also simple (only neutronneutron interaction, no centerofmass correction because of the external field) and can give an interesting insight on shell structure and/or finitesize properties, which is not the case for nuclear matter. Neutron systems are also more favorable for ab initio treatment, due also to the large value of the neutronneutron scattering length that makes the system close to the unitary limit Braaten and Hammer (2006).
For example, by comparing the ab initio quantum MonteCarlo (QMC) calculations with Skyrme HartreeFock calculations, suggestions for the improvement of future Skyrme energy density functionals have been given concerning the density distribution, the spinorbit splittings, the isovector gradient terms of the functionals, and the pairing sector Pudliner et al. (1996); Gandolfi et al. (2011). In another work, the properties of periodically modulated neutron matter were studied by the quantum MonteCarlo method, and this was used to constrain the isovector term of Skyrme functionals Buraczynski and Gezerlis (2016, 2017). In Ref. Bonnard et al. (2018), the ab initio calculations for neutron drops are used to further constrain previous EFTinspired energydensity functionals in nuclear matter, by comparing quantities like the total energies, singleparticle potentials, and density distributions. In a recent work Shen et al. (2018a, b), the effect of the tensor force in neutron drops has been studied by the ab initio relativistic BruecknerHartreeFock (RBHF) theory with the Bonn A interaction Machleidt (1989). A clear evidence of the impact of tensor force on the evolution of spinorbit splittings has been found, and this provides a guide for determining the strength of the tensor force in the nuclear medium Shen et al. (2018a, b).
In these works Shen et al. (2018a, b), only information on the tensor force between neutrons is present, but not information on the interaction between neutrons and protons. This information is essential for a complete study of the tensor force in nuclei, where all the previous studies we have mentioned point to the dominance of the neutronproton component. Therefore, in this work we include the protons into the neutron drops, that is, we study the neutronproton drops with the RBHF theory. Obviously, we are not comparing with any experimental data, but conceptually we deem we can use the results on this ideal system as pseudodata to provide reliable information to constrain a nuclear energy density functional. In this work we aim at fitting a Skyrme energy density functional. We keep this neutronproton system as simple as possible: we do not include the Coulomb interaction between the protons, and neutron and proton masses are the same.
There are two main reasons why we choose to study artificial trapped neutronproton drops instead of realistic finite nuclei. First, concerning the purpose of this work, there is no difference between extracting information on the tensor force from finite nuclei and from artificial neutronproton drops. But, without centerofmass correction and Coulomb interaction, the neutronproton drops are easier to be calculated within RBHF. Secondly, following the spirit of using ideal systems such as nuclear matter and neutrondrops to study different aspects of nuclear systems, the neutronproton drops provide another ideal case to investigate different properties with more freedom than finite nuclei: for instance, one can vary the neutron and proton numbers without worrying if the system is bound or not.
At variance with previous works in which the tensor force is added perturbatively on top of existing Skyrme functionals, such as in Ref. Colò et al. (2007), the tensor force in this work is included selfconsistently with a refit of the whole parameter set, like in the case of the TIJ family Lesinski et al. (2007). But differently from the case of the TIJ family, the tensor terms in this work will be fitted using the pseudodata of the neutronproton drops coming from RBHF calculations. Of course, such pseudodata does not necessarily need to come from RBHF calculations or from a specific ab initio calculation, relativistic or nonrelativistic. In principle, relativistic effects such as nucleonantinucleon excitations can be expressed in nonrelativistic framework in terms of a threebody force Brown et al. (1987), and it has been shown that the results of RBHF using the onebosonexchange Bonn interaction does agree well with the results of nonrelativistic BHF using twobody and threebody forces Sammarruca et al. (2012). For the fit as a whole, we take inspiration from the successful fitting protocol of SAMi RocaMaza et al. (2012, 2013a), which has been shown to display satisfactory results for spin and spinisospin resonances. The newly established Skyrme functional is named SAMiT.
In Sec. II, we give a brief summary of the Skyrme interaction with tensor terms, and HartreeFock equations. The numerical details for the fitting are discussed in Sec. III. Results with the new interaction for nuclear matter and finite nuclei are presented in Sec. IV. Finally, the summary and perspectives for future investigations will be given in Sec. V.
Ii Theoretical Framework
ii.1 Relativistic BruecknerHartreeFock theory
In this Subsection, we will outline the theoretical framework of the RBHF for finite nuclear system. For more detail one can see Refs. Shen et al. (2016, 2017).
We start with a relativistic onebosonexchange interaction which describes the scattering data Machleidt (1989). The Hamiltonian can be expressed as:
(1) 
where the relativistic matrix elements are given by
(2)  
(3) 
The indices run over a complete basis of Dirac spinors with positive and negative energies. The twobody interaction contains the contributions from different mesons labelled by : scalar , vector , and pseudovector . For the isovector mesons and additional isospin matrices are included. In the Bonn interaction Machleidt (1989), a form factor of monopoletype is attached to each vertex and is the meson propagator. Retardation effects were deemed to be small and omitted from the beginning.
In Brueckner’s theory Brueckner et al. (1954); Day (1967), the effective interaction in the nuclear medium, the matrix, is used instead of the bare nucleonnucleon interaction, which has a strong repulsive core and is difficult to be used directly in nuclear manybody theory. The matrix takes into account the shortrange correlations by summing up all the ladder diagrams of the bare interaction and it is obtained by solving the BetheGoldstone equation,
(4) 
where in the RBHF theory are eigenstates of the relativistic HartreeFock equations with the corresponding singleparticle energies, are the antisymmetrized twobody matrix elements (II.1). The intermediate states run over all states above the Fermi surface with . In the above expression, is the starting energy Baranger (1969); Davies et al. (1969); Shen et al. (2017).
The singleparticle motion fulfills the relativistic HartreeFock equation in the external field of a harmonic oscillator (HO):
(5) 
with HO potential
(6) 
where is the nucleon mass and is the HO strength. The selfconsistent singleparticle potential is calculated by the matrix
(7) 
where the index runs over the occupied states in the Fermi sea (nosea approximation).
In the end, the total energy is calculated as
(8) 
where is the expectation value of the external field.
ii.2 Skyrme density functional
In this Subsection, we will outline the theoretical framework of the Skyrme interaction and the SkyrmeHartreeFock theory. For a detailed description of the Skyrme HartreeFock theory and the formulas in spherical nuclei, we refer to Vautherin and Brink (1972).
The Skyrme effective interaction with twobody tensor force is written in the standard form as Vautherin and Brink (1972); Stancu et al. (1977).
(9)  
(10) 
where , is the hermitian conjugate of acting on the left. The spinexchange operator reads , and is the total nucleon density.
The HartreeFock equations can be obtained by the variational method once the total energy of the HartreeFock ground state is written down. The HartreeFock equations read
(11) 
where is the singleparticle energy, is the corresponding wave function, and labels neutrons (protons). The singleparticle potential is a sum of central, Coulomb and spinorbit terms,
(12) 
The detailed expression of the central and Coulomb terms can be found, e.g., in Ref. Chabanat et al. (1998). For the calculation of neutron(proton) drops, we add an external field, which in this work is chosen as the harmonic oscillator potential (6).
The spinorbit term reads Stancu et al. (1977); Reinhard et al. (1995); Sagawa and Colò (2014)
(13) 
where the spinorbit density Chabanat et al. (1998). One should note that, starting from Eq. (9) one would derive . We adopt here a more general form and release this constraint as in SAMi RocaMaza et al. (2012). The parameters and in Eq. (13) include contributions from the (exchange part of the) central force and from the tensor force,
(14) 
They are given by
(15a)  
(15b) 
In the end, the HartreeFock total energy can be calculated in the usual way and reads
(16) 
where is the expectation value of the singleparticle kinetic energy on the orbit , is the rearrangement term. The factor one half in front of the external field arises because the remaining contribution is included in , in full analogy with the kinetic energy which is also a onebody operator.
Iii Numerical details
In the new Skyrme functional established in this work, we include the so called term, including and from the central force, and from the tensor force, as shown in Eq. (14). The HartreeFock equation is solved in coordinate space with a spherical box having size fm and a radial mesh of fm. The fitting protocol is based on the successful SAMi functional RocaMaza et al. (2012), with further constraint on the tensor force from neutronproton drops by RBHF theory: more precisely, the pseudodata of neutron matter from ab initio calculations are replaced by the results of neutron drops from RBHF calculations. The set of data and pseudodata to be fitted are:

the binding energies and charge radii of Ca, Ca, Zr, Sn, and Pb;

the spinorbit splittings of the proton orbit in Ca, proton orbit in Zr, and proton orbit in Pb. This will mainly determine the spinorbit term and ;

the relative change of spinorbit splittings from neutronproton drops to calculated by RBHF theory using Bonn A interaction. Since only the relative change is fitted, it mainly determines the tensor term and but not the spinorbit term and ;
As mentioned in Ref. Shen et al. (2018a), the evolution of the spinorbit splittings by RBHF theory can be used to constrain the tensor force. On the other hand, the spinorbit splittings given by RBHF theory using Bonn A interaction are slightly smaller than experimental data Shen et al. (2016, 2017). For example, proton splitting of O given by RBHF is MeV, in comparison with experimental value, that is MeV. Therefore in our fit, we use the relative change of spinorbit splittings of neutronproton drops to constrain the tensor force, and the experimental data of finite nuclei to constrain the absolute strength of spinorbit splittings.
For the constraint on the isovector channel, we use the total energy of neutron drops calculated by RBHF instead of the total energy of neutron matter as in the case of SAMi. The results of RBHF are in reasonable agreement with other ab initio calculations and show some advantages: for the results of RBHF are similar to the results of AV8’ + IL7, but for RBHF gives more repulsion comparing with AV8’ + IL7 Shen et al. (2018a). This is favourable as AV8’ + IL7 can describe well the light nuclei system but tends to give too much attraction for pure neutron systems Sarsa et al. (2003) (see also green shadowed area and red squares in Fig. 2 and related text below). This is also more consistent for our purpose as the tensor force is constrained from the same RBHF calculations.
The numerical details of RBHF calculations can be found in Refs. Shen et al. (2017, 2018c, 2018b) and briefly summarized here. The calculations are performed in spherical box with box size fm. Notice the box size for RBHF calculations is considerably smaller than that of Skyrme functional calculations, which is 15 fm. In principle, a very large box would be the best choice, but in practice one needs to find a balance between the precision and computational cost. For RBHF calculations, the convergence properties with respect to the box size have been carefully checked in Ref. Shen et al. (2017) for mediumheavy nuclei, and in going from box size fm to fm, the energy per nucleon has been found to change less than MeV. With an additional external HO field, this convergence will be further improved as the singleparticle wave functions decay very fast when goes out of the size of the HO potential. Therefore, for our purpose here, dealing with nucleons at most in a HO trap (the rms radii for a neutron drop is fm), the box size fm is enough. As for Skyrme functional calculations, since heavy nuclei such as Pb are also to be studied, a larger box size is needed.
The initial singleparticle basis is chosen as a Dirac WoodsSaxon basis Zhou et al. (2003), and the energy cutoff is MeV for the positive energy states and MeV for the negative energy states (with a supplementary condition that at least 2 negativeenergy states are included in each block). The singleparticle angular momentum cutoff is 24 . For neutronproton drops with , the twoparticle coupled total angular momentum cutoff is 8 ; for , this is 10 . If not specified otherwise, in all the calculations of this work the external field is chosen as an harmonic oscillator field with MeV. For RBHF, it has been implemented as a vector external field.
The results of neutronproton drops calculated by RBHF theory are shown in Table 1, including the total energies, neutron rms radii, and spinorbit splittings. From the Table one can see the effect of tensor force on the evolution of spinorbit splittings as the number of neutron changes Otsuka et al. (2005); Shen et al. (2018a). For the system with , that is spin saturated ( orbit fully occupied), the tensor force makes no effect; for the system with , that is spin unsaturated (the orbit is occupied while the is empty), the tensor force produces the largest effect so that the spinorbit splittings decrease.
68.6  2.82  5.52  8.73  5.52  8.73  
162.1  3.02  2.91  5.82  2.85  4.93 
Iv Results and discussion
iv.1 Fitted results
In Table 2 we list the data and pseudodata to be fitted, together with the corresponding errors for the fit and the number of data points . The partial contributions to the total from each category are also listed. The error for total binding energy is chosen instead of binding energy per nucleon. This will give more weight to heavier nuclei. Again, as in the fitting protocol of SAMi RocaMaza et al. (2012), the errors of the pseudodata are chosen in such a way that these data are used as a guide for the new functional and shall not affect the accuracy of fitting to experimental data. From the Table, it can be seen that when including tensor, the description of the experimental binding energies and charge radii keeps similar accuracy as for other Skyrme functionals.
Ref.  
MeV  5  Wang et al. (2017)  
fm  5  Angeli and Marinova (2013)  
3  Zalewski et al. (2008) ^{2}^{2}2data compiled in table III of Ref. Zalewski et al. (2008), see caption of this table for further details.  
4  [Tab. 1]  
4  Shen et al. (2018a)  
The parameters and saturation properties of SAMiT functional are listed in Table 3, where the corresponding estimations of the standard deviation Bevington and Robinson (1992) are also given. For the detail of the covariance analysis, see Appendix A. Comparing with SAMi RocaMaza et al. (2012), the strengths of the spinorbit terms for SAMiT are larger (for SAMi, MeV fm, MeV fm). This is understandable as SAMi does not contain a genuine tensor force that is known to reduce the spinorbit splittings of spinunsaturated systems such as Zr or Pb that have been fitted in both functionals. Therefore, the new set including tensor needs larger spinorbit strength to reproduce the same data.
Value  Error  Value  

MeV fm  fm  
MeV fm  MeV  
MeV fm  
MeV fm  
MeV  
MeV  
MeV  
(fixed)  
(fixed)  
MeV fm  
MeV fm  
MeV fm  
MeV fm 
It is worth to discuss in a specific manner the values of the socalled LandauMigdal parameters and which are associated, respectively, with the spin and spinisospin particlehole (ph) interaction. If one only looks at the ground state properties, different sets with different values of and can achieve similar accuracy Cao et al. (2010). But if one wants to have at the same time a good description of the properties of spin and spinisospin excitations like the GamowTeller or spindipole resonance, there will be more restrictions RocaMaza et al. (2012). As suggested in Refs. Suzuki and Sakai (1999); Wakasa et al. (2005), and taken into account by the fit of SAMi, we ought to respect the relations . Here, we obtain slightly different values with respect to SAMi as we are attempting to include more terms in the functional (tensor terms) and to reconcile with different pseudodata (associated with neutronproton drops and neutron drops from RBHF calculations).
iv.2 Ground state
In Table 4 we give the binding energy and charge radius of several doubly magic spherical nuclei calculated by using SAMiT, in comparison with experimental data Wang et al. (2017); Angeli and Marinova (2013) when existing. In most cases the descriptions of binding energy and charge radius are accurate within . Considering that there are free parameters ( and have been fixed), this accuracy is comparable with other commonly seen functionals, nonrelativistic or relativistic ones. For example, as indicated by the fitted , the /data of SAMiT for binding energies and radii are and , respectively (Table 2), while those of SAMi are and with similar data sets RocaMaza et al. (2012).
El.  (MeV)  (MeV)  (fm)  (fm)  

O  16  127.78(33)  127.62  2.774(4)  2.699 
Ca  40  343.74(52)  342.05  3.477(3)  3.478 
48  415.32(50)  415.99  3.515(3)  3.477  
Ni  56  468.73(1.06)  483.99  3.784(4)  
68  591.27(56)  590.41  3.901(4)  
Zr  90  783.35(46)  783.89  4.263(3)  4.269 
Sn  100  812.91(1.18)  824.79  4.480(5)  
132  1100.80(54)  1102.85  4.714(4)  4.709  
Pb  208  1637.81(67)  1636.43  5.479(5)  5.501 
Figure 1 shows the Equation of State (EoS) for both neutron matter and symmetric nuclear matter calculated by the new fitted Skyrme functional SAMiT, in comparison with those of SAMi functional. For the neutron matter, we also show the results of different ab initio calculations: selfconsistent Green’s function method with chiral and forces up to NLO Drischler et al. (2016) (cyan shade); quantum MonteCarlo method with chiral and forces up to NLO Lynn et al. (2016) (green shade); loop expansion around the HartreeFock energy with chiral force up to NLO and force up to NLO Hebeler and Schwenk (2010) (brown shade), see also Figure 1 from Ref. Tews et al. (2017).
Since the fitting protocol of SAMiT is similar to that of SAMi, the symmetric matter equationofstate of the two functionals are almost identical. On the other hand, the symmetry energy at saturation density of SAMiT, MeV, is slightly larger than that of SAMi, MeV. Therefore, the equationofstate of the neutron matter given by SAMiT shows more repulsion, and this is also in better agreement with the constraints from different ab initio calculations. The slope of the symmetry energy at saturation density of SAMiT is MeV, also slightly larger than that of SAMi, MeV.
In Fig. 2, we show the total energy (in units of ) of neutron drops from to in a HO trap ( MeV) calculated by SAMiT, in comparison with results of SAMi, RBHF calculations using the Bonn A interaction Shen et al. (2018a), and QMC calculations using AV8’+UIX (upper bound of the shaded area) and AV8’+IL7 (lower bound of the shaded area) Gandolfi et al. (2011); Maris et al. (2013), QMC using local chiral interaction NLO with two cutoffs fm and fm Tews et al. (2016), coupledcluster theory using chiral interaction NLO Potter et al. (2014). The nocore shell model using same chiral interaction NLO (calculated up to ) coincide with those by coupledcluster theory Potter et al. (2014) and will not be plotted. As it can be expected from the EoS in Fig.1, the total energy of neutron drops given by SAMiT are generally larger than those of SAMi, but smaller than those of RBHF with Bonn A calculations. The error we adopted for the energy of neutron drops in the fitting, for , assumes a clearer meaning by looking at Fig. 2: in fact, the values differing from the results of RBHF with Bonn A by about 0.01 are seen to be acceptable for our purpose.
Comparing with other (nonrelativistic) ab initio calculations using chiral interactions, the energy given by SAMiT (or those by RBHF with Bonn A) are slightly lower, similar as in the comparison in Fig. 1.
To see the subshell structure in neutron drops, we show in Fig. 3 the twoneutron separation energies for the same calculations we have just mentioned. Since those calculations are in a HO external field, they all show clear shell structure for , and . Besides those, the subshell structure at , and a somewhat weaker one can also be seen. For these subshells, SAMiT also shows similar trends as the RBHF calculations. Here the results of both functionals is very similar, with some improvement of SAMiT for and somehow also in the region . The different peak position in this region, for RBHF at (and a minor one at ) and for SAMiT (or SAMi) at is due to different ordering and different occupations of the singleparticle levels: For RBHF, from to , the level is being occupied and the energy suddenly increases much; from to the level continues to be occupied and the energy increases less; from to , the level is being occupied before the level is filled at , therefore another small jump appears at . For SAMiT (or SAMi) functional, from to the level is firstly filled and from to the begins to be occupied. Therefore, the jump of the energy is smooth and the twoneutron separation energy keeps increasing from to . From to the level continues to be filled until fully occupied, the energy increases less and therefore the decreases smoothly.
Figure 4 shows the neutron and proton and spinorbit splittings of neutronproton drops, whose trends are used as pseudodata in the fitting procedure and are mainly responsible for constraining the tensor force. Again, we show the results given by SAMiT, and comparing with those of SAMi functional and the pseudodata of RBHF with Bonn A interaction. The results of SAMiT and SAMi have been shifted so that the first points () meet with those of RBHF. The shifts for SAMiT are MeV () and MeV (); while for SAMi they are MeV () and MeV (). This is expected as the absolute value of spinorbit splittings of SAMiT are fitted to experimental value and they are larger than those of RBHF, as explained in the numerical details. On the other hand, SAMiT reproduces the relative changes, which reflects the effect of tensor force, very well. For the relative change of spinorbit splittings, the results of SAMiT are exactly on top of those of RBHF.
In Refs. Otsuka et al. (2005); Shen et al. (2018a), these trends for the evolution of the spinorbit splitting are explained by the effect of tensor forces. For the spinsaturated system with , there is no tensor force contribution from the neutrons; for the spinunsaturated system with , where only one of the spinorbit partners is occupied () while the other is empty (), the tensor force reduces the spinorbit splittings dramatically. This is shown in the RBHF calculation with Bonn A in Fig. 4, where the left panels depict the tensor force effect in the neutronneutron channel and the right panels correspond to effect in the neutronproton channel. For pure neutron drops within the RBHF calculation, a similar phenomenon has been shown in Ref. Shen et al. (2018a).
In the Skyrme functional, the tensor force in the same particle channel (neutronneutron or protonproton) is controlled by the parameter and the tensor force in the different particle channel is controlled by in Eq. (13). Both and have two origins, one is from the exchange terms of the central interaction ( and ) and another one is from the tensor force ( and ) as is clearly shown in Eq. (14). In Table 5, we list the values of and for SAMiT and SAMi functionals. Without tensor force, the contribution to and in SAMi functional comes only from the central part. One finds that the value of in SAMi is slightly larger than that of SAMiT, and therefore SAMi shows similar trends but with larger magnitude in the left panels of Fig. 4, where only the neutronneutron channel matters. With only the contribution from the central force, which is also very much determined by properties like binding energies and charge radii, SAMi has less freedom to fit to the spinorbit splittings in Fig. 4. As a consequence, the value of in SAMi is much smaller than that of SAMiT, and SAMi does not show the trends of the right panels of Fig. 4, where the neutronproton channel matters.
In conclusion, the tensor terms in the new Skyrme functional are determined reasonably and, more importantly, without ambiguities. Had we put more emphasis on the evolution of the experimental singleparticle states, we would have faced the problem that the coupling with collective vibrations may play an important role there Litvinova and Afanasjev (2011); Afanasjev and Litvinova (2015). Instead, fitting to the pseudodata of ab initio RBHF calculations does not imply this ambiguity.
(MeV fm)  (MeV fm)  

SAMiT  
SAMi RocaMaza et al. (2012) 
As a side remark, we note that from Fig. 4 we can also understand why SAMi gives a smaller value for the spinorbit strength because of the lack of tensor terms: both the new functional and SAMi have been fit using the proton spinorbit splittings in neutron spinunsaturated system (Zr with , where the orbit is occupied and the orbit is empty; Pb with , where the orbit is occupied and the orbit is empty), so that the functional with tensor terms must have a stronger spinorbit strength (or ) as the inclusion of the tensor force will imply a decrease of the SO splittings.
Since the values of are similar in SAMiT and SAMi, we can expect the two functionals show similar trends in the evolution of spinorbit splittings in neutron drops, and this is shown in Fig. 5. Similarly as in Fig. 4, the results of SAMiT and SAMi have been shifted so that the first points meet with those of RBHF: in other words, we put emphasis here on the relative changes rather than the absolute values. The shifts for SAMiT are MeV (), MeV (), MeV (), MeV (); while for SAMi they are MeV (), MeV (), MeV (), MeV (). As mentioned in the numerical settings, the pseudodata of spinorbit splittings by RBHF are used only to constrain the tensor terms and , but not the absolute strength of spinorbit terms. The spinorbit term have been constrained by the experimental spinorbit splittings.
Also in this case, SAMiT shows better agreement with the results of RBHF calculations. With fitting to the trend of the spinorbit splitting in neutronproton drops given by RBHF (as shown in Fig. 4), the trend of the same splitting in neutron drops shown in Fig. 5 can be reproduced automatically. It should also be noted that in detail, certain deviations exist. Using the results of RBHF as a baseline, the decrease of spinorbit splittings by SAMiT is smaller around (), and getting stronger around ( slightly larger while slightly smaller), and even more around ( slightly larger while agree). Similar trend is found for SAMi. This may be understandable considering the fact that the tensor terms in the Bonn A interaction are finiterange oneboson exchange interactions, while the tensor terms in Skyrme functional (10) are zerorange interactions. Another cause may be that the difference of the density distributions can also significantly influence the change of spinorbit splittings. Overall, we do not expect an exact matching with RBHF results but just to achieve the best possible fit.
In Fig. 6, we show the values of and from different Skyrme functionals, including three categories:

Skyrme functionals in which the tensor terms are included on equal footing with other terms and an overall fit is performed (brown squares in the figure), such as Skxta and Skxtb Brown et al. (2006), the TIJ family (I,J from to ) Lesinski et al. (2007), SkPT, SLy4T and SkOT Zalewski et al. (2008), and the present SAMiT (red star in the figure).
Among different Skyrme functionals, the values of SAMiT are close to those of T43, Skxta, and SkO’. The tensor part of Skxta is determined by the nonrelativistic HKT matrix Hosaka et al. (1985), which reproduces the matrix elements obtained from the bare nucleonnucleon Paris interaction Lacombe et al. (1980). It is interesting to see that the tensor force given by the relativistic matrix of Bonn A interaction is similar to the nonrelativistic matrix of Paris interaction. For SkO’, the tensor force is not included and therefore all the contributions to and come from the central part . This place a tight constraint for the central force and the groundstate properties are sacrificed to some extent. For example, the binding energy of Pb given by SkO’ is MeV, while the experimental data is MeV and the one given by SAMiT is MeV (see Table 4).
iv.3 Excited states
In this Subsection, we investigate some of the excited state properties within the selfconsistent HF plus Random Phase Approximation (RPA) Coló et al. (2013) using the new fitted Skyrme functional SAMiT. The basis for the configuration space RPA calculations is chosen by including all occupied states, and unoccupied states with an angular momentum cutoff , and states with increasing values of the quantum number in each block.
First of all, we would like to confirm that SAMiT does perform as well as other functionals for noncharge exchange resonances. With this aim, we show in Figure 7, as representative examples, the strength function associated with the isoscalar giant monopole resonance (GMR) and the isovector giant dipole resonance (IVGDR): the results obtained with SAMiT are compared with results from the SAMi functional RocaMaza et al. (2012) and with the experimental centroid energies Youngblood et al. (1999); Ryezayeva et al. (2002). The strength function is defined as Coló et al. (2013)
(17) 
where is the multipole excitation operator with total angular momentum , is the ground state, and is the excited state with energy . The operators used in GMR and IVGDR calculations are
(18) 
respectively.
The centroid energy of GMR by SAMiT is MeV, in good agreement with the experimental value MeV; our calculation exhausts of the energy weighted sum rule (EWSR) between and MeV. The centroid energy of IVGDR by SAMiT is MeV, slightly larger than the experimental value MeV; we fulfil of the EWSR between and MeV. We have also calculated using SAMiT without tensor and the results are very similar to SAMiT with tensor, therefore we do not plot them in the figure. In this subsection, unless specified otherwise, all calculations without tensor refer to dropping the tensor in the residual interaction of the RPA calculation, but the tensor is kept in HartreeFock.
We now study two different types of nuclear collective excitation where the spin and spinisospin channels of the functional should play a crucial role. In Fig. 8, we show the GTR strength function in Ca, Zr, and Pb calculated by SAMiT with and without tensor, in comparison with experimental data Wakasa et al. (1997); Yako et al. (2009); Krasznahorkay et al. (2001); Akimune et al. (1995); Wakasa et al. (2012) and results of SAMi functional RocaMaza et al. (2012). The operators for the GTR calculation is
(19) 
where is the spin, is the isospin raising (lowering) operator. In all cases, the nonenergy weighted sum rule (NEWSR)
(20) 
is fulfilled by the theoretical calculations.
For Ca, the locations of lowenergy peak and highenergy peak given by SAMiT are MeV and MeV, and they are in good agreement with the experimental values: MeV and MeV. For Zr, the locations given by SAMiT are MeV and MeV, and they also agree well with the experimental values: MeV and MeV. In comparison with SAMi, the highenergy peaks in these two nuclei described by the new functional are similar, and the lowenergy peaks have been improved. For Pb, the peak location given by SAMiT is MeV, also in good agreement with the experimental value MeV.
In the GTR shown in Fig. 8, the tensor force shows some influence but not a huge one. The effects of tensor force, however, are more significant in the case of the SDR. In Ref. Bai et al. (2010), it has been found that tensor correlations have a unique, multipoledependent effect on the three SDR components: the state is being pushed to higher energy, while the and states are being pushed to lower energy. Quantitatively, the residual interaction matrix element in the is the largest, the matrix element for is the next largest, and the effect on is rather small Bai et al. (2010).
In Fig. 9, we show the SDR strength function in the channel for Pb ( and total). We compare, once again, the results obtained by SAMiT with and without tensor with experimental data Wakasa et al. (2012) and the results with the SAMi functional RocaMaza et al. (2012). The operator used for the RPA calculations is
(21) 
Without tensor force, the new functional gives similar results as SAMi. It can be seen that the effect of tensor force is consistent with what has been found in Ref. Bai et al. (2010), namely different channels show different effects: For the , the tensor force improves the description for the data while for does not, and for the influence of the tensor force is very small. As the channel gives the largest contribution to the total SDR, the overall agreement with the experimental data is improved by including the tensor force. The experimental (SAMiT) NEWSR is ( fm) for , ( fm) for , ( fm) for , and fm ( fm) for the total. The NEWSR given by SAMiT is similar compared with SAMi RocaMaza et al. (2012).
The sum rule,
(22) 
is fully exhausted in the calculation.
As shown by the above formula, the sum rule of SDR is related with the neutron and proton rms radius. Since the sum rule given by SAMiT can well reproduce the experimental data, and the proton (charge) radius of Pb is well fitted by SAMiT (see Table 4), it would be interesting to show the neutron skin too. The neutron skin of Pb calculated by SAMiT is fm, to be compared with some recent experimental data fm (from parity violation Abrahamyan et al. (2012)), fm (from measurements with antiprotonic atoms Kłos et al. (2007)), fm (pion photoproduction Tarbert et al. (2014)), fm (electric dipole polarizability RocaMaza et al. (2013b)).
In Refs. Bai et al. (2010, 2011), different strengths of tensor terms have been proposed, on top of the SLy5 functional Chabanat et al. (1998). In connection with the SDR study, it was found in order to improve the strength function while not make the strength function worse, a relatively large value of and is needed (in Ref. Bai et al. (2010) MeV fm and MeV were suggested). Guided by the calculations from RBHF theory with Bonn A interaction, such values for the tensor terms are too strong. How to improve the description of the channel in Pb while not enter in contradiction with ab initio results is an interesting question for future investigations.
Finally, in Fig. 10, we show the SDR strength function of the and channel in Zr calculated by SAMiT with and without tensor, in comparison with experimental data Yako et al. (2006) and SAMi functional RocaMaza et al. (2012). The peak position of has been improved by including the tensor, while for the effect is not significant. The NEWSR calculated by the new functional is fm ( exhausted in the calculation), while the experimental value is fm. The neutron skin of Zr calculated by SAMiT is fm, to be compared with experimental data fm Yako et al. (2006).
V Summary
In this work we have developed a new Skyrme functional SAMiT inspired by the fitting protocol of the successful SAMi functional, with further information on tensor terms provided by ab initio calculations. For that purpose we have calculated for the first time neutronproton drops as predicted by RBHF calculations based on the Bonn A potential. With the evolution of SO splittings as a function of neutron number in the neutronproton drops the tensor terms are well constrained. The new Skyrme functional is then used to investigate the ground and excitedstate properties of nuclei.
For the ground state properties of spherical nuclei, such as binding energy and charge radius, the new functional SAMiT can achieve similar precision as other successful functionals. Therefore, the introduction of strong tensor terms in the Skyrme functional does not decrease the accuracy of describing ground state properties.
By fitting only to the evolution of SO splittings in the neutronproton drops as predicted by RBHF theory, the evolution in the neutron drops by RBHF can be reproduced automatically by the new functional. In other words, the information of the tensor in the neutronneutron channel is consistent in these two systems.
The new functional SAMiT also gives a good description of the excitation energies and sum rules of Giant Resonances. We have investigated the giantmonopole resonance and giantdipole resonance of Pb; the GamowTeller resonance of Ca, Zr, and Pb; the spindipole resonance of Zr and Pb.
The new functional, as in the case of SAMi, respects the empirical hierarchy of the spin and spinisospin LandauMigdal parameters . This has two advantages: i) pushes the model to produce stable results in infinite matter around saturation Pastore et al. (2015). The instability of SAMiT is briefly studied and presented in Appendix B; ii) ensures a better reproduction of spin and spinisospin resonances. For the GTR of Ca and Zr, the lowerenergy and higherenergy peak are both reproduced very well and, to our knowledge, no other interaction of the same type is as accurate as SAMiT in these two cases. The higherenergy peaks of these two cases are well described by both SAMi and SAMiT, but the lowerenergy peaks are improved by SAMiT. For the GTR of Pb, the peak position is also nicely reproduced by SAMiT, similar as SAMi.
The SDR components are also improved comparing with SAMi functional, especially in the case of the channel in Pb which is improved directly by the tensor term. As the NEWSR of SDR is related with the neutron skin, and for Pb and Zr they can be well described by SAMiT, the predictions of neutron skin of these two nuclei by SAMiT are also in good agreement with experimental data.
Acknowledgments
We would like to thank L. G. Cao for cross checking the results; D. Davesne and A. Pastore for discussions and providing results of instability analysis; J. Meng and H. Sagawa for stimulating discussions. This work was partly supported by Funding from the European Union’s Horizon 2020 research and innovation programme under Grant agreement No. 654002.
Appendix A Covariance analysis
In the fit of the new Skyrme functional SAMiT, the is defined as
(23) 
where is the parameters, are observables with data points, ‘theo.’ stands for the calculated values and ‘ref.’ for experimental data or, in this work, pseudodata from RBHF calculations. The corresponding adopted errors are denoted as . See Table. 2 for the detailed information in the fit of SAMiT.
Assuming the reaches a minimum and is a well behaved function of the parameters around the optimal value , it can be expanded and approximated as
(24) 
The curvature matrix is defined as
(25) 
The error (or covariance) matrix is defined as the inverse of the curvature matrix,
(26) 
The error for parameter can be obtained as
(27) 
The correlation between different parameters can be evaluated by the correlation matrix
(28) 
Assuming the observable are smooth functions around the optimal value , one has
(29) 
Within this approximation, one can calculate the covariance between two observables and as
(30) 
The uncertainty of observable can be calculated as
(31) 
The Pearson productmoment correlation coefficient between observables is evaluated as
(32) 
In Fig. 11 we show the Pearson productmoment correlation between different observables given by SAMiT, including: saturation density , saturation energy , effective mass , incompressibility coefficient , centroid energy of GMR in Pb (Fig. 7), symmetry energy , slope parameter of symmetry energy , neutron skin of Pb, centroid energy of IVGDR in Pb (Fig. 7), centroid energy of GTR in Pb (Fig. 8).
Since the first five properties are of isoscalar nature and the second five are of isovector nature, there is a general pattern, as a rule, that the correlations between isoscalar and isovector properties are weaker than those among the same type. For example, the correlations among saturation density, incompressibility, and centroid energy of the GMR are quite strong; the correlations among symmetry energy, slope parameter, neutron skin, and the centroid energy of the IVGDR are also strong.
In principle, the GTR excitation energy should be correlated with the parameter . In the fit of SAMiT, as well as in that of SAMi RocaMaza et al. (2012), this parameter is fixed to a reasonable value so that the GTR can be well reproduced. Since is fixed, this correlation does not exist for SAMiT (or SAMi). From Fig. 11 it can be seen that the GTR centroid energy is correlated with the effective mass. This is understandable as the effective mass will influence the singleparticle level density, which has an important effect on the GTR.
Appendix B Instability analysis
It is known that the Skyrme functional exhibits spurious spin and spinisospin instabilities, which can affect the properties of nuclear groundstate as well as excitedstates Pastore et al. (2015). In Fig. 12, we show the position where instability occurs in the density and momentum transfer plane for (a) SAMiT and (b) SAMi functionals, obtained by solving the linear response function in infinite nuclear matter. The notations here are exactly the same as in Ref. Pastore et al. (2015). Different channels are characterized by quantum numbers of the total spin , projection along the axis, and total isospin . Beside the instability related to the gasliquid phase transition in the channel, the others are unphysical. It can be seen both SAMiT and SAMi show instability near the saturation density, and the tensor term makes the functional more unstable. But as the instability appears at a relatively large value of momentum transfer , it shall not influence much on the discussions of this work.
References
 (1)
 Vautherin and Brink (1972) D. Vautherin and D. M. Brink, Phys. Rev. C 5, 626 (1972).
 Dechargé and Gogny (1980) J. Dechargé and D. Gogny, Phys. Rev. C 21, 1568 (1980).
 Bender and Heenen (2003) M. Bender and P.h. Heenen, Rev. Mod. Phys. 75, 121 (2003).
 Meng (2016) J. Meng, ed., Relativistic Density Functional for Nuclear Structure (World Scientific Pub., 2016).
 Baldo and Burgio (2016) M. Baldo and G. Burgio, Prog. Part. Nucl. Phys. 91, 203 (2016) .
 RocaMaza et al. (2018) X. RocaMaza, G. Colò, and H. Sagawa, Phys. Rev. Lett. 120, 202501 (2018).
 RocaMaza and Paar (2018) X. RocaMaza and N. Paar, Prog. Part. Nucl. Phys. 101, 96 (2018), 10.1016/j.ppnp.2018.04.001.
 Sagawa and Colò (2014) H. Sagawa and G. Colò, Prog. Part. Nucl. Phys. 76, 76 (2014).
 Brown et al. (2010) G. E. Brown, T. T. S. Kuo, J. W. Holt, and S. Lee, eds., The NucleonNucleon Interaction and the Nuclear ManyBody Problem (World Scientific Publishing Co. Pte. Ltd., 2010).
 Lesinski et al. (2007) T. Lesinski, M. Bender, K. Bennaceur, T. Duguet, and J. Meyer, Phys. Rev. C 76, 014312 (2007).
 Schiffer et al. (2004) J. P. Schiffer, S. J. Freeman, J. A. Caggiano, C. Deibel, A. Heinz, C.L. Jiang, R. Lewis, A. Parikh, P. D. Parker, K. E. Rehm, S. Sinha, and J. S. Thomas, Phys. Rev. Lett. 92, 162501 (2004).
 Otsuka et al. (2005) T. Otsuka, T. Suzuki, R. Fujimoto, H. Grawe, and Y. Akaishi, Phys. Rev. Lett. 95, 232502 (2005).
 Otsuka et al. (2006) T. Otsuka, T. Matsuo, and D. Abe, Phys. Rev. Lett. 97, 162501 (2006).
 Brown et al. (2006) B. A. Brown, T. Duguet, T. Otsuka, D. Abe, and T. Suzuki, Phys. Rev. C 74, 061303 (2006) .
 Colò et al. (2007) G. Colò, H. Sagawa, S. Fracasso, and P. Bortignon, Phys. Lett. B 646, 227 (2007).
 Brink and Stancu (2007) D. M. Brink and F. Stancu, Phys. Rev. C 75, 064311 (2007).
 Long et al. (2008) W. Long, H. Sagawa, J. Meng, and N. Van Giai, EPL (Europhysics Lett. 82, 12001 (2008).
 Jiang et al. (2015) L. J. Jiang, S. Yang, B. Y. Sun, W. H. Long, and H. Q. Gu, Phys. Rev. C 91, 034326 (2015).
 Li et al. (2016) J. J. Li, W. H. Long, J. L. Song, and Q. Zhao, Phys. Rev. C 93, 054312 (2016).
 Karakatsanis et al. (2017) K. Karakatsanis, G. A. Lalazissis, P. Ring, and E. Litvinova, Phys. Rev. C 95, 034318 (2017).
 Brink and Stancu (2018) D. M. Brink and F. Stancu, Phys. Rev. C 97, 064304 (2018).
 LópezQuelle et al. (2018) M. LópezQuelle, S. Marcos, R. Niembro, and L. Savushkin, Nucl. Phys. A 971, 149 (2018), 10.1016/j.nuclphysa.2018.01.012.
 Zong and Sun (2018) Y.Y. Zong and B.Y. Sun, Chinese Phys. C 42, 024101 (2018).
 Wang et al. (2018) Z. Wang, Q. Zhao, H. Liang, and W. H. Long, Phys. Rev. C 98, 034313 (2018).
 Colò et al. (1994) G. Colò, N. Van Giai, P. F. Bortignon, and R. A. Broglia, Phys. Rev. C 50, 1496 (1994).
 Niu et al. (2012) Y. F. Niu, G. Colò, M. Brenna, P. F. Bortignon, and J. Meng, Phys. Rev. C 85, 034314 (2012).
 Litvinova and Afanasjev (2011) E. V. Litvinova and A. V. Afanasjev, Phys. Rev. C 84, 014305 (2011).
 Afanasjev and Litvinova (2015) A. V. Afanasjev and E. Litvinova, Phys. Rev. C 92, 044317 (2015).
 Drut et al. (2010) J. E. Drut, R. J. Furnstahl, and L. Platter, Prog. Part. Nucl. Phys. 64, 120 (2010) .
 Day (1967) B. D. Day, Rev. Mod. Phys. 39, 719 (1967).
 Perdew and Kurth (2003) J. P. Perdew and S. Kurth, in A Primer in Density Functional Theory (Springer, Berlin, Heidelberg, 2003) pp. 1–55.
 Dickhoff and Barbieri (2004) W. Dickhoff and C. Barbieri, Prog. Part. Nucl. Phys. 52, 377 (2004).
 Lee (2009) D. Lee, Prog. Part. Nucl. Phys. 63, 117 (2009).
 Liu et al. (2012) L. Liu, T. Otsuka, N. Shimizu, Y. Utsuno, and R. Roth, Phys. Rev. C 86, 014302 (2012) .
 Barrett et al. (2013) B. R. Barrett, P. Navrátil, and J. P. Vary, Prog. Part. Nucl. Phys. 69, 131 (2013).
 Hagen et al. (2014) G. Hagen, T. Papenbrock, M. HjorthJensen, and D. J. Dean, Reports Prog. Phys. 77, 096302 (2014).
 Carlson et al. (2015) J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. E. Schmidt, and R. B. Wiringa, Rev. Mod. Phys. 87, 1067 (2015).
 Hergert et al. (2016) H. Hergert, S. Bogner, T. Morris, A. Schwenk, and K. Tsukiyama, Phys. Rep. 621, 165 (2016).
 Shen et al. (2016) S.H. Shen, J.N. Hu, H.Z. Liang, J. Meng, P. Ring, and S.Q. Zhang, Chinese Phys. Lett. 33, 102103 (2016) .
 Shen et al. (2017) S. Shen, H. Liang, J. Meng, P. Ring, and S. Zhang, Phys. Rev. C 96, 014316 (2017).
 Epelbaum et al. (2009) E. Epelbaum, H.W. Hammer, and U.G. Meißner, Rev. Mod. Phys. 81, 1773 (2009).
 Machleidt and Entem (2011) R. Machleidt and D. Entem, Phys. Rep. 503, 1 (2011).
 Liang et al. (2018) H. Liang, Y. Niu, and T. Hatsuda, Phys. Lett. B 779, 436 (2018).
 Pudliner et al. (1996) B. S. Pudliner, A. Smerzi, J. Carlson, V. R. Pandharipande, S. C. Pieper, and D. G. Ravenhall, Phys. Rev. Lett. 76, 2416 (1996) .
 Gandolfi et al. (2011) S. Gandolfi, J. Carlson, and S. C. Pieper, Phys. Rev. Lett. 106, 012501 (2011) .
 Maris et al. (2013) P. Maris, J. P. Vary, S. Gandolfi, J. Carlson, and S. C. Pieper, Phys. Rev. C 87, 054318 (2013).
 Potter et al. (2014) H. D. Potter, S. Fischer, P. Maris, J. P. Vary, S. Binder, A. Calci, J. Langhammer, and R. Roth, Phys. Lett. B 739, 445 (2014), URL http://www.sciencedirect.com/science/article/pii/S0370269314007436.
 Tews et al. (2016) I. Tews, S. Gandolfi, A. Gezerlis, and A. Schwenk, Phys. Rev. C 93, 024305 (2016), URL https://link.aps.org/doi/10.1103/PhysRevC.93.024305.
 Zhao and Gandolfi (2016) P. W. Zhao and S. Gandolfi, Phys. Rev. C 94, 041302 (2016).
 Shen et al. (2018a) S. Shen, H. Liang, J. Meng, P. Ring, and S. Zhang, Phys. Lett. B 778, 344 (2018a) .
 Shen et al. (2018b) S. Shen, H. Liang, J. Meng, P. Ring, and S. Zhang, Phys. Rev. C 97, 054312 (2018b).
 Bonnard et al. (2018) J. Bonnard, M. Grasso, and D. Lacroix, Phys. Rev. C 98, 034319 (2018).
 Braaten and Hammer (2006) E. Braaten and H.W. Hammer, Physics Reports 428, 259 (2006).
 Buraczynski and Gezerlis (2016) M. Buraczynski, A. Gezerlis, Phys. Rev. Lett. 116, 152501 (2016).
 Buraczynski and Gezerlis (2017) M. Buraczynski, A. Gezerlis, Phys. Rev. C 95, 044309 (2017).
 Machleidt (1989) R. Machleidt, in Adv. Nucl. Phys., Advances in Nuclear Physics, Vol. 19, edited by J. W. Negele and E. Vogt (Springer US, Boston, MA, 1989) pp. 189–376.
 Brown et al. (1987) G. E. Brown, W. Weise, G. Baym, and J. Speth, Comments Nucl. Part. Phys. 17, 39 (1987).
 Sammarruca et al. (2012) F. Sammarruca, B. Chen, and L. Coraggio, and N. Itaco, and R. Machleidt, Phys. Rev. C 86, 054317 (2012) .
 RocaMaza et al. (2012) X. RocaMaza, G. Colò, and H. Sagawa, Phys. Rev. C 86, 031306 (2012) .
 RocaMaza et al. (2013a) X. RocaMaza, G. Colò, and H. Sagawa, Phys. Scr. T154, 014011 (2013a).
 Brueckner et al. (1954) K. A. Brueckner, C. A. Levinson, and H. M. Mahmoud, Phys. Rev. 95, 217 (1954).
 Baranger (1969) M. Baranger, in Nuclear Structure and Nuclear Reactions, Proceedings of the International School of Physics “Enrico Fermi”, Course XL, edited by M. Jean (Academic, New York, 1969), pp. 511–614.
 Davies et al. (1969) K. T. R. Davies, M. Baranger, R. M. Tarbutton, and T. T. S. Kuo, Phys. Rev. 177, 1519 (1969).
 Stancu et al. (1977) F. Stancu, D. M. Brink, and H. Flocard, Phys. Lett. B 68, 108 (1977).
 Chabanat et al. (1998) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 635, 231 (1998).
 Reinhard et al. (1995) P.G. Reinhard, H. Flocard, N. Physics, and G. Reinhard, Nucl. Phys. A 584, 467 (1995).
 Suzuki and Sakai (1999) T. Suzuki and H. Sakai, Phys. Lett. B 455, 25 (1999).
 Wakasa et al. (2005) T. Wakasa, M. Ichimura, and H. Sakai, Phys. Rev. C 72, 067303 (2005).
 Sarsa et al. (2003) A. Sarsa, S. Fantoni, K. E. Schmidt, and F. Pederiva, Phys. Rev. C 68, 024308 (2003) .
 Shen et al. (2018c) S. Shen, H. Liang, J. Meng, P. Ring, and S. Zhang, Phys. Lett. B 781, 227 (2018c).
 Zhou et al. (2003) S.G. Zhou, J. Meng, and P. Ring, Phys. Rev. C 68, 034323 (2003).
 Wang et al. (2017) M. Wang, G. Audi, F. G. Kondev, W. Huang, S. Naimi, and X. Xu, Chinese Phys. C 41, 030003 (2017).
 Angeli and Marinova (2013) I. Angeli and K. Marinova, At. Data Nucl. Data Tables 99, 69 (2013).
 Zalewski et al. (2008) M. Zalewski, J. Dobaczewski, W. Satuła, and T. R. Werner, Phys. Rev. C 77, 024316 (2008).
 Bevington and Robinson (1992) P. R. Bevington and D. K. Robinson, Data Reduction and Error Analysis for Physical Sciences (McGrawHill, New York, 1992).
 Cao et al. (2010) L.G. Cao, G. Colò, and H. Sagawa, Phys. Rev. C 81, 044302 (2010).
 Drischler et al. (2016) C. Drischler, A. Carbone, K. Hebeler, and A. Schwenk, Phys. Rev. C 94, 054307 (2016) .
 Lynn et al. (2016) J. E. Lynn, I. Tews, J. Carlson, S. Gandolfi, A. Gezerlis, K. E. Schmidt, and A. Schwenk, Phys. Rev. Lett. 116, 062501 (2016) .
 Hebeler and Schwenk (2010) K. Hebeler and A. Schwenk, Phys. Rev. C 82, 014314 (2010).
 Tews et al. (2017) I. Tews, J. M. Lattimer, A. Ohnishi, and E. E. Kolomeitsev, Astrophys. J. 848, 105 (2017).
 Dobaczewski et al. (1984) J. Dobaczewski, H. Flocard, and J. Treiner, Nucl. Phys. A 422, 103 (1984).
 Reinhard et al. (1999) P.G. Reinhard, D. J. Dean, W. Nazarewicz, J. Dobaczewski, J. A. Maruhn, and M. R. Strayer, Phys. Rev. C 60, 014316 (1999).
 Goriely et al. (2005) S. Goriely, M. Samyn, J. Pearson, and M. Onsi, Nucl. Phys. A 750, 425 (2005).
 Hosaka et al. (1985) A. Hosaka, K.I. Kubo, and H. Toki, Nucl. Physics, Sect. A 444, 76 (1985).
 Lacombe et al. (1980) M. Lacombe, B. Loiseau, J. Richard, R. Mau, J. Côté, P. Pirès, and R. de Tourreil, Phys. Rev. C 21, 861 (1980).
 Coló et al. (2013) G. Coló, L. Cao, N. Van Giai, and L. Capelli, Comput. Phys. Commun. 184, 142 (2013).
 Youngblood et al. (1999) D. H. Youngblood, H. L. Clark, and Y.W. Lui, Phys. Rev. Lett. 82, 691 (1999).
 Ryezayeva et al. (2002) N. Ryezayeva, T. Hartmann, Y. Kalmykov, H. Lenske, P. von NeumannCosel, V. Y. Ponomarev, A. Richter, A. Shevchenko, S. Volz, and J. Wambach, Phys. Rev. Lett. 89, 272502 (2002).
 Wakasa et al. (1997) T. Wakasa, H. Sakai, H. Okamura, H. Otsu, S. Fujita, S. Ishida, N. Sakamoto, T. Uesaka, Y. Satou, M. B. Greenfield, and K. Hatanaka, Phys. Rev. C 55, 2909 (1997).
 Yako et al. (2009) K. Yako, M. Sasano, K. Miki, H. Sakai, M. Dozono, D. Frekers, M. B. Greenfield, K. Hatanaka, E. Ihara, M. Kato, T. Kawabata, H. Kuboki, Y. Maeda, H. Matsubara, K. Muto, S. Noji, H. Okamura, T. H. Okabe, S. Sakaguchi, Y. Sakemi, Y. Sasamoto, K. Sekiguchi, Y. Shimizu, K. Suda, Y. Tameshige, A. Tamii, T. Uesaka, T. Wakasa, and H. Zheng, Phys. Rev. Lett. 103, 012503 (2009).
 Krasznahorkay et al. (2001) A. Krasznahorkay, H. Akimune, M. Fujiwara, M. N. Harakeh, J. Jänecke, V. A. Rodin, M. H. Urin, and M. Yosoi, Phys. Rev. C 64, 067302 (2001).
 Akimune et al. (1995) H. Akimune, I. Daito, Y. Fujita, M. Fujiwara, M. B. Greenfield, M. N. Harakeh, T. Inomata, J. Jänecke, K. Katori, S. Nakayama, H. Sakai, Y. Sakemi, M. Tanaka, and M. Yosoi, Phys. Rev. C 52, 604 (1995).
 Wakasa et al. (2012) T. Wakasa, M. Okamoto, M. Dozono, K. Hatanaka, M. Ichimura, S. Kuroita, Y. Maeda, H. Miyasako, T. Noro, T. Saito, Y. Sakemi, T. Yabe, and K. Yako, Phys. Rev. C 85, 064606 (2012).
 Bai et al. (2010) C. L. Bai, H. Q. Zhang, H. Sagawa, X. Z. Zhang, G. Colò, and F. R. Xu, Phys. Rev. Lett. 105, 072501 (2010).
 Abrahamyan et al. (2012) S. Abrahamyan, Z. Ahmed, H. Albataineh, K. Aniol, D. S. Armstrong, W. Armstrong, T. Averett, B. Babineau, A. Barbieri, V. Bellini, R. Beminiwattha, J. Benesch, F. Benmokhtar, T. Bielarski, W. Boeglin, A. Camsonne, M. Canan, P. Carter, G. D. Cates, C. Chen, J.P. Chen, O. Hen, F. Cusanno, M. M. Dalton, R. De Leo, K. de Jager, W. Deconinck, P. Decowski, X. Deng, A. Deur, D. Dutta, A. Etile, D. Flay, G. B. Franklin, M. Friend, S. Frullani, E. Fuchey, F. Garibaldi, E. Gasser, R. Gilman, A. Giusa, A. Glamazdin, J. Gomez, J. Grames, C. Gu, O. Hansen, J. Hansknecht, D. W. Higinbotham, R. S. Holmes, T. Holmstrom, C. J. Horowitz, J. Hoskins, J. Huang, C. E. Hyde, F. Itard, C.M. Jen, E. Jensen, G. Jin, S. Johnston, A. Kelleher, K. Kliakhandler, P. M. King, S. Kowalski, K. S. Kumar, J. Leacock, J. Leckey, J. H. Lee, J. J. LeRose, R. Lindgren, N. Liyanage, N. Lubinsky, J. Mammei, F. Mammoliti, D. J. Margaziotis, P. Markowitz, A. McCreary, D. McNulty, L. Mercado, Z.E. Meziani, R. W. Michaels, M. Mihovilovic, N. Muangma, C. Muñoz Camacho, S. Nanda, V. Nelyubin, N. Nuruzzaman, Y. Oh, A. Palmer, D. Parno, K. D. Paschke, S. K. Phillips, B. Poelker, R. Pomatsalyuk, M. Posik, A. J. R. Puckett, B. Quinn, A. Rakhman, P. E. Reimer, S. Riordan, P. Rogan, G. Ron, G. Russo, K. Saenboonruang, A. Saha, B. Sawatzky, A. Shahinyan, R. Silwal, S. Sirca, K. Slifer, P. Solvignon, P. A. Souder, M. L. Sperduto, R. Subedi, R. Suleiman, V. Sulkosky, C. M. Sutera, W. A. Tobias, W. Troth, G. M. Urciuoli, B. Waidyawansa, D. Wang, J. Wexler, R. Wilson, B. Wojtsekhowski, X. Yan, H. Yao, Y. Ye, Z. Ye, V. Yim, L. Zana, X. Zhan, J. Zhang, Y. Zhang, X. Zheng, and P. Zhu (PREX Collaboration), Phys. Rev. Lett. 108, 112502 (2012).
 Kłos et al. (2007) B. Kłos, A. Trzcińska, J. Jastrzȩbski, T. Czosnyka, M. Kisieliński, P. Lubiński, P. Napiorkowski, L. Pieńkowski, F. J. Hartmann, B. Ketzer, P. Ring, R. Schmidt, T. von Egidy, R. Smolańczuk, S. Wycech, K. Gulda, W. Kurcewicz, E. Widmann, and B. A. Brown, Phys. Rev. C 76, 014311 (2007).
 Tarbert et al. (2014) C. M. Tarbert, D. P. Watts, D. I. Glazier, P. Aguar, J. Ahrens, J. R. M. Annand, H. J. Arends, R. Beck, V. Bekrenev, B. Boillat, A. Braghieri, D. Branford, W. J. Briscoe, J. Brudvik, S. Cherepnya, R. Codling, E. J. Downie, K. Foehl, P. Grabmayr, R. Gregor, E. Heid, D. Hornidge, O. Jahn, V. L. Kashevarov, A. Knezevic, R. Kondratiev, M. Korolija, M. Kotulla, D. Krambrich, B. Krusche, M. Lang, V. Lisin, K. Livingston, S. Lugert, I. J. D. MacGregor, D. M. Manley, M. Martinez, J. C. McGeorge, D. Mekterovic, V. Metag, B. M. K. Nefkens, A. Nikolaev, R. Novotny, R. O. Owens, P. Pedroni, A. Polonski, S. N. Prakhov, J. W. Price, G. Rosner, M. Rost, T. Rostomyan, S. Schadmand, S. Schumann, D. Sober, A. Starostin, I. Supek, A. Thomas, M. Unverzagt, T. Walcher, L. Zana, and F. Zehr, Phys. Rev. Lett. 112, 242502 (2014) .
 RocaMaza et al. (2013b) X. RocaMaza, M. Brenna, G. Colò, M. Centelles, X. Viñas, B. K. Agrawal, N. Paar, D. Vretenar, and J. Piekarewicz, Phys. Rev. C 88, 024316 (2013b) .
 Bai et al. (2011) C. L. Bai, H. Q. Zhang, H. Sagawa, X. Z. Zhang, G. Colò, and F. R. Xu, Phys. Rev. C 83, 054316 (2011).
 Yako et al. (2006) K. Yako, H. Sagawa, and H. Sakai, Phys. Rev. C 74, 051303 (2006).
 Pastore et al. (2015) A. Pastore, D. Tarpanov, D. Davesne, and J. Navarro, Phys. Rev. C 92, 024305 (2015).