Nuclear dipole polarizability from mean-field modeling constrained by chiral effective field theory

We construct a new Skyrme interaction Sk$\chi$m$^*$ by fitting the equation of state and nucleon effective masses in asymmetric nuclear matter from chiral two- and three-body forces as well as the binding energies of finite nuclei. Employing this interaction to study the electric dipole polarizabilities of $^{48}$Ca, $^{68}$Ni, $^{120}$Sn, and $^{208}$Pb in the random-phase approximation, we find that the theoretical predictions are in good agreement with experimentally measured values without additional fine tuning of the Skyrme interaction, thus confirming the usefulness of the new Skyrme interaction in studying the properties of nuclei. We further use this interaction to study the neutron skin thicknesses of $^{48}$Ca and $^{208}$Pb, and they are found to be consistent with the experimental data.


I. INTRODUCTION
The nuclear electric dipole polarizability α D , which is proportional to the inverse energy-weighted sum of the dipole response of a nucleus to an external electric field, has recently attracted much attention due to its strong correlation with neutron skin thicknesses, the density dependence of the nuclear symmetry energy [1][2][3][4][5][6][7][8][9], and the properties of neutron stars [10][11][12]. The electric dipole response is dominated by the nuclear isovector giant dipole resonance (IVGDR), which is the oldest known nuclear collective excitation and has been extensively studied both theoretically and experimentally [13][14][15][16]. It is well known from theoretical studies based on various models that the properties of the IVGDR, in which neutrons oscillate against protons in a nucleus, is affected by the density dependence of the nuclear symmetry energy [17][18][19][20][21].
Experimentally, the electric dipole strength distributions in 48 Ca [22], 120 Sn [23], and 208 Pb [24] have recently been measured accurately at the Research Center for Nuclear Physics (RCNP) from proton inelastic scattering experiments at forward angles, while that of 68 Ni has been investigated at GSI using Coulomb excitations in inverse kinematics and measuring the invariant mass in the one-and two-neutron decay channels [25]. Recently, ab initio calculations based on chiral effective field theory interactions have been successful at describing α D in medium-mass nuclei [26,27], but the description of the dipole response of heavy nuclei remains a challenge. Therefore, nuclear energy density functionals are still the most widely-used approach to explore the nuclear equation of state (EOS) from the electric dipole response.
In the present work, we construct a new Skyrme interaction Skχm * by fitting the equation of state and nucleon effective masses in asymmetric nuclear matter predicted by chiral two-and three-body forces as well as the binding energies of finite nuclei. We then employ it to study the dipole response of 48 Ca, 68 Ni, 120 Sn, and 208 Pb within the random-phase approximation (RPA). The calculated electric dipole polarizabilities are in good agreement with existing experimental values. In addition, the neutron skin thickness in 48 Ca and 208 Pb is discussed.
The paper is organized as follows. In Section II, we explain in detail the new Skyrme interaction Skχm * and how it is obtained. In Section III, we describe the random-phase approximation employed in our calculations, and in Section IV we present results and discussions for the isovector dipole response, electric dipole polarizability, and neutron skin thicknesses of 48 Ca and 208 Pb from the new Skyrme mean-field model. We end with a summary in Section V.

