Finite amplitude method applied to giant dipole resonance in heavy rare-earth nuclei

Background: The quasiparticle random phase approximation (QRPA), within the framework of the nuclear density functional theory (DFT), has been a standard tool to access the collective excitations of the atomic nuclei. Recently, finite amplitude method (FAM) has been developed, in order to perform the QRPA calculations efficiently without any truncation on the two-quasiparticle model space. Purpose: We discuss the nuclear giant dipole resonance (GDR) in heavy rare-earth isotopes, for which the conventional matrix diagonalization of the QRPA is numerically demanding. A role of the Thomas-Reiche-Kuhn (TRK) sum rule enhancement factor, connected to the isovector effective mass, is also investigated. Methods: The electric dipole photoabsorption cross section was calculated within a parallelized FAM-QRPA scheme. We employed the Skyrme energy density functional self-consistently in the DFT calculation for the ground states and FAM-QRPA calculation for the excitations. Results: The mean GDR frequency and width are mostly reproduced with the FAM-QRPA, when compared to experimental data, although some deficiency is observed with isotopes heavier than erbium. A role of the TRK enhancement factor in actual GDR strength is clearly shown: its increment leads to a shift of the GDR strength to higher-energy region, without a significant change in the transition amplitudes. Conclusions: The newly developed FAM-QRPA scheme shows a remarkable efficiency, which enables to perform systematic analysis of GDR for heavy rare-earth nuclei. Theoretical deficiency of photoabsorption cross section could not be improved by only adjusting the TRK enhancement factor, suggesting the necessity of beyond the self-consistent QRPA approach, and/or a more systematic optimization of the EDF parameters.


