Linear response strength functions with iterative Arnoldi diagonalization
J. Toivanen, B.G. Carlsson, J. Dobaczewski, K. Mizuyama, R.R. Rodriguez-Guzman, P. Toivanen, P. Vesely
aa r X i v : . [ nu c l - t h ] D ec Linear response strength functions with iterative Arnoldi diagonalization
J. Toivanen, B.G. Carlsson, J. Dobaczewski,
1, 2
K. Mizuyama, R.R. Rodr´ıguez-Guzm´an, P. Toivanen, and P. Vesel´y Department of Physics, University of Jyv¨askyl¨a, P.O. box 35, FIN-40014, Finland Institute of Theoretical Physics, Warsaw University, ul. Ho˙za 69, PL-00681, Warsaw, Poland (Dated: November 7, 2018)We report on an implementation of a new method to calculate RPA strength functions withiterative non-hermitian Arnoldi diagonalization method, which does not explicitly calculate andstore the RPA matrix. We discuss the treatment of spurious modes, numerical stability, and how themethod scales as the used model space is enlarged. We perform the particle-hole RPA benchmarkcalculations for double magic nucleus
Sn and compare the resulting electromagnetic strengthfunctions against those obtained within the standard RPA.
PACS numbers: 21.60.Jz, 71.15.MbKeywords: Random Phase Approximation, Skyrme energy density functional, multipole strength functions
I. INTRODUCTION
The linear response theory (LRT) obtained from thelinearized time-dependent mean field method is an im-portant tool for calculating properties of excited statesof many-fermion systems, such as nuclear giant reso-nances. In its charge-changing version, it can also giveaccess to the beta-decay strengths. This method is espe-cially important in heavy nuclei, where the shell-model orconfiguration-interaction approaches are intractable. Anadvantage is also that LRT does not require the knowl-edge of an interaction and can therefore be used bothwithin density functional theory (DFT) and phenomeno-logical energy density functional (EDF) approaches, giv-ing rise to a set of equations of RPA type. Below, forsimplicity we refer to this method and associated equa-tions simply as RPA method/equations. Strength func-tions obtained in this way probe new aspects of the EDFsand thus have a potential of constraining parameters inphenomenological nuclear EDFs.The purpose of the present study is to present an im-plementation of an efficient RPA algorithm that is basedon the local nuclear EDF. For electronic systems, simi-lar methods have been used since many years, see, e.g.,the recent Ref. [1] for a review, and they also consti-tute parts of standardized computer packages such asGAMESS [2, 3]. There are two essential elements ofthese methods, which are at the heart of their efficiencyand scalability, namely, (i) the RPA equations are solvediteratively and (ii) the RPA matrix does not have to beexplicitly calculated. The second of these elements is par-ticularly important; it is based on the observation thatthe action of the RPA matrix on the vector of RPA am-plitudes can proceed through the calculation of the meanfields corresponding to these amplitudes.In nuclear physics context, probably the first studythat used the concept of mean fields in the RPA methodwas that by P.-G. Reinhard [4]. Iterative solutions ofthe RPA equations were introduced by Johnson, Bertsch,and Hazelton [5], and applied to the case of separable interactions, but in fact these methods can also be ap-plied in more complicated situations, as we show here.Strangely enough, these very efficient methods have notyet been used in practical applications. Only very re-cently, Nakatsukasa et al. [6, 7] have implemented theanalogous approach within the so-called finite amplitudemethod (FAM).Our present implementation pertains to the sphericalsymmetry with neglected pairing correlations – thus itconstitutes only a proof-of-principle study. The real chal-lenge is in solving the quasiparticle RPA (QRPA) prob-lem in deformed nuclei. Although at present, a few imple-mentations that are based on solving the standard QRPAequations already exist [8, 9] or begin to emerge [10, 11],such a route is bound to be blocked by the shear di-mensionality of the problem. On the other hand, as weshow here, methods based on the iterative solutions usingmean fields have much better scalability properties andare potentially very promising.The paper is organized as follows. In Secs. II and IIIwe lay down the essential features of the method by pre-senting the use of mean fields and iterative solution of theRPA method, respectively. Then, in Sec. IV we presentthe method to remove the spurious RPA states, which istailored to be used within the iterative approach. Sec-tions V and VI present the convergence and scalabilityproperties of our method, respectively, and summary andconclusions are given in Sec. VII.
II. RPA FROM LINEARIZED TDHF
To be concise, below we present a less general deriva-tion than the standard method [12, 13] to derive the RPAequations from linearized time-dependent Hartree-Fock(TDHF) equations. For density independent forces orfunctionals with terms quadratic in density, the densitymatrix and mean field of a time-dependent nuclear stateare expressed as ρ ( t ) = ρ + ˜ ρ ω e iωt + ˜ ρ † ω e − iωt , (1) h ( t ) = h + ˜ h ω e iωt + ˜ h † ω e − iωt , (2)where ρ and h are the Hartree-Fock (HF) ground-statedensity matrix and mean field, respectively. Inserting ρ ( t ) and h ( t ) of Eqs. (1) and (2) into the TDHF equa-tion, and keeping only terms linear in the fluctuatingquantities ˜ ρ and ˜ h , we get a linearized TDHF equation,or the RPA equations:¯ hω ˜ ρ ω,mi = ( ǫ m − ǫ i ) ˜ ρ ω,mi + ˜ h ω,mi , (3)¯ hω ˜ ρ ω,im = ( ǫ i − ǫ m ) ˜ ρ ω,im − ˜ h ω,im , (4)where we use the letter m for particle states and i forhole states, and where ǫ m,i are the HF single-particleenergies. The fields ˜ h are the first functional derivativesof the used EDF, evaluated using the density amplitudes˜ ρ of Eq. (1). Density dependence of the used EDF beyondquadratic gives rise to rearrangement fields in ˜ h . Theserearrangement parts of ˜ h must be linearized around ρ to make our RPA equations explicitly first order in ˜ ρ .One way to achieve this is by calculating functionalderivatives of the rearrangement parts of ˜ h with respectto density, which technically makes our mean-field rou-tine differ from the standard HF routines. Since in ourimplementation we use the standard Skyrme forces thathave simple density dependencies of the coupling con-stants, the explicit functional differentiation does notcause any mathematical or performance problems. Hada EDF with more complex density dependence been usedit would have been an advantage to instead use the FAMmethod [6] for linearization.If the matrix elements of ˜ h in Eqs. (3) and (4) areexpanded in terms of the particle-hole (p-h) and hole-particle (h-p) matrix elements of ˜ ρ , we obtain the tra-ditional RPA equations. In this work, we do not con-struct the RPA matrix, but directly solve Eqs. (3) and(4) by calculating the matrix elements of fields ˜ h us-ing a HF mean-field routine that uses the time-reversal-invariance breaking density matrix ˜ ρ . Since the sameroutine is used to evaluate the HF and RPA mean fields,the method is always fully self-consistent [14, 15]. Inthe following equations, we use the standard abbrevia-tions X mi = ˜ ρ ω,mi and Y mi = ˜ ρ ω,im . The density vectorthat contains the p-h matrix elements of ˜ ρ ω is definedas X ω = ( X m ,i , X m ,i , . . . , X m D ,i D ), and similarly forthe vector Y of h-p elements, where D is the number ofallowed p-h configurations. Overlaps of RPA vectors aredefined as h X, Y | X ′ , Y ′ i = (cid:0) X ∗ , Y ∗ (cid:1) (cid:18) X ′ T −Y ′ T (cid:19) (5)and the minus sign results from the RPA norm matrix. III. ITERATIVE SOLUTION OF THE RPAEQUATIONS
The RPA equations (3) and (4) constitute a non-hermitian eigenproblem with non-positive-definite norm. We solve this problem by using an iterative method thatduring each iteration only needs to know the product ofthe RPA matrix and a density vector, that is, the right-hand sides of Eqs. (3) and (4): W kmi = ( ǫ m − ǫ i ) X kmi + ˜ h mi ( X k , Y k ) , (6) W ′ kmi = ( ǫ i − ǫ m ) Y kmi − ˜ h im ( X k , Y k ) , (7)where index k labels iterations and the mean fields˜ h ( X k , Y k ) depend linearly on the density vectors X k and Y k . Expressed through the standard RPA matrices A and B [12], Eqs. (6) and (7) for a positive norm basisvector and for its opposite norm partner vector read: (cid:18) W k + W ′ + k (cid:19) = (cid:18) A B − B ′∗ − A ′∗ (cid:19) (cid:18) X k Y k (cid:19) , (8) (cid:18) W k − W ′− k (cid:19) = (cid:18) A B − B ′∗ − A ′∗ (cid:19) (cid:18) Y k ∗ X k ∗ (cid:19) . (9)In exact arithmetic A = A ′ and B = B ′ and thereforeeither Eqs. (8) or (9) could be used in the iteration pro-cedure with equivalent results. Nevertheless, below weuse them both to stabilize the iteration process.Various iterative methods, which only need to knowthe products of the diagonalized matrix and vectors, ex-ist for non-hermitian matrix eigenvalue equations, andgood examples with pseudocode are shown in Ref. [1]. Wechose the non-hermitian Lanczos method [5] in a modifiedform, because it conserves all odd-power energy weighedsum rules (EWSR) if the starting vector (pivot) of it-eration is chosen correctly. In this work, we start froma pivot vector that has its elements set to the matrixelements of electromagnetic multipole operator, X mi = e √ N h φ m | r p Y JM | φ i i , Y mi = 0 , (10)where p = 2 and J = 0 for the 0 + mode, p = 1 and J = 1for the 1 − mode, and p = 2 and J = 2 for the 2 + mode.The constant N is used to normalize the pivot vectorto unity. For the IS 1 − mode we only present resultsobtained with the operator:( r − h r i r ) Y M , (11)to stay consistent with Refs. [16] and [17]. For the choiceof pivot in (10) one can prove [5] that all odd-powerEWSRs are satisfied throughout the iteration procedure.Because we calculate the RPA matrix-vector productsby using the mean-field method, and not with a precal-culated RPA matrix, we introduce small, but significantnumerical noise to the resulting vectors. If correctivemeasures are not used to remove or reduce this noise,the iteration method fails and produces complex RPAeigenvalues early on in the iteration. We stabilize ouriterative solution method by modifying the method ofRef. [5] in two ways. First, we use the non-hermitianArnoldi method instead of the non-hermitian Lanczosmethod. The advantage of Arnoldi method is that itorthogonalizes each new basis vector against all previousbasis vectors and their opposite norm partners, that is, (cid:18) ˜ X k +1 ˜ Y k +1 (cid:19) = (cid:18) W k + W ′ + k (cid:19) − k X i =1 (cid:18) X i Y i (cid:19) a ik + k X i =1 (cid:18) Y i ∗ X i ∗ (cid:19) b ik , (12) (cid:18) ˜ Y k +1 ∗ ˜ X k +1 ∗ (cid:19) = − (cid:18) W k − W ′− k (cid:19) + k X i =1 (cid:18) X i Y i (cid:19) b ′∗ ik − k X i =1 (cid:18) Y i ∗ X i ∗ (cid:19) a ′∗ ik , (13)where the overlap matrices a ik , b ik , a ′∗ ik , and b ′∗ ik are calcu-lated as in Eq. (5). Again, in exact arithmetic, Eqs. (12)and (13) are equivalent and the lower matrix elements a ′∗ ik and b ′∗ ik in Eq. (13) are exact complex conjugates of theelements of the upper Krylov-space [18] RPA matrices.We assume this to keep the Krylov-space RPA matrix inthe standard form.In the Lanczos method, only the tridiagonal parts ofRPA matrices are calculated, and small changes in ba-sis vectors due to Lanczos re-orthogonalization (whichalways must be used to preserve orthogonality) do notshow up in the constructed RPA matrix. In the Arnoldimethod, these small but important elements outside thetridigonal part improve the stability as compared toLanczos.The norm of the obtained new residual vector inEq. (12) can be either positive or negative. We do notin practice use Eq. (13), which in exact arithmetic wouldduplicate the results of Eq. (12). Instead, we store onlythe positive-norm basis states and use a similar methodas in Ref. [5] to change sign of the norm in case thenorm of the residual vector in Eq. (12) is negative. Thus,explicitly, for the positive norm of the residual vector˜ N k +1 = h ˜ X k +1 , ˜ Y k +1 | ˜ X k +1 , ˜ Y k +1 i , we define the newnormalized positive-norm basis vector as X k +1 mi = 1 p ˜ N k +1 ˜ X k +1 mi , Y k +1 mi = 1 p ˜ N k +1 ˜ Y k +1 mi . (14)If ˜ N k +1 <
0, the new normalized positive norm basisvector is defined as X k +1 mi = 1 p − ˜ N k +1 ˜ Y k +1 ∗ mi , Y k +1 mi = 1 p − ˜ N k +1 ˜ X k +1 ∗ mi . (15)When maximum number of iterations has been madeor the iteration has been stopped, the generated Krylov-space RPA matrix, with dimension d << D , is diagonal-ized with standard methods, that is we solve: (cid:18) a b − b ∗ − a ∗ (cid:19) (cid:18) x k y k (cid:19) = ¯ hω k (cid:18) x k y k (cid:19) . (16)The approximate RPA solutions are then obtained bytransforming the Krylov-space basis vectors, X kmi = d X l =1 (cid:0) X lmi x kl + Y l ∗ mi y k ∗ l (cid:1) , (17) Y kmi = d X l =1 (cid:0) Y lmi x k ∗ l + X l ∗ mi y kl (cid:1) , (18)for all k = 1 , . . . , d , and these vectors are used to evaluatethe strength functions.Standard RPA method that constructs and diagonal-izes the full RPA matrix can ensure that the lower matri-ces in the RPA supermatrix are exact complex conjugatesof the upper matrices. Our mean-field method can havesmall differences in the implicitly used upper and lowerRPA matrices due to finite numerical precision. The con-sequence of this is that we will in general have a ′ ij = a ij and b ′ ij = b ij . This spoils the consistency of Eqs. (12)and (13) and can make the Arnoldi iteration to fail andproduce complex energy solutions.The numerical errors in the matrix-vector products canbe reduced by symmetrization. We thus calculate theRPA fields twice, first using the densities of a positivenorm basis vector ( X k , Y k ), and second using the den-sities of negative norm vector ( Y k ∗ , X k ∗ ). The two re-sulting vectors are subtracted from each other to get thefinal stabilized RPA matrix-vector product, (cid:18) W k W ′ k (cid:19) = 12 W k + − W ′− k ∗ W ′ + k − W k ∗− ! , (19)to be used in Eq. (12). Together, the Arnoldi iterationmethod and symmetrization of matrix-vector productsstabilize our mean-field based iterative diagonalization. IV. TREATMENT OF SPURIOUS RPA MODES
For the discussion of various spurious modes in theRPA method we refer the reader to, e.g., Ref. [13]. In thepresent study, we only consider spherical ground statesneglecting pairing correlations, so the only spurious exci-tation is generated by the total linear momentum. There-fore, the only affected RPA mode is the isoscalar 1 − mode. In traditional RPA calculations that constructand diagonalize the full RPA matrix, the spurious 1 − mode is typically removed after the RPA diagonalization.Often a modified transition operator (11) is used, whichhas the property of h HF | h ˆ F , ˆ P cm i | HF i = 0, as long asthe commutator is evaluated within a complete set of ba-sis states. In a finite model spaces of localized orbitalsthis relation is no more exactly valid, and the correctedoperator does not remove spurious components exactly.To remove the spurious isoscalar 1 − mode from ourphysical RPA excitations we use the same method asin Ref. [6], where the basis vectors are orthogonalizedagainst the spurious translational mode P and its conju-gate ”boost” operator R , which have the form:ˆ P µ = 1 √ X mi (cid:16) i ( φ m ||∇ || φ i ) (cid:2) c † m ˜ c i (cid:3) µ + h . c . (cid:17) , (20)ˆ R µ = 1 √ X mi (cid:16) ( φ m || r || φ i ) (cid:2) c † m ˜ c i (cid:3) µ + h . c . (cid:17) . (21)The spurious RPA vectors ( P , P ∗ ) and ( R , R ∗ ) containthe p-h and h-p matrix elements of Eqs. (20) and (21), re-spectively. Our method differs from that of Ref. [6] in thefact that we orthogonalize our basis during the Arnoldiiteration, which fits naturally with the iterative solutionmethod and guarantees that the obtained approximateRPA excitations have exact zero overlaps with spuriousmodes. This is equivalent to diagonalizing the full RPAmatrix in the subspace orthogonal to the spurious states.In our implementation, each generated new Arnoldi basisvector is orthogonalized as (cid:18) X k Y k (cid:19) phys. = (cid:18) X k Y k (cid:19) − λ (cid:18) PP ∗ (cid:19) − µ (cid:18) RR ∗ (cid:19) , (22)where the overlaps λ and µ are defined as λ = h R, R ∗ | X k , Y k ih R, R ∗ | P, P ∗ i , (23) µ = − h P, P ∗ | X k , Y k ih R, R ∗ | P, P ∗ i . (24)When more symmetries are broken, formulas equiva-lent to Eqs. (22)–(24) can be used to remove spuriouscomponents coming from each broken symmetry of themean field. V. CONVERGENCE OF STRENGTHFUNCTIONS
The iterative Arnoldi method is meaningful for the cal-culation of strength functions only if the number of itera-tions needed for accurate results is significantly less thanthe full RPA dimension. To study how many Arnoldi it-erations we need for good accuracy, we calculated electro-magnetic isoscalar (IS) and isovector (IV) strength func-tions [16] for doubly magic nuclei. All calculations wereperformed by implementing the RPA iterative solutionswithin the computer program HOSPHE [19], which solvesthe self-consistent equations in the spherical harmonic-oscillator (HO) basis. We studied both the convergenceof smoothed strength functions as a function of numberof Arnoldi iterations and as a function of the number ofHO shells.We used the same definitions of the 0 + , 1 − , and 2 + transition operators as in Ref. [17] and the Skyrme func-tional SLy4 of Ref. [20]. The function we used to smooththe strength functions was also the same as in [16], with R box = 20 fm. Because the HF ground state of Snis spherically symmetric, our approximate RPA phononshave good angular momentum. We tested the use of largebasis sets up to 40 HO shells. The HF ground state ener-gies were well converged for all double magic nuclei when25 HO shells were used. Below, we present the resultsonly for
Sn.
A. Convergence as a function of the number ofArnoldi iterations
1. The + strength functions Figure 1 shows the 0 + IS and IV smoothed strengthfunctions for
Sn calculated with 100 Arnoldi iterationscompared with the standard RPA results from Ref. [17].Agreement between the strength functions is excellent.The four panels of Fig. 2 show the convergence of thesmoothed strength functions of Fig. 1 as the number ofArnoldi iterations increases. The panels show differencesof the strength functions calculated at 20 iteration inter-vals. The convergence is quite satisfactory after 100–120iterations.
FIG. 1: The 0 + strength functions in Sn calculated byusing 25 HO shells and 100 Arnoldi iterations for the SLy4functional (solid lines), compared with the standard RPA cal-culation of Ref. [17] obtained for the SkM* functional (dashedlines).FIG. 2: Convergence of the
Sn 0 + strength functions ofFig. 1. Solid lines are for the IS and dashed lines are forthe IV strength functions. Each panel shows the difference oftwo strength functions, one with n iterations and the othercalculated with n −
20 iterations.
2. The + strength functions Figures 3 and 4 show similar results as Figs. 1 and2, but for the 2 + strength functions in Sn. As forthe 0 + case, the IS and IV strength functions from theArnoldi iteration agree very well with the strength func-tions of Ref. [17]. The convergence of strength functionsis as fast as for 0 + ; after 120 iterations the smoothedstrength functions change only by about 5%. We thushave to make only 120 iterations to calculate reason-ably accurate 2 + strength functions for the RPA problemwhose dimension is D = 1020. The large double spikesobserved in Fig. 4 below 10 MeV are due to the lowestRPA phonons, which by the smoothing procedure acquire ≃ FIG. 3: Similar to Fig. 1 but for the 2 + strength functions.All results were calculated for the SLy4 functional.FIG. 4: Similar to Fig. 2 but for the 2 + strength functions.
3. The − strength functions Figure 5 compares the 1 − strength functions of our it-erative method (solid line) with the strength functionsfrom Ref. [17] (dashed line). The solid line shows theIS strength calculated when the generated Arnoldi basisis orthogonalized against spurious mode during iterationand dotted line corresponds to similar iterative calcula-tion without orthogonalization. The low-lying state at0.72 MeV, which has a large overlap with the spurious IS1 − mode disappears when the orthogonalization methodof Eqs. (22)–(24) is used. Also for the 1 − strength func-tion, 100–120 Arnoldi iterations were needed to producereasonably accurate results, see Fig. 6. FIG. 5: Main panel: the 1 − strength functions in Sn,calculated using 100 Arnoldi iterations and with spurious ISmode removed (solid line), and results of the standard RPAfrom Ref. [17] (dashed line). Dotted line shows results of 140Arnoldi iterations without orthogonalization against the spu-rious IS mode. Inset: same as in the main panel, but forthe IV strength functions. All results were calculated for theSLy4 functional.
When no orthogonalization is made against the spu-rious mode, the obtained excitations contain small com-ponents of the spurious mode. This affects the physicalpart of the IS strength distribution, especially around 20–30 MeV. The standard RPA strength function of Ref. [17]has not been corrected for the spuriosity but only thestrength of the lowest-lying state that has a large over-lap with the spurious IS mode has been omitted. At 20–30 MeV, this strength agrees well with our uncorrectedstrength.The orthogonalization method improves the conver-gence of the strength function, because now the 8.3 MeV1 − excitation is lowest in energy and thus converges first.Without orthogonalization against the spurious mode, weneed 140 iterations instead of 100 to get acceptably con-verged strength function. FIG. 6: Similar to Fig. 4 but for the 1 − strength functions.The IV strength-function differences were multiplied by thefactor of 200 fm . B. Convergence in function of the number of HOshells
In Section V A, we showed our strength functions cal-culated with 25 HO shells. This was found to be satisfac-tory, and using more shells did not appreciably changethe obtained strength functions. The only effect of usingmore oscillator shells was that we needed to use slightlymore Arnoldi iterations to produce well converged re-sults. In the case of 40 shells, about twenty more itera-tions were needed, compared to calculations made with25 shells. In Figs. 7, 8, and 9, we show the convergence ofstrength functions for the 0 + , 2 + , and 1 − modes, respec-tively. Each panel shows the difference of two strengthfunctions obtained in the intervals of ∆ N = 4 HO shells,between N of 22 and 38. FIG. 7: Similar as Fig. 2 but for the convergence of the 0 + strength functions as a function of the number of HO shells N . These plots overstress the variations of strength func-tions in the sense that slight shifts of peaks create theoscillating patterns in the difference plots. To illustrate
FIG. 8: Similar as Fig. 4 but for the convergence of the 2 + strength functions as a function of the number of HO shells N .FIG. 9: Similar as Fig. 6 but for the convergence of the 1 − strength functions as a function of the number of HO shells N . this point, in Fig. 10 we show the 1 − strength functionscalculated for N = 22, 26, 30, 34, and 38 HO shells.Poor convergence of the IS surface mode creates someuncertainty in the position and width of the high-energybump. Larger bases should probably be used if convergedresults for this particular mode were required.As noted in Ref. [5], well before the maximum numberof iterations (equal to the RPA dimension D ) is reached,the iteratively generated RPA matrix in the Krylov spacecan become singular. In that case, the stabilized iterationmethod of Eqs. (12)–(19) still protects us from obtain-ing complex RPA eigenvalues, but the condition numberof the Krylov-space RPA matrix approaches infinity, be-cause one or more of its eigenvalues collapse nearly tozero.In the standard method, the RPA matrix is calculatedby using the bare p-h basis states. In our method, weinstead start from the pivot vector of Eq. (20) and theArnoldi iteration then produces the rest of our basis vec-tors composing the Krylov subspace. This subspace is FIG. 10: Similar as Fig. 5 but for the 1 − strength functionscalculated for the numbers of HO shells N = 22, 26, 30, 34,and 38TABLE I: Spherical RPA and QRPA dimensions D as func-tions of the number of HO shells N .RPA QRPA N + − + + − +
10 70 195 261 390 1040 151020 205 555 766 2880 8180 1272025 273 734 1020 5538 15912 2508830 340 915 1271 9470 27420 43630 spanned by the eigenstates of the RPA matrix which hasan overlap with the pivot vector. Thus, in general, theArnoldi iterations can only be continued until this sub-space is exhausted in which case the condition numbergoes to infinity. However, with finite numerical preci-sion this maximum limit of Arnoldi iterations is furtherreduced.In a typical iteration, during the first few iterationsthe condition number of the Krylov-space RPA matrixfluctuates, then approaches a stable plateau, and finallysuddenly goes toward infinity. When that happens, theiteration must be stopped and one must backtrack tothe iteration where the condition number was still ac-ceptable. Therefore, the number of Arnoldi iterationscan depend on the size of the HO basis, and the resultspresented in this section correspond to the numbers ofiterations fixed according to this prescription.
VI. SCALING OF ITERATIVE SOLUTIONMETHOD
We illustrate the benefits of the iterative solution ofthe RPA or QRPA equations over the traditional methodby comparing how the numerical work increases in theiterative method as the HO basis is increased. As can be seen in Table I, the RPA dimensions D for doubly magic spherical nuclei increase almost lin-early with the number of oscillator shells N . This iseasy to understand, because in this case, only the num-ber of particle states increases and the number of holestates always stays constant. Therefore the time to solvethe full RPA eigenproblem in this case scales approxi-mately as N . In the spherical QRPA, the dimensionsscale roughly between N and N , and the full QRPAscales approximately as N or N . The physically inter-esting and computationally challenging calculations arefor deformed nuclei with pairing, and we should thereforecompare the N scaling of iterative and standard QRPAdiagonalizations.In the case of all symmetries of the mean field beingbroken, the QRPA dimension D is: D = 19 [( N + 1)( N + 2)( N + 3)] . (25)This dimension increases very steeply ( N ) as the num-ber of HO shells is increased. For N = 14, the QRPAdimension is D = 1849600, for example. The correspond-ing standard QRPA solution scales as N and is thusuntractable. Therefore, for the QRPA calculations in de-formed nuclei, we must truncate the single-particle space.The best method is to use the two-basis method [21], bywhich one solves the HFB equations in the basis gen-erated by the HF part of the HFB matrix, and truncatethe basis using a cutoff on the obtained pseudo-HF single-particle energies. But even then, the QRPA calculationsscale as the sixth power of the number of useful single-particle states, and are thus prohibitively difficult. FIG. 11: Times to calculate 100 Arnoldi iterations for thespherical QRPA method applied to
Sn as functions of N .Squares and circles show results for the 1 − and 2 + modes,respectively, and lines show cubic fits. In order to illustrate the scaling properties of the it-erative QRPA method, we calculated the correspondingmatrix-vector products with our developmental QRPAcode, where the pairing has been set to zero. The RPAfields ˜ h only depend on the normal RPA density matrixand the calculation of pairing part of the QRPA matrix-vector product is very fast due to a small number of thepairing coupling constants and the simple density depen-dence of typical pairing EDF. Therefore the time to cal-culate the pairing part is negligible. The only significantincrease of running time in spherical QRPA compared toRPA comes from the need to handle higher-dimensionalbasis vectors in the calculation of various overlaps andvector additions during the Arnoldi iteration. Therefore,as we keep all particle-particle and hole-hole RPA ampli-tudes in the calculation, but set them to zero, our tim-ing results accurately reflect the timing of the sphericalQRPA calculation.In Fig. 11, we show the scaling properties of the spher-ical QRPA calculation with iterative Arnoldi method. Itis clear that the scaling of our iterative method is as N ,that is, it is linear with respect to the QRPA dimen-sion D . Of course, the prefactor itself is linearly propor-tional to the number of Arnoldi iterations. However, asdiscussed in the previous section, the Arnoldi iterationmethod cannot in practice go full dimension before thegenerated Krylov-space matrices become singular. Aslong as we are satisfied with a few hundred iterationsat most, the iterative method gives us a vast speed im-provement. In the full RPA or QRPA diagonalization,the calculation and storage of a very large dense RPA orQRPA matrices also takes a considerable additional time– a step that the iterative method avoids completely.In addition to the moment-method based iteration,which is ideal for strength functions, the iterative methodcan also be modified to be suitable for different kinds ofother calculations. If we are interested in a number ofvery well converged lowest RPA eigenmodes, restartedArnoldi methods [22] can be used. These methods usemore iterations than basis states, i.e., after a maximumnumber d of basis vectors is generated, new approxima-tions for the wanted d ′ < d eigenmodes are calculated,and iteration is then continued to generate new improvedbasis states from d ′ + 1 to d again. The restarting canbe made as many times as needed to produce wantednumber of well converged lowest excitations.Methods such as Arnoldi or Lanczos produce conver-gence at the extreme ends of the excitation energy spec-trum. If eigenmodes away from the extremes are lookedfor, shift and invert methods [23, 24] can be used. Thesemethods allow iterative methods to be used to find RPAeigenmodes anywhere inside the RPA excitation spec-trum. VII. SUMMARY AND CONCLUSIONS
We have presented a method to calculate accurate RPAresponse functions by using the iterative Arnoldi diag-onalization related to the sum-rule conserving Lanczosmethod of Ref. [5]. We used strictly the same EDF forthe ground state calculation and RPA excitations. Wehave showed how the Arnoldi method must be stabilizedin order to apply it reliably to the RPA eigenvalue prob-lem. The resulting electromagnetic strength functionsare in good agreement with the standard RPA resultsand are obtained with numerical effort smaller by ordersof magnitude.Our method closely resembles the FAM of Nakatsukasa et al. [6, 7], except that our iterative method is differentand that we use the HO basis instead of the mesh incoordinate space. The FAM and our method both allowthe existing EDF mean-field codes to be used for thecalculation of the RPA or QRPA matrix-vector products.With minor modifications, mostly pertaining to the fullimplementation of the time-odd mean fields, these codescan easily be extended to RPA/QRPA. In particular, ourfuture implementation of the deformed QRPA solutionwill be based on the code HFODD [25].We also implemented the method to remove compo-nents of the spurious RPA modes from the calculatedstrength functions that keeps the physical excitations ex-actly orthogonal against the spurious excitations in anyfinite model space.The smaller numerical effort of the iterative Arnoldimethod, and the fact that in this method one does nothave to calculate and store the RPA or QRPA matrices,allows our method to be applied to the calculation ofelectromagnetic and beta decay strengths and strengthfunctions for deformed heavy nuclei. Work to extend ourformalism and codes to deformed superfluid nuclei is inprogress.
Acknowledgments
We are thankful to J. Terasaki for providing us withnumerical values of the strength functions calculated inRef. [17]. This work was supported by the Academy ofFinland and University of Jyv¨askyl¨a within the FIDIPROprogram and by the Polish Ministry of Science andHigher Education under Contract No. N N 202 328234. [1] S. Tretiak, C.M. Isborn, A.M.N. Niklasson, and M. Chal-lacombe, J. Chem. Phys. , 054111 (2009).[2] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert,M.S. Gordon, J.H. Jensen, S. Koseki, N. Matsunaga,K.A. Nguyen, S. Su, T.L. Windus, M. Dupuis, and J.A.Montgomery, Jr, J. Comput. Chem. , 632 (1992).[5] C.W. Johnson, G.F. Bertsch, and W.D. Hazelton, Com-put. Phys. Commun. , 155 (1999).[6] T. Nakatsukasa, T. Inakura, and K. Yabana, Phys. Rev.C , 024318 (2007). [7] T. Inakura, T. Nakatsukasa, and K. Yabana, Phys. Rev.C , 044301 (2009).[8] S. P´eru and H. Goutte, Phys. Rev. C , 044313 (2008).[9] K. Yoshida and Nguyen Van Giai, Phys. Rev. C ,064316 (2008).[10] J. Terasaki et al. , unpublished.[11] C. Losa et al. , unpublished.[12] P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer-Verlag, Berlin, 1980).[13] J.P. Blaizot and G. Ripka,
Quantum theory of finite sys-tems , MIT Press, Cambridge Mass., 1986.[14] B.K. Agrawal and S. Shlomo, Phys. Rev. C , 014308(2004).[15] Tapas Sil, S. Shlomo, B.K. Agrawal, and P.-G. Reinhard,Phys. Rev. C , 034316 (2006).[16] J. Terasaki, J. Engel, M. Bender, J. Dobaczewski, W.Nazarewicz, and M. Stoitsov, Phys. Rev. C , 034310(2005).[17] J. Terasaki and J. Engel, Phys. Rev. C , 044301 (2006).[18] Y. Saad, Iterative methods for sparse linear systems (2nd ed.) (SIAM, 2003).[19] B.G. Carlsson, J. Dobaczewski, J. Toivanen, and P.Vesel´y, to be published in Computer Physics Commu-nication.[20] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R.Schaeffer, Nucl. Phys.
A627 (1997) 710.[21] B. Gall, P. Bonche, J. Dobaczewski, H. Flocard, and P.-H. Heenen, Z. Phys.
A348 , 183 (1994).[22] R.B. Lehoucq and D.C. Sorensen, SIAM J. Matrix Anal.Appl. , 789 (1996).[23] Y. Saad, Numerical methods for large eigenvalue prob-lems, Algorithms and architectures for advanced scien-tific computing (Manchester University Press, Manch-ester, 1992).[24] Zhongxiao Jia and Yong Zhang, Computers Math. Ap-plic. , 1117 (2002).[25] J. Dobaczewski et al. , Comput. Phys. Commun. ,166 (1997); , 183 (1997); , 164 (2000); , 158(2004); , 214 (2005);180