II. A SKYRME INTERACTION FROM CHIRAL EFFECTIVE FIELD THEORY
Recently, extended Skyrme interactions have been constructed in Ref. [28] by fitting to the asymmetric nuclear matter equation of state from chiral two-and threenucleon forces [29,30]. In contrast to previous works [31,32] imposing constraints on mean-field models from low-density neutron matter, the full density dependence of the asymmetric matter equation of state was used in Ref. [28] for the χ 2 minimization function. In these calculations, the chiral two-body force was treated at nextto-next-to-next-to-leading order (N3LO) in the chiral power counting, while the three-body force was treated at N2LO. Three different choices of the momentum-space cutoff were considered {Λ = 414, 450 and 500 MeV}, and the unknown low-energy constants associated with short-distance dynamics were fitted [33][34][35][36] in each case to nucleon-nucleon scattering phase shifts and deuteron properties (in the case of the nucleon-nucleon interaction) as well as the binding energy and beta-decay lifetime of 3 H (in the case of the three-nucleon force). To fix the gradient contributions to the energy density functional, the extended Skyrme interactions were also fitted to reproduce the binding energies of 7 doubly-closed-shell nuclei: 16 O, 40 Ca, 48 Ca, 56 Ni, 100 Sn, 132 Sn, and 208 Pb. The conventional Skyrme interaction [37] was found to be inadequate to describe the density dependence of the equation of state from chiral effective field theory, leading the authors of Ref. [28] to include an additional densitydependent and momentum-independent term. The three Skyrme interactions (Skχ414, Skχ450 and Skχ500) reported in Ref. [28] predict, however, too large (even larger than the nucleon bare mass m) isoscalar and isovector effective masses m * s,0 and m * v,0 at saturation density, which are given by [37] 1 m * This fact makes Skχ414, Skχ450 and Skχ500 unsuitable for studying the nuclear giant dipole resonances. It is well known that the m * v,0 is related to the enhancement factor κ in the Thomas-Reiche-Kuhn sum rule [38], given by the relation 1 + κ = m/m * v,0 , and thus affects the nuclear dipole response function. We note that in Skyrme functionals the neutron and proton effective masses m * n and m * p at total baryon number density ρ are related to the isoscalar and isovector effective masses m * s (ρ) and m * v (ρ) by [39] 1 where δ = (ρ n − ρ p )/(ρ n + ρ p ) is the isospin asymmetry with ρ n (ρ p ) being the neutron (proton) density.
In the present work, instead of directly fitting the three EOSs predicted by chiral effective field theory as in Ref [28], we treat them as theoretical uncertainties on the equation of state and construct a Skyrme interaction by fitting the central values. As shown later, taking into account the theoretical uncertainties makes it possible to use the conventional Skyrme interaction to reproduce the predicted EOS from chiral effective field theory. To improve the description of nucleon effective masses, we further include in our fit the isoscalar and isovector effective masses m * s,0 /m = 0.82 ± 0.08 and m * v,0 /m = 0.69 ± 0.02 at saturation density that are extracted from the single-particle energies of protons and neutrons in asymmetric nuclear matter, computed recently in Refs. [40,41] from chiral two-and three-body forces. From the momentum-and energy-dependent nucleon single-particle energy e(p) = p 2 2m + ReΣ(p, e(p)), where p is the nucleon momentum and the self-energy Σ is computed at second order in perturbation theory, the nucleon effective mass m * can be extracted by the definition p m * = de dp .
Since the effective masses in the Skyrme-Hartree-Fock model are momentumindependent, we average the momentum-dependent effective masses from chrial effective theory over a momentum range around the Fermi momentum, which accounts for most of the theoretical uncertainties. It is interesting to see that the m * v,0 from chiral effective field theory is consistent with that from an analysis of dipole resonances in 208 Pb based on conventional Skyrme interactions [42].
Our new Skyrme interaction therefore has the usual form: where σ i is the Pauli spin operator, P σ = (1 + σ 1 · σ 2 )/2 is the spin-exchange operator, k = −i(∇ 1 − ∇ 2 )/2 is the relative momentum operator, and k is the conjugate operator of k acting on the left. Following Refs. [43,44], we express the 9 parameters t 0 −t 3 , x 0 −x 3 and α of the Skyrme interaction in terms of 9 macroscopic quantities: m * s,0 , m * v,0 , the nuclear matter saturation density ρ 0 , the energy per particle of symmetric nuclear matter E 0 (ρ 0 ), the incompressibility K 0 , the gradient coefficient G S , the symmetry-gradient coefficient G V , and the magnitude E sym (ρ 0 ) and density slope L of the nuclear symmetry energy at ρ 0 [43,44]. Here the G S and G V are defined by expressing the momentumdependent or finite-range terms of the Skyrme energy density functional H fin as [43] with Consequently, the Skyrme mean-field models have the following 10 parameters: These parameters have clear physical meanings and their ranges are empirically known, which increases significantly the efficiency of the optimization algorithm. The parameters G S , G V and W 0 , which are not constrained by the nuclear equation of state, have been studied [45][46][47] consistently in chiral effective field theory via the density matrix expansion in the Hartree-Fock approximation. However, higher-order terms in the expansion are needed to produce reliable constraints on the associated parameters in Skyrme functionals. In the present work we constrain the values of G S , G V and W 0 by fitting to the binding energies of double-closed-shell nuclei as in Ref. [28]. As in the usual fitting procedure, we optimize the parameters p(p 1 , p 2 , · · · , p N ) by minimizing the weighted sum of the squared deviations from the data as where E i is the binding energy per nucleon in nuclear matter with the nucleon density ρ i = {0.01, 0.04, 0.08, 0.16, 0.20, 0.24, 0.28, 0.32} fm −3 and the isospin asymmetry δ 2 i = {0, 0.25, 0.5, 0.75, 1}; B i is the total binding energy of the ith nucleus; 'Sky' and 'EFT' indicate predictions of the Skyrme-Hartree-Fock model and chiral effective field theory, while 'exp' indicates the experimental data. The adopted error σ B is set to be 1 MeV, while σ E is taken to be max(0.25, σ EFT (ρ)) MeV where σ EFT (ρ) is the theoretical error of the EOS due to the variation of the momentum cutoff used in chiral effective field theory. Since the predictions of chiral effective field theory are very precise at low densities (e.g., the error is 0.02 MeV in neutron matter at ρ = 0.01 fm −3 ), a minimum error of 0.25 MeV for the EOS is introduced to take account of the deficiency of the Skyrme parametrization. The optimization is carried out by employing the POUNDERS algorithm in the distribution of PETSc/Tao [48], which has been successfully used previously to calibrate nuclear energy density functionals [44,49,50].
Once the optimized parameter set p 0 is obtained, the χ 2 near p 0 can be approximately written as with M ij = ∂ 2 χ 2 ∂pi∂pj p0 . The errors of parameters then can be estimated as which is the allowed variation of the parameter p i within the ellipsoid of χ 2 −χ 2 (p 0 ) ≤ 1. Using the covariance matrix Σ = M −1 , the correlation matrix can be estimated as where C ij is the correlation coefficient between parameters p i and p j . The extrapolation error of any observables and the correlation coefficient between them can also be determined (see e.g., [51,52]). The obtained parameter set p 0 , named Skχm * , including its uncertainties, together with the corresponding Skyrme parameters are listed in Tab. I. It is seen that the parameters ρ 0 , E 0 (ρ 0 ), K 0 , E sym (ρ 0 ) and L, which characterize the density-and isospin-dependence of the nuclear EOS, are well constrained by the predicted EOS from chiral effective field theory. If we had removed from the fit all data points shown in Fig. 1 for the equation of state of asymmetric matter, the uncertainty in L would be roughly 60 MeV. Also, the isoscalar gradient coefficient G S = 141.5 ± 7.3 MeV · fm 5 and the spin-orbit coupling strength W 0 = 119.8±4.9(MeV·fm 5 ) are consistent with the values at sub-saturation densities of about 0.03 fm −3 in the energy density functional derived from chiral two-and three-body forces [47]. We show in Fig. 1 the EOSs of asymmetric nuclear matter predicted by the Skχm * interaction together with predictions from chiral effective field theory. One can see that the new interaction can well reproduce the EOSs from chiral effective field theory calculations. In Tab. II, we list the binding energies and charge radii of 16 O, 40 Ca, 48 Ca, 56 Ni, 100 Sn, 132 Sn, and 208 Pb obtained by Hartree-Fock calculations using Skχm * , together with their experimental values. It is seen that our results are overall in very good agreement with the experimental data for the 7 doubly-closed-shell nuclei.   Fig. 2 shows the absolute values of correlation coefficients calculated among the 10 parameters ρ 0 , E 0 (ρ 0 ), K 0 , E sym (ρ 0 ), L, G S , G V , W 0 , m * s,0 and m * v,0 . We find that m * v,0 is almost independent of all other parameters, which reflects the fact that it is poorly constrained by the data from the EOSs and nuclear binding energies used in our fit. Except m * v,0 , most parameters are correlated with each other, but generally the correlations are weak with the absolute value of the Pearson coefficient R less than 0.8. The spin-orbit coupling strength W 0 is strongly correlated with both the isoscalar effective mass m s,0 and the isoscalar gradient coefficient G S , which are related to the momentum-dependent terms in the Skyrme interaction.