I. INTRODUCTION
Collective excitations of the atomic nuclei reflects various properties of the nuclear structure and the underlying interaction between nucleons. Their macroscopic or microscopic description has been a major subject in the nuclear theory [1][2][3][4][5][6]. Recently, self-consistent mean-field models, based on the nuclear density functional theory (DFT), have been intensively applied to the collective excitations in heavy open-shell nuclei, where ab initio models are still not computationally feasible.
The giant dipole resonance (GDR) is a noticeable phenomenon generated by the electric dipole excitation. It is basically understood as a collective oscillation of all the neutrons against all the protons, occupying a major part of the nuclear giant resonances [1][2][3][4]. The GDR plays an essential role in the nuclear photoabsorption reaction, determining the centroid energy and width of the cross section. Nuclear photo-absorption reaction impacts also on the dynamics of various astrophysical scenarios [7]. Therefore, GDR can provide a good testing ground for DFT-based theories to describe the nuclear collectivity, as well as the relevant physical properties QRPA is the separable approximation for the residual interaction [15,[24][25][26]. Such a truncation or approximation, however, may invoke the spurious excitations due to the broken self-consistency between the Hartree-Fock-Bogoliubov (HFB) ground state and the QRPA solution.
The finite-amplitude method (FAM) provides an alternative way to solve the QRPA problem with a significantly reduced computational cost. With this method, the QRPA linear response problem is solved iteratively, by circumventing actual calculation and diagonalization of the QRPA matrix. FAM was originally developed for a computation of the RPA strength function, and soon after it was expanded to cover the QRPA problem within spherical symmetry [27,28]. In Ref. [29], FAM-QRPA was implemented into the axially symmetric Skyrme-HFB solver, based on the harmonic oscillator basis. Up to the date, FAM has been implemented also to the axially symmetric coordinate-space HFB solver [30] and to the relativistic mean-field framework [31,32]. Various applications of the FAM include descriptions of giant and pygmy dipole excitations [33,34], efficient computation of the QRPA matrix elements [35], and evaluation of beta-decay rates, including the proton-neutron pairing correlations [36,37]. The contour integration technique of FAM-QRPA was developed to describe individual QRPA modes [38] and for a fast calculation of the energy-weighted sum rules [39]. In addition to FAM, the iterative Arnoldi method presents an alternative method to solve the QRPA problem iteratively [40]. It was also applied to the multipole excitations with pairing correlations [41,42].
This article is devoted to FAM-QRPA methodology applied to the GDR of the heavy rare-earth nuclei, within the Skyrme EDF framework. Due to the open-shell nature of these nuclei, pairing and deformation properties must be taken into account in systematic study. We do not assume any truncation of the two-quasiparticle model space in the QRPA, nor the Bardeen-Cooper-Schrieffer (BCS) approximation for the pairing, but keep a full self-consistency between the HFB and QRPA. To check the validity of the FAM-QRPA, the results are compared with several sets of experimental data. We also investigate the impact of the Thomas-Reiche-Kuhn (TRK) sum rule enhancement factor on the isovector dipole excitation. Because the TRK sum rule is independent on theoretical models and only the enhancement factor (or equivalently, the isovector effective mass) includes the information on the nuclear structure, the energy-weighted sum rule of GDR is an important quantity which reflects the properties of EDFs [5,8,9]. The sensitivity of GDR to the isovector effective mass is also discussed.
We introduce the basic formalism of the Skyrme EDF and FAM-QRPA in the next section. The results are presented and discussed in Sec. III. Finally, we summarize this article in Sec. IV.

II. FORMALISM
As a starting point, our HFB calculations were done in the Skyrme EDF framework. In order to write the Skyrme energy density, it is convenient to introduce the isoscalar and isovector local densities ρ 0 (r) = ρ n (r) + ρ p (r), ρ 1 (r) = ρ n (r) − ρ p (r), (1) where ρ n and ρ p are the neutron and proton densities. With these densities, the Skyrme energy density for the particle-hole (ph) channel reads as where t = 0 (1) indicates the isoscalar (isovector) components. The time-even part E even , is a functional of the local density ρ, kinetic density τ , and spin-orbit density J , whereas the time-odd part E odd is expressed with the spin density s, current density j, and kinetic-spin density T . The detailed formulation of these quantities can be found in, e.g., Refs. [43,44]. The coupling coefficients, C ρρ 0 , etc., are uniquely related to the well-known (t, x) parameterization of the Skyrme force [44,45]. Also, some of coupling constants can be connected to the properties of the symmetric or asymmetric nuclear matter, which are useful pseudo-observables for the optimization purposes of the Skyrme EDF parameters. These pseudoobservables can be treated as alternative EDF input parameters instead of coupling constants [46].
In HFB calculation for the ground state of even-even nuclei, a time-reversal symmetry is usually assumed, and hence, the time-odd part of the functional does not make contribution to the HFB solution. On the other side, when the time-reversal symmetry becomes broken, like in the case of QRPA, the time-odd part becomes active. If we start from the original Skyrme force, the consequent time-odd part of the EDF has a unique correspondence to the time-even part. In other words, when we fix the coupling coefficients in the time-even part, those in the timeodd part should be automatically determined. In the EDF framework, however, a further generalization can be considered: one may treat the time-odd coefficients independently from the time-even ones. In this work, time-odd part is determined as in the case of Skyrme force. For Coulomb energy density, the direct term is treated in a usual manner and for the exchange part we employ Slater approximation. For the particle-particle (pp) channel, which describes nuclear pairing correlations, we employ a functional of the density-dependent delta pairing (DDDP) energy den-sity. That is, whereρ is the local pairing density and ρ c = 0.16 fm −3 is the nuclear saturation density. In this article, a mixed DDDP (ζ = 1/2) is adopted. The pairing strengths V pair q will be adjusted in Sec. III.

A. Finite amplitude method
The detailed formulation of FAM-(Q)RPA can be found in Refs. [27][28][29]38]. We briefly follow these works to arrange the formalism necessary in this work. First, we assume an external time-dependent field, inducing a polarization on the HFB ground state. This external field isF where a † µ and a ν are the quasiparticle creation and annihilation operators, respectively, and η is an infinitesimal real parameter. In this article,F is assumed to be independent of ω, and restricted to have the form of the one-body operator. That is, where c † k and c l are the original particle creation and annihilation operators. The expressions of F 20 µν and F 02 µν in terms of the Bogoliubov transformation can be found e.g. in Refs. [5,28]. Time evolution of quasiparticles is described by the time-dependent HFB equation, where the deviation from the static HFB solution is represented as The quantities needed to obtain the multipole transition strength are the FAM amplitudes, X νµ (ω) and Y νµ (ω), at the excitation energy ω. Since the external field induces density oscillations atop of the static HFB density, the self-consistent Hamiltonian also contains an induced oscillation:Ĥ(t) =Ĥ HFB + δĤ(t), where By solving Eq. (8) up to the first order in η yields socalled FAM equations It is worthwhile to note that, by using the expressions of δH 20 µν (ω) and δH 02 µν (ω) in terms of X µν (ω) and Y µν (ω), one can transform Eq. (11) into the matrix form of where A and B are the well-known QRPA matrices [5]. Notice that Eq. (12) yields the standard matrix form of QRPA when external field is set to zero. The solution of Eq. (12), however, would require to compute the QRPA matrices A and B which generally have large dimension, leading to a substantial CPU time requirement. The essential trick of the FAM-QRPA is to keep Eq. (11), and to solve the FAM amplitudes iteratively with respect to the response of the self-consistent Hamiltonian. This allows to circumvent the large numerical cost of matrix QRPA.
The response of the self-consistent Hamiltonian, δH 20 µν (ω) and δH 02 µν (ω), can be expressed in terms of the induced fields: with the well-known HFB matrices, U and V . Originally the induced FAM-QRPA fields, δh, δ∆ and δ∆, were calculated by applying numerical functional derivatives. In Ref. [47], on the other side, these fields were obtained through explicit linearization of the Hamiltonian, in order not to mix the densities with different magnetic quantum numbers K. Thanks to this explicit linearization, the infinitesimal parameter η is no longer needed, and the induced fields can be formulated in the similar manner as the HFB fields. That is, where h and ∆ are the linearized fields with respect of perturbed densities. These densities can be expressed as The procedures that provide h and ∆ for the HFB solution can be also utilized for the linearized fields, h and ∆ , with a minor modification. For an iterative solution of the FAM amplitudes, the Broyden method was utilized to obtain a rapid convergence [48,49].
By using the FAM-QRPA amplitudes obtained through the iteration, the multipole transition strength distribution is expressed as where i > 0 denotes the summation over the states with positive QRPA energies Ω i > 0, and the response function is given by S(F; ω) = tr[f ρ f ] [28,47]. In order to prevent the FAM-QRPA strength from diverging at ω = Ω i , we employ a small imaginary part in the energy, ω → ω γ = ω + iγ, corresponding to a Lorentzian smearing of Γ = 2γ [28]. The explicit formulation of this smeared strength can be found in Ref. [38]: The contour integration technique is worth to be mentioned: discrete QRPA amplitudes or various multipole sum rules can be obtained from S(F; ω γ ) with a suitable selection of the integration contour on a complex (ω, γ)plane [38,39].
We use following external fields to compute the electric isovector dipole (IVD) strength dB(D K ; ω)/dω, (17) with D K (r) = rY 1K (r). In actual calculation, we replace this operator as Indeed, for an even-even axial nucleus, D K and D −K yields an identical transition strength.

