Numerical Convergence of Electromagnetic Responses with the Finite-Amplitude Method

The response of a nucleus to an electromagnetic probe is a key quantity to simulate photabsorption or photodeexcitation processes. For large calculations at the scale of the entire mass table, this response can be estimated by linear response theory. Thanks to the introduction of the finite-amplitude method (FAM), calculations are computationally efficient. In this paper, we investigate in more details the convergence of FAM calculations of the response function as a function of the parameters controlling the numerical implementation of the theory. We show that the response is much less sensitive to the details of the single-particle basis than, e.g., Hartree-Fock-Bogoliubov calculations.


Introduction
Photonuclear processes play central roles as probes to nuclear structure and in applications of nuclear physics to nuclear technology, isotope production or astrophysics.Much work is thus devoted to building reliable and comprehensive data libraries of properties such as γstrength functions [1].Notwithstanding fully phenomenological approaches based on explicit fit to experimental data, large-scale theoretical predictions of such quantities typically rely on the quasiparticle random-phase approximation (QRPA) and linear response theory [2][3][4][5][6].Because of their unfavorable computational scaling, most of these calculations have often been restricted to spherical symmetry and/or even-even nuclei.
With the introduction of the finite-amplitude method (FAM), the computational cost of the linear response in deformed nuclei became considerably lower and is now somewhat comparable to the one of a deformed Hartree-Fock-Bogoliubov (HFB) calculation [7,8].This has enabled studies of nuclear multipole responses in heavy, deformed nuclei [9][10][11][12] as well as more accurate estimates of nuclear moments of inertia and the collective inertia tensor [13,14].The FAM was originally introduced and tested with non-relativistic Skyrme functionals and later also implemented with relativistic functionals [15].
In most of these early applications, the uncertainties of the results induced by numerical inputs such as the basis parameters, quasiparticle cutoff and smearing width were not fully discussed, with the exception of Refs.[10] and [12] which presented results for two different basis sizes.Such uncertainties should be quantified before embarking on any large-scale calculations of the response function.The goal of this paper is thus to perform a more systematical study of the convergence of FAM calculations with Skyrme functionals as a function of LLNL-PROC-855567 various numerical features of the implementation.We briefly recall the theoretical formalism in Section 2 before presenting a selection of the results in Section 3.

Theoretical Framework
We use the general framework of the energy density functional (EDF) to compute the electromagnetic response of nuclei.The theory is presented in great details in textbooks and review articles, most notably [16][17][18][19][20].In the following, we only recall the most essential formulas.

Hartree-Fock-Bogoliubov Theory with Skyrme Functionals
The ground-state many-body wavefunction |Φ⟩ of the nucleus is computed at the Hartree-Fock-Bogoliubov (HFB) approximation.We recall that |Φ⟩ is then a product state of quasiparticle annihilation operators that are related to any single-particle basis by a (unitary) Bogoliubov, or canonical, transformation W characterized by the matrices U and V, where β, β † are the quasiparticle operators and c, c † the single-particle operators.The one-body density matrix ρ = V * V T and pairing tensor κ = V * U T can both be expressed as functions of the Bogoliubov transformation, and the total energy is a functional of ρ, κ and κ * .In practice, the matrices U and V are determined by requiring that the total energy E[ρ, κ, κ * ] be a minimum with respect to all independent matrix elements ρ i j , κ i j , and κ * i j .The solution to the resulting HFB equation fully determines the Bogoliubov matrices U and V.
In this paper, we choose the functional form E[ρ, κ, κ * ] to be a Skyrme functional in the particle-hole channel (terms proportional to ρ) and a simple zero-range, density-dependent functional in the particle-particle channel (terms proportional to κ and κ * ).As is well known, the zero range of the pairing functional requires applying a cutoff on the number of quasiparticles used to compute the density matrix and pairing tensors.

Linear Response Theory with the Finite-Amplitude Method
Linear response theory quantifies the effect of applying a small one-body perturbation operator F(t) to the ground state described by the many-body wavefunction |Φ⟩ and is most naturally derived as the small-amplitude limit of time-dependent Hartree-Fock-Bogoliubov theory (TDHFB) where R(t) is the time-dependent generalized density and Ĥ(t) the time-dependent HFB matrix.Here, one will assume that where the frequency ω plays the role of the excitation energy in linear response theory.Under the effect of the perturbation operator, the generalized density R0 associated with the HFB ground state |Φ⟩ becomes, to the first order, R(t) = R0 + δ R(ω)e −iωt + δ R(ω) † e iωt , and correspondingly Ĥ(t) = Ĥ0 + δ Ĥ(ω)e −iωt + δ Ĥ(ω) † e iωt .In the Bogoliubov (quasiparticle) basis, LLNL-PROC-855567 the perturbed density δ R(ω), perturbed mean field δ Ĥ, and perturbation operator F(ω) take the form The standard, matrix equation of linear response reads where the matrices A and B together form the QRPA matrix and are obtained from the second derivatives of the energy with respect to the generalized density R. Solving ( 6) directly can become computationally prohibitive in the case of heavy deformed nuclei where the number of relevant quasiparticles can be of the order of 10 3 , leading A and B to be complex-valued dense matrices of size of about 10 6 .

