Recursive method for computing matrix elements for two-body interactions

A recursive method for the efﬁcient computation of two-body matrix elements is presented. The method consists of a set of recursion relations for the computationally demanding radial integral and adds one more tool to the set of computational methods introduced by Horie and Sasaki [H. Horie and K. Sasaki, Prog. Theor. Phys. 25 , 475 (1961)]. The neutrinoless double-β decay will serve as the primary application and example, but the method is general and can be applied equally well to other kinds of nuclear structure calculations involving matrix elements of two-body interactions.


I. INTRODUCTION
Computation of matrix elements of two-body interactions is an inevitable task for many different kinds of nuclear structure calculations. For small numbers of matrix elements these computations can be done in a more or less straightforward way. In realistic situations, however, the amount of needed matrix elements is usually very large and this starts to pose serious demands for the needed computing resources. In this case more sophisticated methods for matrix element calculations are called for.
In the present work the single-particle wave functions are taken to be the eigenfunctions of the isotropic harmonic oscillator. The usual method to compute the two-body matrix elements in this case consists of the Moshinsky transformation and the associated transformation brackets. The remaining radial integral is reduced from two to just one dimension and can be calculated numerically at least. The Moshinsky transformation works very well for central interactions but for noncentral interactions it is not so convenient. This is a setback because for the neutrinoless double-β decay (0νββ), for example, one has to calculate matrix elements of a tensor force. The contribution of the tensor part to the total nuclear matrix element is small but for precise computations one has to take it into account. Because of these difficulties, in the present work the two-body matrix elements are calculated by using the Fourier-transform method of Horie and Sasaki, introduced in [1], followed by the application of the usual methods of spherical tensor algebra. The computationally most demanding part in this approach is the calculation of the involved radial integrals. To compute these integrals efficiently, a recursive method is introduced. In this method the radial integral is first expressed as a sum of simpler integrals. After this, a set of three-index recursion formulas for these integrals can be derived.
In this article we take the neutrinoless double-β decay to serve as the primary example and application of the presented recursive method, but the method itself is general and may be applied to other types of transition operators as well. The first rudimentary applications of the same method to neutrinoless double-β decay were done already quite some time ago [2].
This article is organized as follows. In Sec. II the Fourier transform method for computing two-body matrix elements for central and tensor interactions is reviewed. The radial integral is extracted from the two-body matrix element and a recursive method for its computation is introduced. In Sec. III the new method is applied to the case of 0νββ decay and a computational algorithm for the nuclear matrix element calculation is sketched. In Sec. III we also compare the efficiency of our new 0νββ code (based on the sketched algorithm) to our old 0νββ code and to the code presented in Ref. [3]. Sec. IV summarizes the present work.
The functions j λ are the spherical Bessel functions and the harmonic-oscillator wave functions involved are defined through the associated Laguerre polynomials L (α) n as with the normalization Next we use the method of Horie and Sasaki [1] to cast the radial integral (6) into the form where the J integral is defined as and the additional coefficients as where q i = (l + l + 2s i − λ i )/2 and ν = 2/b. The quantities a s (nl,n l ) originate from a Taylor expansion of a product of two associated Laguerre polynomials and have the following explicit expression: It is worth noting here that our definition (13) for the coefficients a s (nl,n l ) is different from that found from Ref. [1]. We see that calculation of the radial integral (6) boils down to the computation of the J integrals (11). We note in addition, that for many nuclear structure calculations the principal oscillator index n has a relatively small range of values. In that case the double summation in Eq. (10) does not actually contain very many terms. Our next task is to derive recursion formulas for the J integral. The standard way to find such formulas is by using an appropriate generating function. By applying the properties of the spherical Bessel functions and the harmonic-oscillator wave functions we get, after some algebra, the following set of recursion formulas: We notice that in the recursion formulas (14)-(16) the orbital angular momentum indices l and l change in a coherent manner. This is fine for the Fermi and GT type of interactions for which we need J integrals with l = l only. For the tensor interaction the situation is more complicated because we need also the integrals with l = l . However, one can deduce from the form of the matrix element (5) that we only need to consider cases with l = l or |l − l | = 2. Relations (14)-(16) can now be used to recursively generate all the possible J integrals. The needed seed points are 054308-2 of the form J λ (0l,0l ). The actual calculation procedure is highlighted in the next section, where the above derived recursion formulas are applied to the case of neutrinoless double-β decay.

