Finite amplitude method in linear response TDDFT calculations

The finite amplitude method is a feasible and efficient method for the linear response calculation based on the time-dependent density functional theory. It was originally proposed as a method to calculate the strength functions. Recently, new techniques have been developed for computation of normal modes (eigenmodes) of the quasiparticle-random-phase approximation. Recent advances associated with the finite amplitude method are reviewed.


Introduction
The nuclear energy density functional methods are extensively utilized in studies of nuclear structure and reaction [1]. For the linear response calculations, the quasiparticle-random-phase approximation (QRPA), based on the time-dependent density-functional theory [2,3,4,5], is a standard theory [6]. However, for realistic energy functionals, it demands both complicated coding and large-scale computational resources. To resolve these issues, there have been several developments [7,8,9], including the finite amplitude method (FAM) [10].
The FAM allows us to avoid explicit evaluation of complex residual fields, which significantly reduces the necessary coding effort for realistic energy functionals. In addition, the use of iterative solution of the FAM equations also reduces the computational task and the memory requirement. The FAM was first adopted for the linear response calculations for the electric dipole mode, using the Skyrme energy functionals without the pairing correlations [11,12]. The method was soon extended to superfluid systems [13]. Then, the FAM was adopted in the code hfbtho to study monopole resonances in superfluid deformed nuclei with axially symmetry. Mario Stoitsov played a leading role in this project, which started during his visit to RIKEN in 2010 [14]. Later, there have been further developments in the FAM with relativistic and non-relativistic frameworks [15,16,17] The FAM was originally proposed for a feasible method to calculate the response function of a given one-body operator. Recently, new methodologies have been developed, for calculating discrete eigenstates using the FAM. In this article, we review some of these recent developments.

The finite amplitude method
In this section, we briefly illustrate the essential idea of the FAM. Let us start from the Hartree-Fock-Bogoliubov (HFB) equations, where H, E i > 0, and Φ i = U i V i are the HFB Hamiltonian including the particle-number cranking term (−µN ), the quasiparticle energies, and the quasiparticle states associated with the HFB ground state, respectively.
are conjugate to Φ i with negative quasiparticle energies [6]. They are orthonormalized as The generalized density matrix R [6] at the ground state can be written in a simple form as In the small-amplitude (linear) approximation, the density fluctuation δR can be decomposed into normal modes δR (n) , In the linear order, the matrix elements of δR (n) between Φ i andΦ j are denoted by X The other matrix elements of δR vanish in the quasiparticle basis: Note that δR (n) are in general non-Hermitian, nevertheless δR(t) is Hermitian. The density fluctuation δR (n) induces an residual field in the HFB Hamiltonian, H = H + δH (n) , where δH (n) linearly depends on δR (n) . The essential idea of the FAM is that, instead of explicitly expanding δH (n) in the linear order in δR (n) , we compute δH (n) by the finite difference using a small parameter η as The Hamiltonian H[R+ ηδR (n) ] should be evaluated at the density, i , using the following quasiparticles For a given set of forward and backward amplitudes, {X (n) , Y (n) }, the calculation of Eq. (7) is relatively easy because all we need to calculate are the one-body quantities with two-quasiparticle indices, such as δH In contrast, the QRPA matrices have four-quasiparticle indices [6], A ij,kl and B ij,kl . The calculation of these matrix elements demands significant coding effort.

Calculation of the QRPA matrices
In this section, we recapitulate the matrix-FAM (m-FAM) proposed in Ref. [15]. This provides a simple numerical method to calculate the QRPA matrices using the principle of the FAM, Eqs. (7), (8), and (9).
The QRPA equation in the matrix form is given by H Z n = ω n N Z n , where The most demanding part is the calculation of the matrix elements, A ij,kl and B ij,kl , which are formally defined by the derivative of the HFB Hamiltonian with respect to the density. However, since the FAM allows us to evaluate H Z for a given vector Z, these matrix elements are provided by the FAM in the following way. Let us define the "forward" unit vectorê kl as X ij = δ ik δ jl and Y ij = 0, and the "backward" oneẽ kl as Y ij = δ ik δ jl and X ij = 0. Namely, these unit vectors arê Then, it is trivial to see that the upper component of (H e n ) ij and (Hẽ n ) ij are identical to A ij,kl and B ij,kl , respectively.
On the other hand, using the FAM, (Hê kl ) up reads where Φ † i δHΦ j in the right hand side can be computed according to Eq. (7). Here, R + ηδR is given by the ground-state quasiparticles, Ψ † i = Φ † i and Ψ ′ i = Φ i , except for i = k and l, Following the same procedure with the "backward"ẽ kl , we can calculate B ij,kl . The numerical coding of the present method is extremely easy. All we need to calculate is the HFB Hamiltonian at the density R + ηδR which is defined by Eq. (14). After constructing the QRPA matrix, the QRPA normal modes of excitation are obtained by diagonalizing the QRPA matrix [6].

Iterative FAM with a contour integral in the complex frequency plane
In this section, we present a contour integral technique combined with the FAM developed in Ref. [17].
The iterative FAM (i-FAM) with a complex frequency ω provides a solution of the linear response equation (H − ωN ) Z(ω) = −F, where F is the one-body external field [10,13].
Here, we assumeF is a Hermitian operator. The QRPA response function is calculated as [10] This can be expressed in terms of the QRPA normal modes as This shows that the transition strength for the n-th normal mode, | n|F |0 | 2 , is a residue at ω = ω n . Therefore, if we choose the contour C n in the complex ω-plane that encloses ω = ω n > 0, it reads The corresponding QRPA normal modes, X (n) and Y (n) , are also given by the contour integral as [17] X (n) The contour integral method with i-FAM is complementary to the m-FAM. In the present approach, the eigenmodes are obtained by solving the i-FAM equations for complex frequencies combined with the contour integral. In the m-FAM, we do not resort to an iterative algorithm to solve the FAM equations, but we need to diagonalize the QRPA matrix at the end. For a small model space, the m-FAM has a significant advantage, however, the computational task of the m-FAM strongly depends on the size of the QRPA matrix D, which typically scales as D 3 .