Scaling of iterative solution method

We illustrate the benefits of the iterative solution of the RPA or QRPA equations over the traditional method by comparing how the numerical work increases in the iterative method as the HO basis is increased.

Table 1: Spherical RPA and QRPA dimensions $D$ as functions of the number of HO shells $N_{0}$.

		RPA			QRPA
$N_{0}$	$0^+$	$1^-$	$2^+$	$0^+$	$1^-$	$2^+$
$10$	$70$	$195$	$261$	$390$	$1040$	$1510$
$20$	$205$	$555$	$766$	$2880$	$8180$	$12720$
$25$	$273$	$734$	$1020$	$5538$	$15912$	$25088$
$30$	$340$	$915$	$1271$	$9470$	$27420$	$43630$

As can be seen in Table 1, the RPA dimensions $D$ for doubly magic spherical nuclei increase almost linearly with the number of oscillator shells $N_{0}$. This is easy to understand, because in this case, only the number of particle states increases and the number of hole states always stays constant. Therefore the time to solve the full RPA eigenproblem in this case scales approximately as $N_{0}^3$. In the spherical QRPA, the dimensions scale roughly between $N_{0}^2$ and $N_{0}^3$, and the full QRPA scales approximately as $N_{0}^6$ or $N_{0}^9$. The physically interesting and computationally challenging calculations are for deformed nuclei with pairing, and we should therefore compare the $N_{0}$ scaling of iterative and standard QRPA diagonalizations.

In the case of all symmetries of the mean field being broken, the QRPA dimension $D$ is:

$$D=\frac{1}{9}[(N_{0}+1)(N_{0}+2)(N_{0}+3)]^2\,.$$ (24)

This dimension increases very steeply ($N_0^{6}$) as the number of HO shells is increased. For $N_{0}=14$, the QRPA dimension is $D=1849600$, for example. The corresponding standard QRPA solution scales as $N_0^{18}$ and is thus untractable. Therefore, for the QRPA calculations in deformed nuclei, we must truncate the single-particle space. The best method is to use the two-basis method [21], by which one solves the HFB equations in the basis generated by the HF part of the HFB matrix, and truncate the basis using a cutoff on the obtained pseudo-HF single-particle energies. But even then, the QRPA calculations scale as the sixth power of the number of useful single-particle states, and are thus prohibitively difficult.

Figure 11: Times to calculate $100$ Arnoldi iterations for the spherical QRPA method applied to $^{132}{\rm Sn}$ as functions of $N_{0}$. 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 iterative QRPA method, we calculated the corresponding matrix-vector products with our developmental QRPA code, where the pairing has been set to zero. The RPA fields $\tilde h$ only depend on the normal RPA density matrix and the calculation of pairing part of the QRPA matrix-vector product is very fast due to a small number of the pairing coupling constants and the simple density dependence of typical pairing EDF. Therefore the time to calculate the pairing part is negligible. The only significant increase of running time in spherical QRPA compared to RPA comes from the need to handle higher-dimensional basis vectors in the calculation of various overlaps and vector additions during the Arnoldi iteration. Therefore, as we keep all particle-particle and hole-hole RPA amplitudes in the calculation, but set them to zero, our timing results accurately reflect the timing of the spherical QRPA calculation.

In Fig. 11, we show the scaling properties of the spherical QRPA calculation with iterative Arnoldi method. It is clear that the scaling of our iterative method is as $N_{0}^3$, that is, it is linear with respect to the QRPA dimension $D$. Of course, the prefactor itself is linearly proportional to the number of Arnoldi iterations. However, as discussed in the previous section, the Arnoldi iteration method cannot in practice go full dimension before the generated Krylov-space matrices become singular. As long as we are satisfied with a few hundred iterations at most, the iterative method gives us a vast speed improvement. In the full RPA or QRPA diagonalization, the calculation and storage of a very large dense RPA or QRPA 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 method can also be modified to be suitable for different kinds of other calculations. If we are interested in a number of very well converged lowest RPA eigenmodes, restarted Arnoldi methods [22] can be used. These methods use more iterations than basis states, i.e., after a maximum number $d$ of basis vectors is generated, new approximations for the wanted $d' < d$ eigenmodes are calculated, and iteration is then continued to generate new improved basis states from $d'+1$ to $d$ again. The restarting can be made as many times as needed to produce wanted number of well converged lowest excitations.

Methods such as Arnoldi or Lanczos produce convergence at the extreme ends of the excitation energy spectrum. If eigenmodes away from the extremes are looked for, shift and invert methods [23,24] can be used. These methods allow iterative methods to be used to find RPA eigenmodes anywhere inside the RPA excitation spectrum.