A. Benchmark calculation
The HFB calculations were done by using SkM* Skyrme parameterization at the ph-channel [50]. This set of parameters has been confirmed to be stable in the linear response calculation for the infinite nuclear matter [51]. Since SkM* lacks tensor terms, corresponding time-odd terms were also excluded. For the pp-channel, the pairing strengths for neutrons and protons were adjusted to reproduce empirical pairing gaps of 156 Dy: to 60 MeV. We use computer code hfbtho, which is an HFB solver based on the harmonic oscillator (HO) basis within the axial symmetry [52]. In this work, the HO basis consisted of 20 shells. The imaginary part of ω γ for the FAM strength was set to γ = 0.5 MeV, corresponding to smearing width of Γ = 1.0 MeV, unless otherwise stated.
We would like to emphasize that, in contrast to the standard solution of the QRPA by matrix diagonalization (MQRPA), no truncations on the QRPA quasiparticle model space are imposed in our FAM-QRPA scheme. The only cutoffs employed are the number of HO shells and the pairing window, thus, self-consistency between the HFB and QRPA is fully maintained. Figure 1 shows a benchmark result of FAM-QRPA applied to the GDR. Here we plotted the photoabsorption cross sections obtained with the IVD strengths for 144,154 Sm: The quadrupole matter distribution deformations of the HFB ground states are β = 0 and 0.317 for 144 Sm and 154 Sm, respectively. The IVD strengths of the K = 0 and |K| = 1 modes split in 154 Sm due to the ground state deformation, whereas those are identical for spherical 144 Sm. In both cases, FAM-QRPA shows a good agreement with the experimental data [53]: the typical frequency and width of the GDR can be well reproduced with our model parameters, with the smearing width of Γ = 1.0 MeV. The plateau distribution for the deformed 154 Sm can be understood as a product of the split between K = 0 and |K| = 1 modes [13,14]. Our result also shows a good consistency to that of Ref. [13], in which MQRPA within the coordinate-space representation was adopted.
We have computed FAM strength function within the MPI parallelized scheme, where each part of the strength function was distributed on a separate core, similarly as in Ref. [47]. This scheme achieves a remarkable efficiency, enabling us to compute deformed heavier systems, where MQRPA is available only with a truncation of the model space. Typically, a computation of both of the K-modes took about 1500 CPU hours in a multicore Intel Sandy Bridge 2.6-GHz processor system. Our survey of the GDR has been performed for eveneven rare-earth nuclei from Gd (Z = 64) to W (Z = 74) isochains. Because several sets of experimental data are available [54][55][56], they can provide a more systematic check for FAM-QRPA GDR results. In Tables I and II, we summarize the ground state properties of computed nuclei. The HFB calculation with SkM* concludes that all the nuclei considered here have rather stable prolate deformation.
Our results from FAM-QRPA are summarized in Figs. 2 and 3, in which the photoabsorption cross sections are compared with experimental data, where available. We emphasize here that experimental data in Refs. [56][57][58] does not correspond to total photoabsorption cross section, but the photo-neutron cross sections of (γ, n) and (γ, 2n) type of reactions. Thus, it includes only a part of total photoabsorption cross section due to smaller number of output channels. Generally, we find a reasonable agreement between the FAM-QRPA and the experiments. Typical frequencies of GDR are fairly well reproduced throughout the rare-earth isotopes heavier than Sm. The width and the plateau top of the distribution are well understood as a product from the splitting of K = 0 and |K| = 1 modes, corresponding to the prolate deformation commonly found on their ground states. For 152 Gd and several isotopes of W, the width of the The experimental data sets A, B (photoabsorption), and C (neutron yield) are taken from Refs. [54], [55], and [56], respectively.
GDR is graphically narrower than other nuclei, as expected due to smaller prolate deformation. In our HFB calculations, the proton pairing collapses for 172−178 Yb. This collapse itself, however, does not make a significant impact on the GDR, as the GDR strength distributions look similar irrespective of the proton pairing collapse. Although the pairing could affect the GDR indirectly through the ground state properties (mainly deformation), their changes are small among the rare-earth nuclei, as shown in the present study.
There is an observed deficiency in the calculated photoabsorption cross sections at the region of heavier rareearth isotopes, namely for Z ≥ 70. For example, the calculated photoabsorption cross section of 174 Yb underestimates the experimental data of Ref. [54] at region ω = 12 − 17 MeV in which GDR becomes noticeably strong. A similar kind of GDR deficiency with Skyrme EDFs was reported in Ref. [59].
In Fig. 3, we have also plotted results for the Yb and Hf isotopes by using a smaller smearing width: Γ = 0.5 MeV. With this smaller width, we can reproduce cross section of the neutron yield up to ω ≈ 13 MeV, which includes the first peak of the experimental data (Color online) The same as Fig. 2 but for Yb, Hf, and W isotopes. The experimental data sets A, B(photoabsorption), D(neutron yield), and E(neutron product) are taken from Refs. [54], [55], [57], and [58], respectively. For Yb and Hf, the results obtained with the smearing width of Γ = 0.5 MeV are also plotted with dotted lines.
of Refs. [57,58]. The second peak of the neutron yield cross section may be attributable to the opening of twoneutron emission channel: for 174 Yb, for example, that peak locates at ω = 15 MeV, which is just above the two-neutron separation energy. With Γ = 0.5 MeV, however, we could not achieve a complete improvement of aforementioned deficiency, and the total photoabsorption cross section remains underestimated. In Fig. 4, we also show that the narrower width of Γ = 0.3 MeV leads computed photoabsorption cross sections to overshoot the experimental values at the peak positions, whereas a discrepancy at other frequencies remains. Consequently, the GDR deficiency found here is not improved by simply changing the smearing parameter Γ. Further systematic experiments of photoabsorption cross section, with improved accuracy, would be helpful by providing a more complete testing ground for theoretical models.
Theoretical deficiency of the GDR may be connected to the essential properties of the model. In order to remedy this deficiency, one could consider e.g. beyond QRPA effects or systematic adjustment of the EDF parameters. These are, however, beyond the scope of this article. Alternatively, we discuss a role of the TRK sum rule enhancement factor κ IVD , and its rule in the isovector dipole excitation [8,9]. This quantity can be related to the isovector effective mass of the infinite nuclear matter (INM), which can be used as an input parameter to define the Skyrme EDF parameters [46]. Because there has been some ambiguity about the empirical value of this parameter, knowledge of its effect on GDR will be also profitable for the future optimization of the EDF parameters.