III. ISOVECTOR DIPOLE RESPONSE AND RANDOM-PHASE APPROXIMATION
The dipole response of a nucleus can be calculated in the random-phase approximation using the dipole oper-atorF where r i is the nucleon's radial coordinate and Y 1M (r i ) is the corresponding spherical harmonic function. In the RPA method, the isovector dipole strength function is evaluated as where |0 is the RPA ground state, and |ν is the RPA excited state with E ν being its energy. Defining the moments of the strength function as the electric dipole polarizability α D can then be calculated according to In the present work, we employ the Skyrme-RPA code by Colo et al. [53] to calculate the electric dipole polarizabilities of 48 Ca, 68 Ni, 120 Sn, and 208 Pb. As pointed out  (14) 208 Pb 19.6(6) 19.87 (29) in Ref. [5], to compare the RPA results with experiment, the experimentally extracted values of α D for 120 Sn [23] and 208 Pb [24] need to be corrected by subtracting the contributions from quasi-deuteron excitations [54,55], while that of 68 Ni [26] should be modified by including the corrections from the extrapolated low-energy and high-energy regions. The corrected experimental data are shown in Table III.  48 Ca, 68 Ni, 120 Sn, and 208 Pb obtained from RPA calculations using Skχm * . For comparison, experimental data [22][23][24][25] are shown as black solid circles.
In Fig. 3, we show the isovector dipole transition strength functions for 48 Ca, 68 Ni, 120 Sn, and 208 Pb obtained from the RPA calculations using the Skyrme interaction Skχm * . Here the curves are obtained by smearing each of the discrete RPA peaks by a Lorentzian function of 2 MeV width, which is taken to reproduce the experimental width of the IVGDR. For comparison, the measured strength functions for these nuclei are shown as black solid circles, and the experimental peak energies E x = 18.9, 17.1, 15 and 13.43 MeV [22][23][24][25]56] of the IVGDR in 48 Ca, 68 Ni, 120 Sn, and 208 Pb are indicated by the arrows. It is seen that the Skχm * interaction gives dipole strength functions which are consistent with experimental measurements, except the predicted second peak at higher energy that is not found in experimental data, especially for 120 Sn. Such an unphysical highenergy peak is expected to become a long high-energy tail if more complex configurations beyond the 1p − 1h states are included in the RPA calculation [57]. For 208 Pb, the predicted peak energy E x from Skχm * at 13.6 MeV is in particularly good agreement with the experimental value. This is due to the m * v,0 = 0.694m, or equivalently κ = 0.44, predicted by chiral effective field theory, which is roughly consistent with the value κ = 0.4 determined by fitting the IVGDR peak energy in 208 Pb using conventional Skyrme interactions [42].
The relation between m * v,0 or κ and E x can be understood from the approximate relation [18,21] While the inverse energy weighted sum rule m −1 [Eq. (17)] is related to the electric dipole polarizability via Eq. (18), the energy weighted sum rule m 1 depends inversely on the isovector effective mass m * v,0 at saturation density in the conventional Skyrme mean-field model [58]. The latter can be understood by the fact that the m 1 can be expressed as [53] where the factor κ is inversely proportional to m * v,0 . For 68 Ni, panel (b) of Fig. 3 shows that there is a socalled pygmy dipole resonance around 9.55 MeV in the experimental strength function that carries 2.8(5)% of the energy weighted sum rule m 1 [25]. Due to the smearing of the RPA states in the calculation, the pygmy resonance is not seen and appears instead as a low-energy tail in the theoretical strength function. The RPA strength in the energy region of 0-11 MeV carries 2.7% of the total energy weighted sum rule and is similar to that carried by the pygmy dipole resonance. Its contribution to the dipole polarizability is 0.38 fm 3 and thus underestimates the experimental value of 1.13 fm 3 [5,25] as for most Skyrme-type interactions. Fig. 4 exhibits the electric dipole polarizabilities of 48 Ca, 68 Ni, 120 Sn, and 208 Pb predicted by the Skχm * mean-field model together with the experimental results [22][23][24][25] (detailed numerical values are listed in Table III). We find that RPA calculations using the Skχm * interaction can overall well reproduce the experimental data. For 120 Sn the prediction of Skχm * is slightly larger than the experimental value. This could be due to the fact that, in the present work, we do not include pairing correlations, which reduce the electric dipole polarizability in the Sn isotope chain [5]. In the case of 48 Ca, we observe that the results are consistent with recent ab initio calculations [26]  effective field theory. We would like to point out that we have calculated the electric dipole polarizability by using discrete RPA peaks because the 1p − 1h RPA cannot give rise to the spreading width of nuclear giant resonances observed in experiments. Including the spreading width, which can be done by taking into account the coupling of 1p − 1h states with more complicated multiparticle-multi-hole configurations, is expected to reduce the calculated α D slightly. This can be seen by smearing the RPA peaks using the Lorentzian function with a certain width. In the case of only one Lorentzian function of width Γ, the electric dipole polarizability is at most reduced by [5] ∆α According to this formula, using the experimental value of the peak energy (13.43 MeV) and width (4.07 MeV) of the IVGDR [56], we estimate that the reduction should be less than about 0.45 fm 3 (∼ 2%) for 208 Pb. Such a small correction does not influence our conclusion. The neutron skin thickness, which has been demonstrated to be strongly correlated with the electric dipole polarizability, is an important probe of the density slope of the nuclear symmetry energy as well as the EOS of pure neutron matter. Therefore, we also employ the Skχm * mean-field model to predict the neutron skin thicknesses of 48 Ca and 208 Pb, and the obtained results are given by ∆r np ( 48 Ca) : 0.167 ± 0.004 fm, ∆r np ( 208 Pb) : 0.170 ± 0.005 fm.

V. SUMMARY
We have constructed a new Skyrme interaction Skχm * by fitting the EOSs and nucleon effective masses in asymmetric nuclear matter predicted by chiral effective field theory together with the binding energies of selected closed shell doubly-magic nuclei. The new model was employed to study the isovector dipole response of 48 Ca, 68 Ni, 120 Sn and 208 Pb. We have found that the new interaction can well reproduce the experimental data on the peak energy of the giant dipole resonance and the electric dipole polarizability. We have further calculated the neutron skin thicknesses of 48 Ca and 208 Pb. The predicted neutron skin thickness of 208 Pb from our study is also consistent with the experimental values [59][60][61]. Our results thus confirm the usefulness of the new Skyrme interaction Skχm * in studying the isovector properties of nuclei in regimes where ab initio calculations with chiral nuclear forces may not be feasible.