Finite-Amplitude Method
The finite-amplitude method (FAM) provides an alternative way to determine the amplitudes X and Y characterizing the perturbed density and scales much more favorably than the direct matrix approach [7,8].The FAM yields a system of coupled equations where δH 20 and δH 02 are functionals of the perturbed densities given by the Bogoliubov transformation of δR(ω) as Given the initial HFB solution characterized by the matrices U and V and the matrix elements F 20 and F 02 of the perturbation operator F, equations ( 7)-( 8) can be solved iteratively.
After the X and Y amplitudes have been determined, the actual excitation strength can be readily obtained as with the response function S (ω) given by In practice, calculations are performed at complex values of the frequency ω as ω → ω + i Γ 2 , where Γ is the smearing width.This is not only mathematically motivated by the need to avoid the poles of the response function on the real axis, but also corresponds, physically, to folding the solutions to (6) with a Lorentzian distribution.In the following we will choose Γ = 1.0 MeV if it is not explicitly specified. LLNL-PROC-855567

Quasiparticle Cutoff
Additional complications emerge when performing practical calculations with zero-range pairing forces because of the ultra-violet divergence of the pairing tensor [20][21][22].At the HFB level, all quasiparticles µ such that E µ > E cut are discarded.Depending on the value of the actual cutoff energy and of the parameters of the HO basis, this implies that the U and V matrices may not be square but of dimension (N, N qp ) with N the size of the basis and N qp ≤ N the actual number of active quasiparticles.In such a case, the one-body density matrix reads How to handle quasiparticle truncation at the FAM level is less clear.Writing (9) explicitly, we have for example The question here is whether summations over both µ and ν should be truncated up to N qp , only one of them, or none.Some insight can be gained by examining the time-dependent evolution of the quasiparticle operator β µ (t) of TDHFB, Quasiparticle index µ should obey µ ≤ N qp to be consistent with the definition of the onebody density matrix (15).However, the expression for δβ µ (t) is not informative about any additional truncation on ν.Therefore, one may be allowed to let the row indices of the matrices X and Y span the whole quasiparticle space while their column indices must be truncated.We realize that this is an a minima argument that should call for further clarifications, perhaps by extending the analysis the behavior of the HFB densities ρ and κ in homogeneous nuclear matter at the limit of large momenta; see discussion in Chapter 4 of [20].In the rest of this proceeding, we will adopt the "1/2" truncation (only column indices are truncated).

Benchmark Results
In this section, we present a series of calculations aimed at benchmarking the convergence of FAM calculations of the response function as a function of different numerical input parameters.All calculations were performed in the 240 Pu nucleus with the SLy4 parameterization of the Skyrme force [23].In the pairing channel, we use V (n) 0 = −219.973MeV and V (p) 0 = −285.524MeV as in [24].HFB calculations were performed with the code HFBTHO [25] (axial symmetry) using a quadrature grid of N GH = N GL = 40 Gauss-Hermite and Gauss-Laguerre integration points and N Leg = 80 Gauss-Legendre integration points [26].
In this work, we considered both electric (EL) and magnetic (ML) operators.We looked at the convergence patterns of both isoscalar and isovector operators and, for each operator, at the K = 0 and K > 0 components.Results were very similar in all cases, so we will only show results for the isoscalar, electric quadrupole (E2) transition with K = 0 defined by the operator where the Y lm (θ, φ) are spherical harmonics [27] and we took e as the regular electric charge. LLNL-PROC-855567

Convergence of the Response Function with Respect to the HO Basis
Both the HFB and FAM equations are solved in the harmonic oscillator (HO) basis.Basis functions are characterized by an oscillator length b 0 , which is equivalent to an oscillator frequency ω 0 = ℏ/(mb 2 0 ).In realistic calculations of deformed nuclei, the basis itself is deformed, or stretched.For axially-deformed basis, this means that the two frequencies ω z and ω ⊥ are different, leading to different oscillator lengths b z b ⊥ .In HFBTHO, the ratios of the two frequencies is related to the basis deformation parameter β 2 ; see [26] for details.The oscillator is thus characterized by two parameters, b 0 and β 2 .The total number of basis functions is determined by the number of shells N sh .In the left panel of Fig. 1, we show how the isoscalar, electric quadrupole excitation strength changes with the size of the basis as determined by the total number of shells N sh .For each calculation, the oscillator length is b 0 = 2.3 fm and the basis is spherical, β 2 = 0.It is somewhat surprising to see that the first peak at ω ≈ 5 MeV is hardly changed, the height of the second peak at ω ≈ 10 MeV varies by less than 15%, and most of the variability of the results is in the small third peak at around ω ≈ 12 MeV.This rather weak dependency of the results with the number of shells should be contrasted with the rather strong dependency of the HFB energy; see the right panel of Fig. 1.