C. Energy-weighted sum rule
To discuss the sensitivity of GDR to the model parameters, we investigate the energy-weighted sum rule (EWSR), defined as In terms of the transition matrix elements, it can be rewritten as It is well known that, by applying the Thouless theorem [60], the EWSR based on the QRPA can be replaced with the expectation value of the double commutator of the HFB ground state [61,62]. For the present case this reads as where κ IVD is the enhancement factor due to the momentum dependence of the effective interaction. For the Skyrme force, it can be given as For the IVD mode, the EWSR has the same value for K = 0 and 1 cases, even if the ground state is deformed. Note also that, for INM, κ NM = 2m(C ρτ 0 − C ρτ 1 )ρ c / 2 is obtained.
Before going to applications, we check the validity of Eq. (22) in a generalized EDF framework [39]. When the EDF is formally generalized, and has no correspondence with respect of the underlying effective force, Thouless theorem is not guaranteed to remain valid. Because we employed the Skyrme EDF combined with the mixed DDDP, the EWSR from actual QRPA calculations can deviate from Eq. (22). In Ref. [39], the authors showed that the Thouless theorem provides a reasonable approximation to the EWSR of the isoscalar/isovector monopole and quadrupole modes, even when Skyrme EDF lacks exact correspondence with respect to the underlying effective interaction but still holds the local gauge invariance. Here we give a similar test for the IVD mode.
For 174 Yb, the energy-weighted sum rule, integrated from the transition strength function up to ω = 50 MeV, yields m 1 (D K ) = 283.9 and 283.8 e 2 fm 2 MeV for K = 0 and K = 1 modes, respectively. Because of the deformation and the resultant splitting of K = 0 and K = 1 strengths, there is a small difference between two values. The contour integration technique of the complex-energy FAM, developed as an efficient tool to compute the sum rules in Ref. [39], yields m 1 (D K=0 ) = 290.1 e 2 fm 2 MeV with an integration contour radius of 200 MeV. On the other side, the double commutator procedure of Eq. (22) gives m 1 (D K=0,±1 ) = 287.7 e 2 fm 2 MeV, with enhancement factor of κ IVD = 0.397, when using HFB proton and neutron densities.
We find the value from the double commutator method being consistent to those from the FAM-QRPA calculation. Consequently, Thouless theorem can provide a reliable approximation of the IVD sum rule. The computational cost for the double commutator procedure is drastically lighter compared to the FAM-QRPA, since it requires information only about the HFB state.
We have summarized the EWSR values and enhancement factors of Eqs. (22) and (23) in Tables I and II. It is clearly shown that the enhancement factor for the TRK sum rule is almost constant at this region of the nuclear chart. The ratio of two enhancement factors, is approximately 0.74 for rare-earth systems calculated here with SkM* Skyrme parameterization. This is simply due to similar value obtained from the density integration of Eq. (24), with only a limited variation on the proton and neutron density profiles with respect of N and Z.
In Fig. 5, we plot the enhancement factor κ IVD and the total HFB energy E HFB as a function of the axial deformation parameter β for 174 Yb. These are obtained from the HFB calculation with a constraint on β. Approximately, the minimum of the total HFB energy corresponds to the maximum of κ IVD . This can be understood mainly in terms of the symmetry energy, which favors a large overlap between neutrons and protons. Due to other ingredients, especially the Coulomb energy, the true ground state and maximum of κ IVD do not exactly coincide. This curve shows, however, that the IVD enhancement factor is not very sensitive to the details of the ground state deformation. Next, we investigate if the deficiency against the experimental photoabsorption data found in 174 Yb could be improved by changing κ NM parameter.