III. APPLICATION: RADIAL INTEGRALS FOR 0νββ DECAY
As an example of the presented method we shall sketch a computational algorithm for calculating the radial integrals and two-particle matrix elements for the case of neutrinoless double-β decay with light-neutrino exchange. The 0νββ decay (Z,A) −→ (Z + 2,A) + 2e − is a nuclear process considered today as one of the best probes for beyond-standard-model physics [4,5]. Assuming the light-neutrino exchange as the dominant decay mechanism the 0νββ half-life can be expressed as where g A is the axial-vector coupling, G 0ν the phase-space factor of the final-state leptons, and M 0ν the nuclear matrix element which contains all the nuclear structure information related to the decay process. The quantity η 2 x contains the beyond-standard-model physics in the form of neutrinomixing matrix elements and neutrino mass eigenstates.
The nuclear matrix element has the following decomposition in terms of the Fermi (F), Gamow-Teller (GT), and tensor (T) contributions The individual parts M (0ν) K can be written as [6] M (0ν) K = J π ,k,J pp nn where k labels the different intermediate states for a given multipole J π . The operators O K inside the two-particle matrix element have the forms of Eqs. (1), (2), and (3), with the radial dependencies defined by the decay mechanism.
The simplest case in Eq. (19) is the Fermi matrix element, for which we have where we have the neutrino potential for the light neutrino exchange Here R A = 1.2A 1/3 fm is the nuclear radius, M i c 2 (M f c 2 ) is the ground-state mass energy of the initial (final) nucleus, and E k is the energy of the nuclear state k of the intermediate nucleus.
The Jastrow form of the short-range correlation function f SRC (r) is defined as whereâ,b, andĉ are constants which have particular values in different parametrizations of the short-range correlation effects. For the neutrino potential (21) we still need the dipole form factor where the vector mass is V = 843 MeV. For more details about the 0νββ decay mechanisms see, for example, [5].
The jj -coupled two-body matrix element in Eq. (19) can be expressed in terms of the LS-coupled elements which are given by Eq. (5). The matrix elements (5) can be computed as soon as we know the value for the radial integral with the two-body interaction given by Eq. (20). In the nuclear matrix elements (19) we have a summation over all the single-particle indices pp ,nn constrained by the chosen valence space and angular-momentum selection rules. This suggests that the amount of needed two-body matrix elements (5) and thus the amount of different radial integrals (and thus the amount of needed J integrals) is very large.
Our computation proceeds now as follows: We set up a recursion lattice with indices (n,(l,l ),n ) and use the relations (14)-(16) to generate all of the needed J integrals at once. To do this, we need to know the recursion upper bounds for the n,n and for the double index (l,l ). These have to be determined in a case-by-case basis. The only interaction-specific things we need to know are the seed-point values. For the Fermi-type interaction these are of the form J 0 (0l,0l). Let us compute these next.
Our J integral can be related to the Talmi integral I l (b) which is given by the expression Now the mentioned relation reads The remaining Talmi integrals can be computed by applying the method of Horie and Sasaki [1] where the integral over the radial coordinate ρ is performed analytically leaving a one-dimensional integral over the neutrino momentum q 054308-3 where p = √ 2bq and I (p) = N 2 0l {Î (1,1,2l,p) −Î (2ĉ,1 + α,2l,p) +Î (ĉ 2 ,1 + 2α,2l,p) −Î (2βĉ 2 ,1 + 2α,2l + 2,p) +Î (2ĉβ,1+α,2l +2,p)+Î (ĉ 2 β 2 ,1+2α,2l +4,p)}, with the abbreviations α = 2âb 2 and β = 2bb 2 . The remaining integrals to be computed are of the form where the coefficients B and m originate from Eq. (29). After making the abbreviations γ 2 = b 2 /2B and m/2 = s and expanding the associated Laguerre polynomial, we get for the final quantities to be computed expressions like where the coefficients c(s,k) are given by Eq. (26). The remaining integrals in Eq. (32) cannot be calculated analytically and the (time consuming) numerical integration becomes our only option. The potential function U (q) in Eq. (32) contains an energy dependence in the form of an energy denominator. This means that when we change the intermediate state [the nuclear matrix elements (19) contain a summation over all of them] we always have to recalculate the integrals (32). Numerically this might be too time consuming so let us try to isolate the energy dependence in the expression (32). By making a partial-fraction decomposition for U (q) we get where C i = C i ( k , V ) are real coefficients and k = E k − (M i c 2 + M f c 2 )/2 contains the energy dependence E k of the intermediate state. Inserting the decomposition (33) into Eq. (32) we see that the energy dependence is isolated to the last integral with coefficient C 5 . In the rest of the terms, energy dependence factors out from the integrals. Thus, these integrals need to be calculated only once. After that they can be stored for later use. The energy-dependent part becomes where the remaining integral can be reduced to the computation of the exponential integral Ei(x) as follows: As an alternative to the explicit sum formula (35), one can apply the following recursion relation where Values for the exponential integral can be calculated efficiently by applying the rational Chebyshev approximations [7,8]. FORTRAN subroutines based on these methods can be found, for example, from [9]. This concludes our example about the recursive computation of two-body matrix elements for neutrinoless double-β decay. The remaining Gamow-Teller and tensor matrix elements are more complicated to handle but the same procedure can be used to compute them also. Figure 1 illustrates the efficiency of our new 0νββ code which is based on the computational method sketched above. It shows two nuclear matrix element calculations for the nucleus 76 Ge with two different single-particle bases [The quasirandom-phase approximation (QRPA) was used to model the nuclear structure in the example calculations of this paper, but the application of the method is not limited to this model because the nuclear model used enters into the matrix element (19) only via the transition densities]. The first one is a toy calculation with only four oscillator basis states, the second calculation is considered as realistic with up to ten single-particle states included. In the toy calculation the speed-up factor when going from the old code to the new one is about 26. Of course a more interesting case is the realistic calculation for which we get an improved computing speed by a factor of about 35. In the calculations of Fig. 1 we have computed only the Fermi and Gamow-Teller parts of the total nuclear matrix element (18). The old code uses the Moshinsky transformation method to compute the two-body matrix elements. The involved radial integrals are computed numerically in two dimensions. Much CPU time is saved by applying a polynomial interpolation to compute radial integrals with different energy denominators. Figure 2 shows the performance of our new code when computing two-body matrix elements (TBMEs). For the toy calculation (with four oscillator basis states) the total number of needed TBMEs is 5946 (Fermi + GT). A massive increase in the computational burden occurs when we switch from a toy calculation to a realistic calculation (ten oscillator basis states). In the realistic case we need a total of 805 264 TBMEs as shown in Fig. 2. For the toy calculation we get a total CPU time of 41.51 s, thus an average time of 6.98 × 10 −3 s for a single TBME. In the realistic case the average CPU time per single TBME is greatly reduced, being about 2.62 × 10 −4 . This reduction can be traced back to the calculation of the J -integral initial values given by Eq. (25). To compute the initial values (for the Fermi-type elements) we must first numerically calculate the integrals related to the first eight partial fraction terms in Eq. (33). However, as already mentioned, we need to compute these integrals only once and after that we can store them for later use. The numerical calculation is a relatively time consuming process and gives more weight to the single TBME CPU time in the toy calculation where the total number of required TBMEs is still quite small. The CPU time consumption values of Fig. 2 can be compared with those given for similar calculations in [3]. In [3] the (antisymmetrized) TBMEs are computed by first applying the Horie and Sasaki method [1] and then integrating the remaining integrals numerically over the neutrino momentum. With this method, the following average CPU time consumption values for a single TBME are obtained: 6.69 × 10 −2 s (calculations with the nucleus 48 Ca) and 9.69 × 10 −2 s (calculations with the nucleus 82 Se). The CPU used for these computations was a single-core Intel Xeon X5560 with 2.8 GHz frequency. Comparing these figures with our average time consumption value of 2.62 × 10 −4 s clearly shows the computing power of the applied recursive algorithm: average CPU time for single TBME goes down by two orders of magnitude.

IV. CONCLUSIONS
In this work the Horie and Sasaki method for computing the radial integrals and two-body matrix elements using harmonic-oscillator wave functions was extended. The radial integral was first represented by a sum of simpler objects called J integrals. Then, the harmonic-oscillator generating function and the properties of the harmonic oscillator itself were applied to derive a set of very simple recursion relations for the J integral. All possible J integrals can now be generated recursively starting from the seed values of a form 054308-5 J (0l,0l ). The radial dependence of the two-body interaction was not specified in the derivation of the recursion formulas so they have a very general nature and can be applied to many different nuclear structure calculations where harmonicoscillator single-particle wave functions are used as basis states. Neutrinoless double-β decay was used as an example to sketch an algorithm to be used for the computation of the radial integrals and nuclear matrix elements. Finally, the efficiency of our new code, based on the sketched method, was compared with our older 0νββ code and an independent calculational procedure of Ref. [3]. A significant reduction in the computation time was witnessed when going from the mentioned computational schemes to the new code, thus confirming the value of the developed method.