LLNL-PROC-855567
In the left panel of Fig. 2, we show a similar study for the variations of the response with the oscillator length b 0 .Here, in each calculation the number of shells is N sh = 20 and the basis is also spherical.Although there are some small differences in the height of each of the peaks, the response function is pretty insensitive to the value of b 0 .The right panel of the figure shows that the variations of the HFB energy with b 0 do not exceed 700 keV at N sh = 20.Similar conclusions apply when varying the deformation β 2 of the basis.

Impact of Quasiparticle Cutoff
As mentioned earlier, HFB and QRPA calculations with zero-range pairing forces require introducing a cutoff on the number of quasiparticles used to define the densities (15) and (16).In Fig. 3, we illustrate the dependence of the response function on the choice of the cutoff.
Given a cutoff value E cut , a consistent calculation requires first calibrating the strengths V (n)  0 and V (p) 0 of the pairing force on some experimental data.Here we perform this calibration at the HFB level by employing the same prescription as in [24].In practice, for each E cut we refitted the pairing strengths so that the HFB pairing gaps match the three-point oddeven mass differences ∆ (3) in 232 Th.After this calibration, the pairing channel of the FAM calculations is well defined and we can compute the perturbed densities (16) according to the prescription discussed in Section 2.2.2.In the left panel of Fig. 3, we show results for the same isoscalar electric quadrupole excitation strength as before for four values of the cutoff.Not surprisingly, calculations are undistinguishable from one another since recalibrating the pairing strengths imply that the underlying HFB solutions are nearly identical by construction.For completeness, we also show in the right panel the non-consistent case where the cutoff of the HFB solution is fixed, here at E cut = 60 MeV, while the cutoff of the FAM calculation (16) itself varies.In other words a single identical HFB calculation is used as the basis for all FAM calculations.Even in this not very realistic scenario, the response depends only weakly on the cutoff.

Dependence of the Response Function on the Smearing Width
Finally, we discuss the dependency of the results on the smearing width Γ of the FAM response.As mentioned earlier, this parameter is introduced because the response function LLNL-PROC-855567 (13) has poles corresponding to QRPA eigenmodes: computing the response at complex ω values makes calculations more efficient [7].At the same time, it is equivalent to computing the response ( 14) by summing explicitly over the solutions of ( 6) and folding the result with a Lorentzian distribution of width Γ.The net result is that Γ becomes a parameter of the calculation and its choice will have an impact on the results.From a computational perspective, large values of Γ make calculations faster and more stable.From a physics standpoint, semi-classical models of the giant dipole response are typically based on a fit of the photoabsorption cross section with a Lorentzian shape [28].Surveys of experimental results suggest a rather large variation of Γ with the mass number A of the nucleus [1].In Fig. 4 we show how the isoscalar electric quadrupole excitation strength changes with Γ.As expected, larger values of Γ smooth out all the features of the response while smaller values make its fine structure more visible.Out of all the input parameters we considered in this work, Γ is by far the one with the largest impact on the results.

Conclusions
We have investigated the convergence patterns of the electromagnetic response function as a function of the numerical parameters controlling the calculation.Overall, the dependency of the results HO basis parameters is weaker than that of the underlying HFB calculations themselves.The response function is also insensitive to the quasiparticle cutoff -when the latter is implemented consistently at both HFB and FAM levels.Only the smearing width Γ used to define the complex frequency ω of the response has a visible impact on the results.

Figure 1 .
Figure 1.Left: Isoscalar electric quadrupole excitation strength dB(ω)/dω in 240 Pu as a function of the excitation energy ω, for different values of the total number of shells in the basis.Right: Variation of the HFB energy across the same range of oscillator shells.

Figure 2 .
Figure 2. Left: Isoscalar electric quadrupole excitation strength dB(ω)/dω in 240 Pu as a function of the excitation energy ω, for different values of the oscillator length (in fm).Right: Variation of the HFB energy across the same range of oscillator lengths.

Figure 3 .
Figure 3. Left: Isoscalar electric quadrupole excitation strength dB(ω)/dω in 240 Pu as a function of the excitation energy ω, for different quasiparticle cutoff energies.For each E cut , the strengths V (n) 0 and V (p) 0 of the pairing force have been refitted to keep the pairing gaps of 232 Th at their reference values.Right: Same as in the left panel, only a single HFB solution with E cut = 60 MeV is used for all four FAM calculations; see text for additional details.