D. Sensitivity to enhancement factor
In the remaining part of this section, we investigate the sensitivity of GDR to the EDFs with different κ NM values. The Skyrme parameterizations suited to this purpose can be found in Ref. [44], where the authors optimized Skyrme parameters by assuming systematic constraints on various INM properties. The unconstrained parameterization, SV-min, was optimized without constraints, and SV-bas parameterization was a base starting point for INM parameter variation. The parameterizations with a variation on κ NM were introduced as SV-kap60 and SV-kap20. Because these constrained parameterizations were otherwise optimized exactly in the same manner as SV-bas, we can check the effect of κ NM on the GDR in a systematic manner. The exact value of κ NM as an input parameter for each functional is present in Table III. Note that the larger value of the enhancement factor corresponds to the lighter isovector effective mass (κ NM = m/m * v − 1) [9,44]. Except for the interaction parameterization employed, the numerical conditions are the same as in the previous calculations.
In Table III, the EWSR of 174 Yb obtained with SkM*, SV-kap20, SV-bas, and SV-kap60 are summarized with the corresponding enhancement factors. As naturally expected from the definitions of κ NM and κ IVD [9,39,44], the EWSR increases with the enhancement factor of INM, in other words, with the reduction of the isovector effective mass. The HFB ground states computed with three SV-functionals are similarly deformed. Thus, the density integration of Eq. (24) is also similar for all three functionals, yielding the common ratio of κ IVD /κ NM ∼ = 0.74.
In Fig. 6, IVD transition strengths obtained with SV-kap20, SV-bas, and SV-kap60 are plotted. One can clearly find that an increment of the enhancement factor leads to a shift of the IVD strength towards higher energy region, as pointed out in Ref. [44] for the GDR of doublymagic 208 Pb. Conversely, the shifted strength trivially yields an enhanced value due to its energy-weight in the EWSR. That is, if the isovector effective mass becomes lighter, the corresponding excitation energies becomes higher. This can be qualitatively explained within an analogy to the single-particle energies of the HO potential, which have the same curvature, but different particle masses. From our calculations, the form of the ground state density is found to be insensitive to κ NM . Thus, in this kind of case, the excitation energies are mainly determined by the effective mass for the collective motion. The energy interval between the K = 0 and K = 1 peaks is not sensitive to κ NM , due to similar deformation parameters listed in Table III. It is also noticeable that the total amplitude of the transition strength function is not significantly changed for different parameterizations, in both K = 0 and K = 1 cases. Thus, only the position of the peak is sensitive on κ NM . From these results, we expect that the mean frequency of GDR to be a suitable observable in order to constrain isovector effective mass parameter during EDF optimization.
The photoabsorption cross sections for SV-EDFs are plotted in Fig. 7. The shift of GDR structure to the high energy region is again observed by the increase of κ NM . The SV-kap20 functional reproduces the experimental data well up to ω = 11 MeV, whereas the general structures matches best to SV-bas result, which corresponds to κ NM = 0.4.
On the other side, however, the GDR deficiency against the experimental data still remains: calculated total photo-absorption cross section cannot be improved by changing the enhancement factor. In order to improve current EDF models, a parameter optimization, combined with an input data on the GDR position and magnitude, may help the situation. Yb calculated with SV-kap20, SV-bas, and SV-kap60 parameterizations. The shown experimental data is the same as in Fig. 3.

