{'fraction_non_alphanumeric': 0.05526376946316523, 'fraction_numerical': 0.015338136184057634, 'mean_word_length': 4.0853226187662495, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite amplitude method for the quasi-particle-random-phase approximation

We present the finite amplitude method (FAM) for superfluid systems. A Hartree-Fock-Bogoliubov code may be transformed into a code of the quasi-particle-random-phase approximation (QRPA) with simple modifications. This technique has advantages over the conventional QRPA calculations, such as coding feasibility and computational cost. We perform the fully self-consistent linear-response calculation for a spherical neutron-rich nucleus 174Sn, modifying the HFBRAD code, to demonstrate the accuracy, feasibility, and usefulness of the FAM.

Recently, there has been a renewed interest in the solution of the RPA problem [17][18][19]. In Ref. [17], the finite amplitude method (FAM) was proposed as a feasible method for a solution of the RPA equation. The FAM allows to calculate all the induced fields as a finite difference, employing a computational program of the static mean-field Hamiltonian. Recently, the FAM has been applied to the electric dipole excitations in nuclei using the Skyrme energy functionals [18]. There has been also a calculation making use the iterative Arnoldi algorithm for a solution of the RPA equation [19]. These newly developed technologies in conjunction with fast solving algorithms for linear systems open the possibility to explore systematically the nuclear excitations over the entire nuclear chart.
So far, these new techniques [17][18][19] have been developed for solutions of the RPA without the pairing correlations. It is well known, however, that almost all but magic nuclei display superfluid features. Therefore, a further improvement is highly desirable to make these methods applicable to the QRPA equations including correlations in the particle-particle and hole-hole channels. The purpose of the present paper is to generalize the FAM to superfluid systems, which enables us to perform a QRPA calculation utilizing a static Hartree-Fock-Bogoliubov (HFB) code with minor modifications. Our final goal would be the construction of a fast computer program for the fully self-consistent and triaxially deformed QRPA. This paper is a first step toward the goal, to present the basic equations of the FAM for the QRPA and show the first result for spherical nuclei. We use the spherically symmetric HFB code called HFBRAD [20] to be converted into the QRPA code.
This paper is organized as follows: In Sec. II, the QRPA equation is derived as the small-amplitude limit of the time-dependent HFB (TDHFB) equations. In Sec. III, we obtain the FAM formulae for the calculation of the induced fields. In Sec. IV, we summarize all the relevant formulae for practical application of the FAM. In Sec. V, we apply the FAM to the HFBRAD and compare the result with that of another self-consistent calculation. Sec. VI is devoted to the conclusions.

II. SMALL AMPLITUDE LIMIT OF THE TDHFB
In this section, we recapitulate the basic formulation of the TDHFB and its small-amplitude limit. In general, we will follow the notation in Ref. [1] unless otherwise specified. We also use = 1 in the following equations.
We start from the energy functional E[ρ, κ, κ * ] which is a functional of the density matrix and pairing tensor.
where |Φ is the HFB state. The single-particle Hamiltonian h and the pairing potential ∆ are obtained with a variation of the energy functional with respect to ρ and κ * , respectively.
The Bogoliubov quasi-particles, (a µ , a † µ ), have a linear connection to the bare particles, Here, the index k indicates the adopted basis such as the harmonic oscillator states or the coordinate space. The quasi-particles a µ are chosen so as to diagonalize the HFB Hamiltonian [1].
(3) Here, the normal ordering is assumed.
In a similar manner, the time-dependent quasiparticles a † µ (t) are characterized by the time-dependent wave functions (U (t), The time evolution of the quasi-particles under a one-body external perturbation F (t) are determined by the following TDHFB equation.
where the TDHFB Hamiltonian is given by Here and hereafter, the constant shift is neglected, since it does not play any role in the TDHFB equation (4). h(t) and ∆(t) become time-dependent, since they depend on the densities, ρ(t) = V * (t)V T (t) and κ(t) = V * (t)U T (t) = −U (t)V † (t), which are time-dependent. Note that the static quasi-particles correspond to a quasistatic solution of Eq. (4), a µ (t) = a µ e iEµt , with F (t) = 0.
Let us assume that the nucleus is under a weak external field of a given frequency ω.
where A † µν ≡ a † µ a † ν and B µν ≡ a † µ a ν . A small real parameter η is introduced for the linearization. In the smallamplitude limit, the second term (B-part) in Eq. (7) can be omitted, because it doesn't contribute in the linear approximation. The Bogoliubov transformation of the external fields (F 20 µν (ω) and F 02 µν (ω)) is given in Appendix A 2.
The external perturbation F (t) induces a density oscillation around the ground state with the same frequency ω. The density oscillation, then, produces the induced fields in the single-particle Hamiltonian, h(t) = h 0 +δh(t) and in the pair potential, ∆(t) = ∆ + δ∆(t). Thus, the Hamiltonian, Eq. (5), is decomposed into the static and oscillating parts; H(t) = H 0 + δH(t).
The time-dependent quasi-particle operators are decomposed in a similar manner: where δa µ (t) can be expanded in the quasi-particle creation operators: It should be noted that δa µ can be expanded only in terms of the creation operators, because the annihilation operators in the right-hand side of Eq. (11) simply represent the transformation among themselves, a µ (t) = ν C µν (t)a ν , and do not affect ρ and κ. The amplitudes, X and Y , must be anti-symmetric to satisfy the fermionic commutation relation, {a µ (t), a ν (t)} = 0. Keeping only linear terms in η, Eq. (4) becomes (12) Substituting Eqs. (6)-(11) into Eq. (12), we obtain the linear-response equations: In Eq. (13), setting the frequency complex, ω → ω+iγ/2, we can introduce a smearing with a width γ.
Expanding δH 20 (ω) and δH 02 (ω) in terms of the forward and backward amplitudes, X and Y , we obtain a familiar expression of the equation [1]: This matrix formulation requires us to calculate the QRPA matrix elements, A µν,µ ′ ν ′ and B µν,µ ′ ν ′ . This is a tedious task and their dimension, which is equal to the number of two-quasi-particle excitations, becomes huge especially for deformed nuclei. Instead, in the FAM [17], we keep the form of Eq. (13) and calculate the induced fields δH 20 (ω) and δH 02 (ω) using the numerical differentiation. We explain this trick in the next section.

III. FINITE AMPLITUDE METHOD FOR THE INDUCED FIELDS
The expressions for δH 20 and δH 02 in Eq. (13) are given in Eqs. (A8) and (A9), respectively. Thus, we need to calculate δh(ω) and δ∆ (±) (ω) for given X and Y . We perform this calculation following the spirit of the FAM [17].
From Eqs. (10) and (11), we obtain the timedependent quasi-particle wave functions: where First, let us discuss how to obtain δh(ω). The timedependent single-particle Hamiltonian h(t) depends on the densities which are determined by the wave functions (U (t), V (t)). Therefore, h(t) can be regarded as a functional of wave functions as Here, it should be noted that the phase factors, e iEµt in Eq. (15), do not play a role. This is because h is a functional of densities, ρ, κ, and κ * , which are given by products of one of (U, V ) and one of the complex conjugate (U * , V * ), such as ρ = V * V T and κ = V * U T . Therefore, the time-dependent phases in Eq. (15) are always canceled, thus can be omitted. Now, we take the small-amplitude limit, keeping only the linear order in η.
Here, δh(ω) can be obtained using Eqs. (16) and (17), expanding up to the first order in η and collecting terms proportional to e −iωt , as The calculation of the derivatives, such as ∂h kl /∂U * k ′ µ , is a tedious task and requires a large memory capacity for their storage in the computation. In the FAM, we avoid this explicit expansion, instead write the same quantity as follows: This is the FAM formula for the calculation of δh(ω).
All we need in the computer program is a subroutine to calculate the single-particle Hamiltonian as a function of For the pair potential, basically, the same arguments lead to the FAM formulae for δ∆ (±) . The time-dependent pair potential ∆(t) can be written as (23) Here, δ∆ (+) and δ∆ (−) are independent, since ∆(t) is non-Hermitian in general. δ∆ (+) can be written in the same form as Eq. (21).
whereŪ * η ,V * η , U η , and V η are given by Eq. (22). The expression for δ∆ (−) is also obtained from Eq. (23), collecting terms proportional to e iωt . It is given by the same expression as Eq. (24), However, (Ū * η ,V * η ; U η , V η ) here are different from Eq. (22) and given bȳ The essential trick of the FAM is to calculate the induced fields, δh(ω) and δ∆ (±) , according to Eqs. (21), (24), and (25) with a small but finite parameter η. Of course, the η 2 and higher-order terms bring some numerical errors, but they are negligible. Therefore, for given X and Y , we are able to calculate these induced fields, by using the static HFB code with some minor modifications. δH 20 (ω) and δH 02 (ω) of Eq. (13) in the quasi-particle basis can be calculated with Eqs. (A8) and (A9), respectively. Then, we may solve the QRPA linearresponse equation (13) to obtain the self-consistent amplitudes, X and Y , utilizing an iterative algorithm (See Sec. IV).

A. Induced fields in terms of densities
Although the basic formulae of the FAM has been provided in Sec. III, we may need to modify them in the practical implementation of the FAM. For instance, some HFB codes, such as HFBRAD, contain subroutines to calculate mean fields as functions of densities, not of wave functions. In this subsection, we rewrite Eqs. (21), (24), and (25) in terms of densities.
The density δρ(t) is written up to linear order in η as where This can be written in the FAM form: whereV * η and V η are given in Eq. (22). The pair tensor κ(t), which is non-Hermitian, can be expressed in a similar manner.
Here, κ (±) can be given in the explicit form as and in the FAM form as whereV * η and U η are given in Eq. η . Now, let us present how to obtain the induced fields in terms of the densities. In general, h(t) and ∆(t) may depend on ρ, κ, and κ * .
(34) In order to obtain the induced fields, all we need to do is to replace ρ by ρ η defined in Eqs. (29), and κ by κ (±) η in Eq. (33), as follows: where the terms of the second and higher orders in η are neglected.

IV. SUMMARY OF THE FINITE AMPLITUDE METHOD
Here we provide a summary of the FAM for the QRPA linear-response calculation for a prompt application. Later, we discuss applications of the FAM to the Skyrme functionals, however, the FAM formulated in this and previous sections is applicable to any kind of energy density functional (mean-field) models.

A. Numerical procedure
The aim is to solve the linear-response equation (13) for a given external field F . In order to obtain the forward and backward amplitudes, X and Y , we resort to an iterative algorithm. Namely, we start from the initial guess for (X, Y ) = (X (0) , Y (0) ) ≡ x (0) , and calculate δh(ω) and δ∆ (±) (ω) according to the formulae, (21), (24), and (25). Then, they are converted into δH 20 (ω) and δH 02 (ω), using Eqs. (A8) and (A9), respectively. In this way, we can evaluate the left and right hand sides of Eq. (13) for a given (X, Y ).
Since Eq. (13) is equivalent to Eq. (14), it is a linear algebraic equation for the vector x ≡ (X, Y ), in the form of A x = b. Many different algorithms are available for the solution of linear systems. In this paper, we resort to a procedure based on Krylov spaces called generalized conjugate residual (GCR) method [21]. Within these kinds of methods, a succession of approximate solutions ( x (0) , x (1) , x (2) , · · · ) converging to the exact one is obtained by the iteration. The GCR algorithm consists in a series of steps each containing the operation of the matrix A on a given vector, and sums and scalar products of two vectors. For the given x = (X, Y ), A x is equal to the left hand side of Eq. (13). Therefore, the quantity A x can be calculated without the explicit knowledge of the QRPA matrix itself.
Here, we summarize the formulae. The linear response equation is given by A Denoting h and ∆ collectively as H ≡ (h, ∆), the induced fields δH are calculated by the FAM formulae, for the calculation of δh(ω) and δ∆ (+) . For δ∆ (−) , they areŪ * The final result does not depend on the parameter η, as far as it is in a reasonable range. The choice of η is discussed in Sec. V.

B. Calculation of the strength function
Using the solution (X, Y ), we can calculate the strength function following the same procedure as Ref. [17].
(39) Here, S(F ; ω) is obtained from the solution (X, Y ). For the operator in the form of Eq. (A10), we may calculate S(F ; ω) as For the operator in the form of Eq. (A13), we have For both cases, in the two-quasi-particle basis, Eqs. (40) and (41) can be written in the unified expression.
where F 20 and F 02 are given by Eqs. (A11) and (A12) for the former case, and by Eqs. (A14) and (A15) for the latter.

V. APPLICATION OF THE FAM TO THE HFBRAD
In order to assess the validity of the FAM, we install the FAM in the HFBRAD code [20]. It has to be noted that the formalism of the HFBRAD is slightly different from the one used in this paper which follows the notations in Ref. [1]. In particular, the wave functions (ϕ 1µ , ϕ 2µ ), the pairing tensorρ, and the pair potentialh are defined in a different manner; ϕ 1µ (k) = U kµ , ϕ 2µ (k) = Vk µ ,ρ kl = κ kl , andh kl = ∆ kl , wherek is the time-reversal state of k. A detailed discussion on the difference among the two notations can be found in Ref. [25].
The HFBRAD [20] is a well known code which solves the HFB in the radial coordinate space assuming the spherical symmetry. It has been designed to provide fast and reliable solutions for the ground state of spherical even-even nuclei. For these nuclei, the time-odd densities are identically zero and thus they have not been implemented in the code. In order to render the QRPA fully self-consistent, we have to add the time-odd terms in the calculation of the induced fields. This task can be simplified for a case of the presence of spherical and spaceinversion symmetry, such as in the case of monopole excitations. For this case, the only time-odd terms with non-zero contribution are those due to the current density [22], moreover the only non-vanishing component of the current density is radial.
We calculate the strength function of the isoscalar monopole for a neutron-rich nucleus, 174 Sn. To check the self-consistency by looking at the spurious component, we also calculate the strength of the nucleon number operator. Both operators are given by the form of Eq. (A10) with f kl = k|r 2 |l for the isoscalar monopole operator and f kl = δ kl for the number operator.
In order to obtain the strength function, first, we have to solve the HFB equations to construct the ground-state wave functions (U, V ). It is accomplished by using the HFBRAD code. The parameters of the present calculation are adjusted to the values used by Terasaki and co-workers in Ref. [6]; The box size is R box = 20 fm, the quasi-particle energy cutoff is E c qp = 200 MeV, the maximum angular momenta of the quasi-particle states are j n max = 21/2 for neutrons, and j p max = 15/2 for protons. We use the Skyrme functional with the SkM* parameter set [23] in the ph-channel and a delta interaction of the volume type with the strength V 0 = −77.5 MeV fm 3 for the pp-and hh-channels. The next step is solving the linear-response equation for a given external field of the frequency ω. At first, we build the induced fields, δh(ω) and δ∆ (±) (ω), starting from a guess choice of the QRPA amplitudes (X (0) , Y (0) ), according to Eq. (38). In the present calculation, we choose either X (0) = Y (0) = 0 or the values of X and Y at the previous energy ω calculated. We resort to the iterative algorithm of the GCR method to solve the equation (13). We include all the two-quasi-particle states (µν) within the HFB model space defined above (E µ(ν) ≤ 200 MeV). The two-quasi-particle space amounts to 12,632 states for J π = 0 + . Note that this number becomes much larger if we treat deformed systems. We set the accuracy of the convergence to be ǫ < 10 −5 , where The number of iterations needed depends on ω; at low energies, about 50-60 iterations are enough to reach the convergence, while, close to the central peak at 12 MeV, more than 300 iterations are needed.
We studied the convergence quality of the solutions as a function of the parameter η used for the numerical derivative. This is shown in Table I. If η is too big (η ≥ 10 −4 ) the derivative of the FAM becomes inaccurate and the linearity of the procedure is partially broken. The residue A x − b reaches a plateau where increasing the number of iterations cannot improve it anymore. For 10 −5 ≤ η ≤ 10 −9 , the calculations converge well and the resulting strength function is stable. If η becomes smaller than 10 −10 , the numerical precision limits are reached and the GCR procedure can no longer obtain the required precision. Therefore, we may conclude that the parameter η in the range of 10 −5 ≤ η ≤ 10 −9 is appropriate to obtain the induced fields accurately. Although the constant value η = 10 −8 is adopted in this paper, we may use a more sophisticated choice, such as the ω-dependent η values [17,18], We report the strength function of the isoscalar monopole mode. To smear the strengths at discrete eigenenergies, we add an imaginary term to the energy: ω → ω + iγ/2, where γ = 1.0 MeV. This procedure is almost equivalent to smearing the strength function with a Lorentzian function with a width equal to γ. The calculated energy-weighted strengths are summed up to 300 MeV and we found that they exhaust 99.6 % of the theoretical sum-rule value given by 2 m A r 2 . In Fig. 1, we compare our results (solid red curve) with the one in Ref. [6] (dashed green curve). The selfconsistent result obtained by Terasaki et al. [6] also employs the HFB solutions calculated with the HFBRAD. However, in Ref. [6], the QRPA matrix is calculated in the canonical-basis representation and an additional truncation of the two-quasi-particle space is introduced for the construction of the QRPA matrix . In contrast, we introduce no additional truncation for our FAM calculation. We compare our results with the one of the cutoff (iii) in Ref. [6] which takes into account the highest number of states for the construction of the QRPA matrix; all the proton quasi-particles up to 200 MeV and the neutron canonical levels with occupancy v 2 > 10 −16 .
In the first two peaks at E ∼ 5 and 8.5 MeV, the two curves are almost perfectly overlapping. The peaks between 11 MeV and 18 MeV occur at the same energy for the two calculations while their height is slightly different. The bump close to zero energy resulting in our calculations has to be attributed to the presence of a spurious mode. To check the position of the spurious mode related to the pairing rotation of the neutrons, we included in Fig. 1 the transition strength associated to the number operator, by the blue dashed line. The spurious mode is well localized close to zero energy.
The present result demonstrates the accuracy and usefulness of the FAM for the superfluid systems. Even if the two codes include some differences in the truncation of the two-quasi-particle space, the similarity of the results is very satisfying

VI. CONCLUSIONS
The finite amplitude method for the QRPA has been presented. The basic idea is identical to the original FAM [17], that we resort to a numerical differentiation to calculate the induced fields and then solve the linear-response equation with an iterative algorithm such as the GCR. With the FAM, a HFB code with simple modifications can be turned into a QRPA code. Especially, it is very easy to construct the QRPA code which has the same symmetry of the parent HFB one whose subroutines are used to perform the numerical derivative. All the terms present in the TDHFB calculation, including the timeodd mean fields, should be taken into account to construct fully self-consistent codes. This requires us some effort to update the original HFB code. Still, the necessary task for coding the FAM is much less than that for the explicit calculation of the QRPA matrix elements for realistic energy functionals. In addition, it does not require a large memory capacity, since we do not construct the QRPA matrix. We have built a fully self-consistent QRPA code using the HFBRAD [20]. The iterative algorithm, for which we adopted the GCR method in this paper, may be replaced by a better algorithm in future. The resulting strength functions of the isoscalar 0 + mode of 174 Sn show high similarity with the fully self-consistent calculations in Ref. [6]. Thus, this paper showed the first application of the FAM for superfluid systems and demonstrated the usefulness of the FAM for the construction of the QRPA code by modifying existing HFB codes. non-Hermitian in general. The Bogoliubov transformation as in Eq. (A7), then, leads to F 20 and F 02 in Eq. (7), In case that F (ω) has a form of pairing-type the same calculation provides F 20 and F 02 by