IV. SUMMARY
We have performed a systematic analysis of the GDR in heavy rare-earth elements, including neutron-rich and deformed isotopes. The calculations were handled in the recently developed parallel FAM-QRPA scheme, in combination with the Skyrme EDF, without any additional truncations on the two-quasiparticle model space. This scheme enabled to perform fully self-consistent QRPA calculations efficiently and free from the spurious effects due to the broken self-consistency.
The mean energy and width, as well as the plateau shape of the photoabsorption cross section, have been fairly well reproduced for nuclei considered. However, some deficiency in the calculated total photoabsorption cross section was seen for Z ≥ 70 isotopes.
We also investigated the behavior of GDR by changing the TRK enhancement factor, connected to the isovector effective mass. It is clearly shown that an increment of the enhancement factor shifts the GDR distribution towards the higher-energy region, corresponding to the lighter isovector effective mass.
The deficiency of GDR total photoabsorption cross section, with respect to the experimental data, remained to be an open question. This deficiency is noticeable at the region of isotopes heavier than Er. This situation could not be improved by tuning smearing width or κ NM . Several possibilities of further improvements can be proposed. The first is to expand the framework to cover the dynamics beyond the QRPA, as well as the other multipole degrees of freedom. Especially, the octupole softness of systems could play a role at higher than the RPA level. Another direction of progress is to perform a more systematic optimization of EDF parameters, and to use GDR data on deformed nuclei as an input. Espe-cially, the tangled effect of the symmetry energy and its slope with the isovector effective mass is expected to be important [9,15,63].
Our FAM-QRPA scheme could be also employed to investigate the low-lying excitations or the pygmy strength. For the excitation energies as well as the partial sum rules of these resonances, not only the particle-hole part but also the pairing part of the EDF is expected contribute notably. Especially, for the low-lying excitations of loosely bound nuclei, the HFB method has an advantage over the BCS method for the treatment of the pairing correlations, especially for the nuclei close to the neutron-drip line [64,65]. For this purpose, FAM-QRPA embedded into the coordinate-space HFB solver could be a good choice of the method [30].