1 Introduction↩︎

Charge-changing nuclear response functions are important in nuclear-structure and weak-interaction phenomenology. They enter the description of \(\beta\) decay and double-\(\beta\) decay, and provide information about collective excitations and weak rates in nuclei across the chart of nuclides [1], [2]. For medium-mass and heavy open-shell nuclei, especially in deformed systems, the quasiparticle random-phase approximation (QRPA) built on nuclear energy-density functionals provides a practical microscopic framework for describing such excitations and transition strengths [1], [3].

In principle, QRPA strength functions can be obtained from the explicit QRPA matrix or from a full eigendecomposition of the linear-response problem. In practice, however, the matrix dimension grows rapidly with the size of the quasiparticle space, and explicit construction becomes prohibitively expensive for heavy and deformed nuclei. The finite amplitude method (FAM) addresses this difficulty by computing the QRPA response through induced fields, without assembling the full QRPA matrix explicitly [4], [5]. The FAM has become an effective tool for self-consistent linear-response calculations, including applications to superfluid and deformed nuclei and, in particular, to charge-changing transitions in axially deformed systems [6]–[8].

Despite their usefulness, conventional FAM-based calculations still proceed one frequency at a time. For each complex frequency on the chosen contour, one must solve a new linear system. This feature is well suited to pointwise evaluation of the response, but it can become costly when the goal is to reconstruct the strength function over a broad energy interval. That observation motivates Krylov-based approaches that capture spectral information more globally. These have appeared in iterative Arnoldi methods for strength functions, in spectral-density approximation techniques, and in recent polynomial and kernel-based approaches to QRPA response problems [9]–[13].

In this paper, we develop a symmetric Lanczos framework for approximating charge-changing QRPA strength functions. Starting from the FAM linear-response equation, we derive equivalent spectral formulations and a reduced eigenvalue problem involving the matrices \(MK\) and \(KM\). These reformulations make it possible to apply a symmetric Lanczos process with respect to the appropriate inner product and to approximate the smeared strength function over a wide frequency interval from a single Krylov run, rather than from separate solves at each frequency. Our goal is to retain the main efficiency advantages of FAM-style matrix-free calculations while exploiting the broader spectral information available through Lanczos projection.

The paper is organized as follows. Section 2 develops the QRPA formulations of the strength function, including the eigendecomposition-based and reduced-eigenproblem representations. Section 3 presents the frequency-by-frequency solvers and the symmetric Lanczos approximation. Section 4 reports numerical results for \(^{112}\)Sn and \(^{150}\)Nd, together with runtime comparisons.

1.0.0.1 Notation.

Throughout the paper, matrices and vectors may be real or complex, depending on the formulation under consideration. For a matrix \(A\in\mathbb{C}^{m\times n}\), we write \(A^\dagger\) for the Hermitian conjugate, \(A^T\) for the transpose, and \(A^*\) for the element-wise complex conjugate. The identity matrix is denoted by \(I\), with its size understood from context. When block matrices are used, their dimensions are likewise determined by the surrounding equations. We reserve \(\omega\) for the physical excitation energy and write \(\omega_\gamma=\omega+i\gamma\) when a finite Lorentzian smearing width \(\gamma>0\) is introduced. Unless stated otherwise, vectors are written as columns, and the Euclidean inner product is understood by default; when the analysis requires a different metric, such as the \(K\)-inner product used in the Lanczos formulation, we say so explicitly.

2 QRPA Formulation of the Strength Function↩︎

This section summarizes the QRPA equations and spectral reformulations that underlie the strength-function calculations developed later. Although the derivation passes through the linear response equation and the response function, the quantity of interest throughout is the charge-changing QRPA strength function. The iterative and Lanczos-based algorithms used to evaluate it are deferred to Section 3.

2.1 The FAM Equations and the Strength Function↩︎

The nuclear linear response, a function of frequency \(\omega\), is genrated by a weak external time-dependent field \(\hat{F}(t)\) of the form \[\label{eq:Ft} \hat{F}(t) = \eta \left(\hat{F}(\omega)e^{-i\omega t} + \hat{F}^\dagger(\omega)e^{i\omega t}\right) \,,\tag{1}\] where \(\eta\) is a small real parameter. Neglecting two-quasiparticle operators that do not contribute to the QRPA, the field can be written in the quasiparticle basis as \[\label{eq:Fw} \hat{F}(\omega) = \frac{1}{2}\sum_{\mu\nu}\left(F_{\mu\nu}^{20}(\omega)\hat{\alpha}_\mu^\dagger\hat{\alpha}_\nu^\dagger+F_{\mu\nu}^{02}(\omega)\hat{\alpha}_\mu\hat{\alpha}_\nu\right),\tag{2}\] where \(F_{\mu\nu}^{20}(\omega)\) and \(F_{\mu\nu}^{02}(\omega)\) are defined, e.g., in Ref.[1], and \(\hat{\alpha}_\mu^\dagger\) and \(\hat{\alpha}_\mu\) are quasiparticle creation and annihilation operators, respectively. For charge-changing applications of the finite-amplitude method in deformed nuclei, we refer in particular to [5]–[7].

In the QRPA, the time evolution of the quasiparticle operators under the external field \(\hat{F}(t)\) is governed by the time-dependent Hartree–Fock–Bogoliubov (TDHFB) equation: \[\label{eq:TDHFB95eq} i\partial_t\hat{\alpha}_\mu(t) = \left[\hat{H}(t)+\hat{F}(t), \hat{\alpha}_\mu(t)\right],\tag{3}\] where the quasiparticle operators acquire small-amplitude oscillations of the form \[\begin{align} \hat{\alpha}_\mu(t) &= \left(\hat{\alpha}_\mu+\delta\hat{\alpha}_\mu(t)\right)e^{iE_\mu t}, \tag{4} \\ \delta\hat{\alpha}_\mu(t) &= \eta\sum_\nu\hat{\alpha}_\nu^\dagger\left(X_{\nu\mu}(\omega)e^{-i\omega t}+Y_{\nu\mu}^*(\omega)e^{i\omega t}\right). \tag{5} \end{align}\] Here \(E_\mu\) is the \(\mu^\textrm{th}\) quasiparticle energy. \(X_{\nu\mu}(\omega)\) and \(Y_{\nu\mu}(\omega)\) are the forward and backward amplitudes, respectively. The TDHFB Hamiltonian \(\hat{H}(t)\) can be written as \[\label{eq:TDHFB} \hat{H}(t) = \hat{H}_0 + \delta\hat{H}(t),\tag{6}\] where \[\label{eq:HFB} \hat{H}_0 = \sum_\mu E_\mu\hat{\alpha}_\mu^\dagger\hat{\alpha}_\mu\tag{7}\] and the induced Hamiltonian is \[\label{eq:dHt} \delta\hat{H}(t) = \eta \left(\delta\hat{H}(\omega)e^{-i\omega t} + \delta\hat{H}^\dagger(\omega)e^{i\omega t}\right)\tag{8}\] with \[\label{eq:dHw} \delta\hat{H}(\omega) = \frac{1}{2}\sum_{\mu\nu}\left(\delta H_{\mu\nu}^{20}(\omega)\hat{\alpha}_\mu^\dagger\hat{\alpha}_\nu^\dagger + \delta H_{\mu\nu}^{02}(\omega)\hat{\alpha}_\mu\hat{\alpha}_\nu\right)\tag{9}\] representing a small-amplitude oscillation, and expressions for \(\delta H_{\mu\nu}^{20}(\omega)\) and \(\delta H_{\mu\nu}^{02}(\omega)\) given, e.g., in Ref.[1]. Substituting these definitions into the TDHFB equation 3 yields the FAM equations [5], [7]: \[\begin{align} \left(E_\mu+E_\nu-\omega\right)X_{\mu\nu}(\omega)+\delta H_{\mu\nu}^{20}(\omega) &= -F_{\mu\nu}^{20}(\omega), \tag{10} \\ \left(E_\mu+E_\nu+\omega\right)Y_{\mu\nu}(\omega)+\delta H_{\mu\nu}^{02}(\omega) &= -F_{\mu\nu}^{02}(\omega). \tag{11} \end{align}\] The quantity of interest is the strength function, defined by \[\label{eq:strength} \frac{dB(\omega;\hat{F})}{d\omega} = -\frac{1}{\pi}\text{Im}\left\{S(\omega;\hat{F})\right\},\tag{12}\] where \[\label{eq:response95amplitude} S(\omega;\hat{F}) = \sum_{\mu<\nu}\left(F_{\mu\nu}^{20}(\omega)^*X_{\mu\nu}(\omega) + F_{\mu\nu}^{02}(\omega)^*Y_{\mu\nu}(\omega)\right).\tag{13}\] We refer to \(S(\omega;\hat{F})\) in Eq. 13 as the QRPA response function for the external field \(\hat{F}\); it is the auxiliary quantity from which the strength function is obtained. In practice, we assume that the external weak field is frequency-independent, i.e., \(F_{\mu\nu}^{20}(\omega)=F_{\mu\nu}^{20}\) and \(F_{\mu\nu}^{02}(\omega)=F_{\mu\nu}^{02}\).

2.2 Obtaining the Strength Function by Solving Linear Equations↩︎

The FAM equations can also be written in a compact matrix form. This reformulation is useful for two reasons. First, it isolates the frequency dependence. Second, it makes it possible to apply standard linear-algebra tools directly to the QRPA strength-function problem. Expanding the induced Hamiltonian matrix elements \(\delta H_{\mu\nu}^{20}(\omega)\) and \(\delta H_{\mu\nu}^{02}(\omega)\) in terms of the QRPA matrices \(A\) and \(B\) [1] via \[\begin{align} \delta H_{\mu\nu}^{20}(\omega) &= -(E_\mu+E_\nu)X_{\mu\nu}(\omega) + \sum_{\mu'<\nu'}\left(A_{\mu\nu,\mu'\nu'}X_{\mu'\nu'}(\omega)+B_{\mu\nu,\mu'\nu'}Y_{\mu'\nu'}(\omega)\right), \tag{14} \\ \delta H_{\mu\nu}^{02}(\omega) &= -(E_\mu+E_\nu)Y_{\mu\nu}(\omega) + \sum_{\mu'<\nu'}\left(B_{\mu\nu,\mu'\nu'}^*X_{\mu'\nu'}(\omega)+A_{\mu\nu,\mu'\nu'}^*Y_{\mu'\nu'}(\omega)\right) \tag{15} \,, \end{align}\] we recast the FAM equations into the linear response matrix equation: \[\label{eq:lre} \left( \begin{bmatrix} A & B \\ B^* & A^* \end{bmatrix} - \omega \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix} \right) \begin{bmatrix} X(\omega) \\ Y(\omega) \end{bmatrix} = - \begin{bmatrix} F^{20} \\ F^{02} \end{bmatrix} \,,\tag{16}\] where \(X(\omega),Y(\omega),F^{20},F^{02}\in\mathbb{C}^n\) are the vectorized representations of the strict upper triangular parts \((\mu<\nu)\) of the respective matrices. One can evaluate the strength function on a frequency grid, \(\omega_i\), \(i=1,2,\ldots,n_\omega\), within an appropriate frequency range by solving 16 and using 12 and 13 for each \(\omega_i\). Note that, when \(\omega_i\) is near one of the eigenvalues of the matrix pencil \[\left(\begin{bmatrix} A & B \\ B^* & A^* \end{bmatrix}, \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix} \right) \,, \label{eq:matrixpen}\tag{17}\] 16 becomes ill-conditioned or singular. To avoid these singularities in practical calculations, we replace the real frequency by the complex value \(\omega \to \omega_\gamma := \omega + i\gamma\) with a smearing width \(\gamma > 0\). The strength function then becomes, \[\label{eq:strength2} \frac{dB(\omega;\hat{F})}{d\omega} = \lim_{\gamma\to0^+}-\frac{1}{\pi}\text{Im}\left\{S(\omega_\gamma;\hat{F})\right\} = \lim_{\gamma\to0^+}-\frac{1}{\pi}\text{Im}\left\{\begin{bmatrix} F^{20} \\ F^{02} \end{bmatrix}^\dagger \begin{bmatrix} X(\omega_\gamma) \\ Y(\omega_\gamma) \end{bmatrix}\right\} \,.\tag{18}\]

2.3 Representing the Strength Function through Eigendecomposition↩︎

An alternative way to compute the strength function is to first perform an eigendecomposition of the QRPA matrix pencil 17 , and then rewrite the strength function in terms of the QRPA eigenmodes. This approach allows us to evaluate the strength function at any frequency \(\omega\) once the eigendecomposition of 17 is obtained. The special matrix structure of 17 yields a corresponding structure in its eigendecomposition, as described below.

If the QRPA matrices satisfy \(A^\dagger=A, B^T=B\) and the matrix \(\begin{bmatrix} A & B \\ B^* & A^* \end{bmatrix}\) is positive-definite, there exists an eigendecomposition whose eigenvalues are real and come in pairs. We refer the reader to [12] and [14] for a proof.

By Theorem 1, the matrix in 16 is invertible for any \(\omega_\gamma\) away from the real poles, and in particular for every \(\gamma>0\): \[\label{eq:lre95inv} \begin{align} &\left(\begin{bmatrix} A & B \\ B^* & A^* \end{bmatrix} - \omega_\gamma \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix} \right)^{-1} = \begin{bmatrix} X & Y^* \\ Y & X^* \end{bmatrix} \begin{bmatrix} (\Omega-\omega_\gamma I)^{-1} & 0 \\ 0 & (-\Omega-\omega_\gamma I)^{-1} \end{bmatrix} \begin{bmatrix} X & Y^* \\ Y & X^* \end{bmatrix}^{-1} \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix}\,. \end{align}\tag{19}\] For compactness in the following derivation, we define \[\mathcal{U}:= \begin{bmatrix} X & Y^* \\ Y & X^* \end{bmatrix}, \qquad \mathcal{F}:= \begin{bmatrix} F^{20} \\ F^{02} \end{bmatrix}\,.\] Using the vectorized notation introduced above and substituting the inverse 19 into Eq. 13 yields \[\label{eq:strength3} \begin{align} S(\omega_\gamma;\hat{F}) &= \mathcal{F}^\dagger \begin{bmatrix} X(\omega_\gamma) \\ Y(\omega_\gamma) \end{bmatrix} \\ &= -\mathcal{F}^\dagger \left(\begin{bmatrix} A & B \\ B^* & A^* \end{bmatrix} - \omega_\gamma \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix} \right)^{-1} \mathcal{F} \\ &= -\left(\mathcal{U}^\dagger\mathcal{F}\right)^\dagger \begin{bmatrix} (\Omega-\omega_\gamma I)^{-1} & 0 \\ 0 & (-\Omega-\omega_\gamma I)^{-1} \end{bmatrix} \left( \mathcal{U}^{-1} \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix} \mathcal{F} \right) \\ &= -\left(\mathcal{U}^\dagger\mathcal{F}\right)^\dagger \begin{bmatrix} (\Omega-\omega_\gamma I)^{-1} & 0 \\ 0 & (\Omega+\omega_\gamma I)^{-1} \end{bmatrix} \left(\mathcal{U}^\dagger\mathcal{F}\right)\,. \end{align}\tag{20}\] This implies that \[\label{eq:Im95strength} \begin{align} -\frac{1}{\pi}\text{Im}\left\{S(\omega_\gamma;\hat{F})\right\} = \left( \begin{bmatrix} X & Y^* \\ Y & X^* \end{bmatrix}^\dagger \begin{bmatrix} F^{20} \\ F^{02} \end{bmatrix} \right)^\dagger \begin{bmatrix} \frac{\gamma/\pi}{(\omega-\Omega)^2+\gamma^2} & 0 \\ 0 & \frac{-\gamma/\pi}{(\omega+\Omega)^2+\gamma^2} \end{bmatrix} \left( \begin{bmatrix} X & Y^* \\ Y & X^* \end{bmatrix}^\dagger \begin{bmatrix} F^{20} \\ F^{02} \end{bmatrix} \right)\,, \end{align}\tag{21}\] where we have used the corresponding identity for the diagonal matrix appearing above. It follows, after taking the limit \(\gamma\to0^+\), that the strength function 18 can be expressed as a sum of Dirac delta functions located at the QRPA poles, \[\label{eq:response3} \frac{dB(\omega;\hat{F})}{d\omega} = \sum_{i=1}^n |T_i(\hat{F})|^2\delta(\omega-\Omega_i) - \sum_{i=1}^n |\tilde{T}_i(\hat{F})|^2\delta(\omega+\Omega_i),\tag{22}\] where \(T_i(\hat{F})\) and \(\tilde{T}_i(\hat{F})\) are defined as \[\label{eq:Ti} \begin{bmatrix} T(\hat{F}) \\ \tilde{T}(\hat{F}) \end{bmatrix} = \begin{bmatrix} X & Y^* \\ Y & X^* \end{bmatrix}^\dagger \begin{bmatrix} F^{20} \\ F^{02} \end{bmatrix}.\tag{23}\] In the eigenmode representation of Eq. 22 , the projected quantities \(T_i(\hat{F})\) and \(\tilde{T}_i(\hat{F})\) are often referred to as QRPA transition amplitudes to the positive- and negative-frequency branches, respectively.

In Figure 1, we plot example QRPA transition amplitudes and the corresponding strength function. The magnitudes of \(T_i(\hat{F})\) and \(\tilde{T}_i(\hat{F})\) do not necessarily match, so the resulting strength function is not, in general, antisymmetric about the origin.

2.4 Reformulation as a Reduced Eigenvalue Problem↩︎

In practice, the QRPA matrices \(A\) and \(B\) and the external weak fields \(F^{20},F^{02}\) are often real. In that case, the dimension of the eigenvalue problem can be reduced by two, by using the products \(MK\) and \(KM\), where \(M=A+B\) and \(K=A-B\). Under the positivity assumptions inherited from Theorem 1, both \(M\) and \(K\) are symmetric positive-definite. Related reductions also appear in [1], [11]. We summarize this construction here for completeness.

The real counterpart of the eigendecomposition ?? is \[\label{eq:eigendecomp95real} \begin{bmatrix} A & B \\ B & A \end{bmatrix} \begin{bmatrix} X & Y \\ Y & X \end{bmatrix} = \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix} \begin{bmatrix} X & Y \\ Y & X \end{bmatrix} \begin{bmatrix} +\Omega & 0 \\ 0 & -\Omega \end{bmatrix},\tag{24}\] with \(X,Y\in\mathbb{R}^{n\times n}\). In addition, the normalization constraints ?? ,?? indicate that \[\label{eq:normalization4} X^TX-Y^TY = I, \quad X^TY-Y^TX=0, \quad XX^T-YY^T=I, \quad XY^T-YX^T = 0.\tag{25}\] These conditions imply that \[\label{eq:normalization5} (X+Y)^{-1} = (X-Y)^T \quad\text{and}\quad (X-Y)^{-1} = (X+Y)^T.\tag{26}\] Note that the eigendecomposition 24 is equivalent to \[\label{eq:eigendecomp95real2} \begin{bmatrix} A & B \\ -B & -A \end{bmatrix} \begin{bmatrix} X & Y \\ Y & X \end{bmatrix} = \begin{bmatrix} X & Y \\ Y & X \end{bmatrix} \begin{bmatrix} +\Omega & 0 \\ 0 & -\Omega \end{bmatrix},\tag{27}\] and the matrix \[\label{eq:Casida} \begin{bmatrix} A & B \\ -B & -A \end{bmatrix}\tag{28}\] is often referred to as the Casida matrix.

One can rewrite the eigenvalue problem for the Casida matrix, \[\label{eq:evp} \begin{bmatrix} A & B \\ -B & -A \end{bmatrix} \begin{bmatrix} x_i \\ y_i \end{bmatrix} = \Omega_i \begin{bmatrix} x_i \\ y_i \end{bmatrix},\tag{29}\] by employing a unitary similarity transformation, \[\label{eq:J} J = \frac{1}{\sqrt{2}} \begin{bmatrix} I & I \\ I & -I \end{bmatrix}\tag{30}\] to obtain \[\label{eq:tran95evp} \begin{bmatrix} 0 & K \\ M & 0 \end{bmatrix} \begin{bmatrix} x_i+y_i \\ x_i-y_i \end{bmatrix} = \Omega_i \begin{bmatrix} x_i+y_i \\ x_i-y_i \end{bmatrix},\tag{31}\] where \(K:=A-B\) and \(M:=A+B\). Then we have \[\label{eq:evp95KM} \begin{align} K(x_i-y_i) = \Omega_i (x_i+y_i), \\ M(x_i+y_i) = \Omega_i (x_i-y_i), \end{align}\tag{32}\] and \[\label{eq:evp95KM2} \begin{align} MK(x_i-y_i) = \Omega_i^2 (x_i-y_i), \\ KM(x_i+y_i) = \Omega_i^2 (x_i+y_i). \end{align}\tag{33}\] In other words, \[\label{eq:MK95decomp} M = (X-Y)\Omega(X-Y)^T \quad\text{and}\quad K = (X+Y)\Omega(X+Y)^T,\tag{34}\] and \[\begin{align} MK &= (X-Y)\Omega^2(X-Y)^{-1}, \tag{35} \\ KM &= (X+Y)\Omega^2(X+Y)^{-1}, \tag{36} \end{align}\]

These identities show that the eigendecomposition 27 of dimension \(2n\) can be replaced by an eigenvalue problem of dimension \(n\) for either the matrix \(MK\) or the matrix \(KM\). In addition, once we have computed an eigenvector \(x_i-y_i\) of \(MK\), we can recover \(x_i+y_i\) from \[\label{eq:xipyi} x_i+y_i = \frac{1}{\Omega_i}K(x_i-y_i) \,,\tag{37}\] without having to find an eigendecomposition of \(KM\). Once the vectors \(x_i+y_i\) and \(x_i-y_i\) are known, the original QRPA eigenvectors \(x_i\) and \(y_i\) follow immediately.

The reduced formulation also gives a direct spectral representation of the strength function. To see this, we insert the unitary similarity matrix \(J\) 30 into \[\begin{bmatrix} X & Y \\ Y & X \end{bmatrix}^T \begin{bmatrix} F^{20} \\ F^{02} \end{bmatrix}\] as \[\begin{bmatrix} X & Y \\ Y & X \end{bmatrix}^T JJ \begin{bmatrix} F^{20} \\ F^{02} \end{bmatrix}.\] This yields \[\label{eq:Ri} \begin{bmatrix} R(\hat{F}) \\ \tilde{R}(\hat{F}) \end{bmatrix} := \frac{1}{2} \begin{bmatrix} (X+Y)^T(F^{20}+F^{02}) + (X-Y)^T(F^{20}-F^{02}) \\ (X+Y)^T(F^{20}+F^{02}) - (X-Y)^T(F^{20}-F^{02}) \end{bmatrix},\tag{38}\] so that the strength function can be written as \[\label{eq:strength95Lanczos} \frac{dB(\omega;\hat{F})}{d\omega} = \sum_{i=1}^n |R_i(\hat{F})|^2\delta(\omega-\Omega_i) - \sum_{i=1}^n |\tilde{R}_i(\hat{F})|^2\delta(\omega+\Omega_i).\tag{39}\] This representation will serve as the starting point for the Lanczos approximation developed in the next section.

The reduced formulation above is the key structural ingredient for the numerical methods developed in the next section. In particular, it provides the reduced operator on which the symmetric Lanczos approximation is built.

3 Computational Methods for the Strength Function↩︎

This section describes the numerical strategies used to evaluate the QRPA strength function in practice. We first summarize two frequency-by-frequency approaches based on the linear response equation and then derive the symmetric Lanczos approximation based on the reduced eigenvalue problem.

3.1 Frequency-by-Frequency Solvers↩︎

3.1.1 Finite Amplitude Method (FAM)↩︎

A standard approach to solving Eqs. 10 and 11 is the iterative finite amplitude method (FAM) [5], [7]. For each fixed frequency \(\omega\), FAM treats the response equations as a fixed-point problem for the amplitudes \(X_{\mu\nu}(\omega)\) and \(Y_{\mu\nu}(\omega)\). A basic iteration proceeds as follows:

1. Use the current amplitudes \(X_{\mu\nu}(\omega)\) and \(Y_{\mu\nu}(\omega)\) to compute \(\delta H_{\mu\nu}^{20}(\omega)\) and \(\delta H_{\mu\nu}^{02}(\omega)\).

2. Add the external field contributions, \[\begin{align} \delta H_{\mu\nu}^{20}(\omega) &\leftarrow \delta H_{\mu\nu}^{20}(\omega) + F_{\mu\nu}^{20}, \\ \delta H_{\mu\nu}^{02}(\omega) &\leftarrow \delta H_{\mu\nu}^{02}(\omega) + F_{\mu\nu}^{02}. \end{align}\]

3. Update the amplitudes according to \[\label{eq:FAM95updates} \begin{align} X_{\mu\nu}(\omega) &\leftarrow -\left(E_\mu+E_\nu-\omega\right)^{-1}\delta H_{\mu\nu}^{20}(\omega), \\ Y_{\mu\nu}(\omega) &\leftarrow -\left(E_\mu+E_\nu+\omega\right)^{-1}\delta H_{\mu\nu}^{02}(\omega). \end{align}\tag{40}\]

The FAM is effective because it avoids constructing the full QRPA matrix explicitly. However, the plain fixed-point iteration is not guaranteed to converge for all frequencies: as with any stationary fixed-point scheme, convergence requires the associated iteration map to be sufficiently contractive. In practice, convergence can become slow or fail near poorly conditioned response points, especially close to QRPA poles. For this reason, practical FAM calculations commonly use acceleration schemes such as Broyden mixing [15]. Even with such acceleration, the full iteration must be repeated at every frequency point in the target interval.

3.1.2 GMRES Solver↩︎

Because Eq. 16 is a linear system, Krylov-subspace methods such as the generalized minimal residual method (GMRES) [16] provide a natural alternative to fixed-point FAM iterations. For a given complex frequency \(\omega_\gamma\), GMRES seeks the approximation in the associated Krylov subspace that minimizes the residual norm. In our calculations, this typically leads to significantly fewer iterations than plain FAM updates.

To build the Krylov basis, the matrices \(A\) and \(B\) need not be formed explicitly. Each GMRES iteration only requires the matrix–vector product \[\label{eq:lre95matvec} \begin{bmatrix} X(\omega_\gamma) \\ Y(\omega_\gamma) \end{bmatrix} \to \left( \begin{bmatrix} A & B \\ B & A \end{bmatrix} - \omega_\gamma \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix} \right) \begin{bmatrix} X(\omega_\gamma) \\ Y(\omega_\gamma) \end{bmatrix}.\tag{41}\] Evaluating this product is equivalent to computing the left-hand side of the FAM equations 10 and 11 . As with the FAM, however, a separate GMRES solve is required for every frequency value of interest. In Section 4.1.2, we describe how this matvec is assembled in practice and present the preconditioner used to accelerate GMRES convergence.

3.2 Lanczos Approximation of the Strength Function↩︎

We now use the reduced formulation to derive a Lanczos approximation for the strength function. Instead of fully diagonalizing 35 (or 36 ), which can still be expensive even after reduction, we approximate the exact spectral representation 39 by means of the Lanczos method [17]. From a numerical perspective, the appeal of this approach is that one Lanczos run can capture the dominant spectral information needed over a broad frequency interval. Closely related Krylov approaches have also been used for linear-response strength calculations in nuclear systems; see, for example, [9].

Although the matrix \(MK\) is not symmetric in the Euclidean inner product, it is self-adjoint with respect to the \(K\)-inner product, i.e., \[\label{eq:Kinner} \langle x, MKy\rangle_K = x^TKMKy = (MKx)^TKy = \langle MKx, y\rangle_K,\tag{42}\] for any \(x,y\in\mathbb{R}^n\). This observation allows us to apply a symmetric Lanczos process in the \(K\)-inner product rather than working with a nonsymmetric algorithm in the Euclidean inner product. An \(m\)-step Lanczos process generates \[\label{eq:Lanczos} MKQ_m = Q_mT_m+f_me_m^T\tag{43}\] where \[\label{eq:Qm95Lanczos} Q_m^TKQ_m = I \quad\text{and}\quad Q_m^TKf_m = 0.\tag{44}\] The projected matrix \(T_m\in\mathbb{R}^{m\times m}\) is tridiagonal and has the eigendecomposition \[T_mY_m = Y_m\Theta_m.\] The eigenvalues \(\{\theta_i\}_{i=1}^m\) of the tridiagonal matrix \(T_m\) approximate the squared excitation energies \(\{\Omega_i^2\}_{i=1}^n\). The corresponding excitation energies are therefore approximated as \[\label{eq:eigs95appr} \Omega_i\approx\sqrt{\theta_i}.\tag{45}\] The eigenvectors \((X-Y)\in\mathbb{R}^{n\times n}\) of \(MK\) are approximated by \(Z_m\in\mathbb{R}^{n\times m}\), defined as \[\label{eq:MK95Ritzvecs} Z_m:=Q_mY_m.\tag{46}\] These Ritz vectors, however, are normalized through the relation \[\label{eq:normal1} Z_m^TKZ_m=I,\tag{47}\] which is different from the normalization \[\label{eq:reduced95normalization} (X-Y)^TK(X-Y)=\Omega\tag{48}\] that holds for the eigenvector \(X-Y\) of \(MK\) as in 34 .

To relate this normalization to the exact QRPA normalization, suppose for a moment that the reduced eigenproblem were solved exactly (for example, through an \(n\)-step Lanczos process) and that we had obtained \[\label{eq:MK95evp} MKZ=Z\Theta\tag{49}\] where \[\label{eq:normal2} Z^TKZ=I \quad and \quad \Theta = \Omega^2.\tag{50}\] Comparing the normalization 50 with the normalization 48 , we immediately obtain \[\label{eq:X-minus-Y} X-Y=Z\Omega^{1/2} = Z\Theta^{1/4}.\tag{51}\] Moreover, since \[X+Y=K(X-Y)\Omega^{-1},\] we have \[\label{eq:X-plus-Y} X+Y = KZ\Omega^{-1/2} = KZ\Theta^{-1/4}.\tag{52}\]

Thus, once the Ritz pairs \((\Theta_m,Z_m)\) have been computed from an \(m\)-step Lanczos run, the strength function is approximated by \[\label{eq:lanczos95strength95approx} \frac{dB(\omega;\hat{F})}{d\omega} \approx \sum_{i=1}^m |S_i(\hat{F})|^2\delta(\omega-\sqrt{\theta_i}) - \sum_{i=1}^m |\tilde{S}_i(\hat{F})|^2\delta(\omega+\sqrt{\theta_i}),\tag{53}\] where the vectors \(S(\hat{F})\) and \(\tilde{S}(\hat{F})\) have components \(S_i(\hat{F})\) and \(\tilde{S}_i(\hat{F})\), respectively, and are defined as \[\label{eq:Si} \begin{bmatrix} S(\hat{F}) \\ \tilde{S}(\hat{F}) \end{bmatrix} = \frac{1}{2} \begin{bmatrix} (KZ_m\Theta_m^{-1/4})^T(F^{20}+F^{02}) + (Z_m\Theta_m^{1/4})^T(F^{20}-F^{02}) \\ (KZ_m\Theta_m^{-1/4})^T(F^{20}+F^{02}) - (Z_m\Theta_m^{1/4})^T(F^{20}-F^{02}) \end{bmatrix}.\tag{54}\]

One might worry that replacing the full \(n\)-term spectral sum by only \(m\) Ritz contributions is too severe an approximation, or that the apparent agreement is driven primarily by replacing each Dirac delta peak by a smooth kernel such as a Lorentzian. The smoothing is certainly important for comparing strength profiles on a finite-frequency grid, but the quality of the approximation also depends on whether the projected transition amplitudes \(S_i(\hat{F})\) and \(\tilde{S}_i(\hat{F})\) capture the relevant QRPA strength carried by the external field. In the numerical section of this paper, we therefore assess the smoothed strength functions through pointwise agreement, KL divergence, and runtime comparisons.

3.2.1 Comparison with Frequency-by-Frequency Solvers↩︎

The conceptual difference between the Lanczos approach and the FAM/GMRES solvers is worth emphasizing. The Lanczos method approximates the spectral representation of the strength function, whereas the FAM and GMRES approximate the response function at one prescribed frequency. Note that the strength function defined as in 12 , \[\frac{dB(\omega;\hat{F})}{d\omega} = -\frac{1}{\pi}\text{Im}\left\{S(\omega;\hat{F})\right\},\] where \[S(\omega;\hat{F}) = \sum_{\mu<\nu}\left(F_{\mu\nu}^{20}(\omega)^*X_{\mu\nu}(\omega) + F_{\mu\nu}^{02}(\omega)^*Y_{\mu\nu}(\omega)\right),\] requires the amplitudes \(X_{\mu\nu}(\omega)\) and \(Y_{\mu\nu}(\omega)\) for a fixed frequency \(\omega\). Consequently, both the FAM and GMRES must be run separately for each frequency at which the strength function is requested.

On the other hand, the strength function written in terms of the eigendecomposition as in 22 , \[\frac{dB(\omega;\hat{F})}{d\omega} = \sum_{i=1}^n |T_i(\hat{F})|^2\delta(\omega-\Omega_i) - \sum_{i=1}^n |\tilde{T}_i(\hat{F})|^2\delta(\omega+\Omega_i),\] where \[\begin{bmatrix} T(\hat{F}) \\ \tilde{T}(\hat{F}) \end{bmatrix} = \begin{bmatrix} X & Y^* \\ Y & X^* \end{bmatrix}^\dagger \begin{bmatrix} F^{20} \\ F^{02} \end{bmatrix},\] shows that \(\omega\) is decoupled from the transition-amplitude calculation. Because the Lanczos procedure approximates the eigenvalues \(\Omega_i\) together with the QRPA transition amplitudes \(T(\hat{F})\) and \(\tilde{T}(\hat{F})\), a single Lanczos run can be used to approximate the strength function over an entire frequency interval. The FAM and GMRES are therefore attractive when only a small number of frequencies are needed, whereas the Lanczos method becomes advantageous when one seeks a broad spectral scan or a full response profile.

3.2.2 Lanczos Initialization↩︎

The symmetric Lanczos process is applied to the reduced matrix \(MK\) with the \(K\)-inner product. We initialize the recursion with the normalized vector associated with the external field combination \[f_{+} := F^{20}+F^{02}, \qquad q_1 = \frac{f_{+}}{\|f_{+}\|_{K}}, \qquad \|v\|_{K}:=(v^TKv)^{1/2}.\] This choice is natural in the reduced formulation because the combination \(X+Y\), which is reconstructed from \(KZ_m\Theta_m^{-1/4}\), couples directly to \(f_{+}\) in Eq. 54 . It also respects the symmetry structure of the axially-deformed calculation. In our implementation, the quasiparticle basis and QRPA matrices inherit selection rules associated with the intrinsic projection quantum number, so \(A\), \(B\), and the reduced operator \(MK\) decompose into invariant symmetry sectors. Because \(F^{20}\) and \(F^{02}\) are generated from an external transition operator with a specified intrinsic component, \(f_{+}\) already lies in the sector relevant to the desired response.

This symmetry-adapted support is important in practice. An unrestricted random starting vector generally contains components in other decoupled sectors and therefore samples spectral information unrelated to the chosen strength function. A random vector restricted to the same nonzero pattern and sign structure as \(f_{+}\) removes most of this mismatch and can reproduce the overall profile well, but it still lacks the physical matrix-element magnitudes of the external field. Starting from \(f_{+}\) therefore both selects the correct symmetry sector and preserves the field-dependent weights needed for the detailed strength profile. We also retain the projections associated with \(f_-:=F^{20}-F^{02}\) during the Lanczos calculation, since Eq. 54 shows that both \(f_{+}\) and \(f_-\) enter the reconstructed transition amplitudes.

4 Numerical Results↩︎

In this section we assess the performance of the symmetric Lanczos method for approximating QRPA strength functions. We first discuss several practical considerations, including the regularization of the Dirac-delta distribution and the construction of matrix–vector products (matvecs) for GMRES and Lanczos. We then present calculations for two realistic nuclei: the medium-mass nucleus \(^{112}\)Sn and the heavier rare-earth nucleus \(^{150}\)Nd. The numerical discussion is organized in stages: first, we compare GMRES with the conventional iterative FAM (IFAM) to establish GMRES as an efficient frequency-by-frequency reference calculation; next, we compare the Lanczos approximation directly with this GMRES reference using strength profiles, pointwise errors, KL divergences, and runtimes.

All iterative solvers used in this section were implemented through the Python package pynfam [8]. This package provides a workflow around the Fortran codes hfbtho and pnfam for charge-changing HFB+QRPA calculations [7], [18]. For each nucleus, hfbtho is first used to obtain the HFB ground state through an axial deformation scan, after which the lowest-energy solution is taken as the reference state for the subsequent pnfam calculation. The latter computes the linear response for the chosen external field on the prescribed energy contour, from which the corresponding strength function is constructed.

Unless otherwise noted, we use common numerical settings across the two nuclei and vary only the nucleus-dependent inputs; both \(^{112}\)Sn and \(^{150}\)Nd are computed with 16 harmonic-oscillator shells. The underlying mean field is generated with SkM\(^*\), a standard Skyrme energy-density functional widely used in self-consistent mean-field and QRPA calculations [19]. In the pnFAM calculations we do not explicitly include the \(J^2\) terms; that is, the spin-current contributions associated with that sector of the Skyrme functional are not turned on by hand in our setup. Detailed discussions of the individual nuclei and their corresponding response functions are deferred to the subsections below.

4.1 Practical Considerations↩︎

4.1.1 Regularization↩︎

Because no conventional function possesses the properties of the Dirac-delta function, it is often defined through limits or the theory of distributions. A common way to define the delta function is as the limit of a Gaussian function \[\label{eq:Gaussian} g_\sigma(t) = \frac{1}{(2\pi\sigma^2)^{1/2}} e^{-\frac{t^2}{2\sigma^2}},\tag{55}\] or as the limit of a Lorentzian function \[\label{eq:Lorentzian} L_\gamma(t) = \frac{\gamma/\pi}{t^2+\gamma^2}.\tag{56}\] As the width parameters \(\sigma\) and \(\gamma\) of the Gaussian and Lorentzian, respectively, approach zero, the corresponding functions become infinitely tall and infinitely narrow, resembling the Dirac-delta function.

For density-of-states (DOS) problems, replacing the delta function with such an approximate function is known as regularizing the spectral density, and it is shown in [10] that the width parameter controls the resolution of the DOS curves. In this paper, we adapt this regularization strategy for our strength function calculations, replacing the delta functions with Lorentzian functions.

4.1.2 Matvecs for GMRES and Lanczos↩︎

The equivalence between the FAM equations 10 ,11 and the linear response equation 16 indicates that we can obtain the matvec \[\label{eq:matvec} \begin{bmatrix} X(\omega_\gamma) \\ Y(\omega_\gamma) \end{bmatrix} \mapsto \left( \begin{bmatrix} A & B \\ B & A \end{bmatrix} - \omega_\gamma \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix} \right) \begin{bmatrix} X(\omega_\gamma) \\ Y(\omega_\gamma) \end{bmatrix}\tag{57}\] by performing the following steps for each \(\mu<\nu\).

1. Use the amplitudes \(X_{\mu\nu}(\omega_\gamma)\) and \(Y_{\mu\nu}(\omega_\gamma)\) to compute \(\delta H_{\mu\nu}^{20}(\omega_\gamma)\) and \(\delta H_{\mu\nu}^{02}(\omega_\gamma)\).

2. Add \(\delta H_{\mu\nu}^{20}(\omega_\gamma)\) to \(\left(E_\mu+E_\nu-\omega_\gamma\right)X_{\mu\nu}(\omega_\gamma)\) and add \(\delta H_{\mu\nu}^{02}(\omega_\gamma)\) to \(\left(E_\mu+E_\nu+\omega_\gamma\right)Y_{\mu\nu}(\omega_\gamma)\).

4.1.2.1 GMRES.

For implementing GMRES, we use the diagonal preconditioner \[\label{eq:precond} P = \begin{bmatrix} \left(E_\mu+E_\nu-\omega_\gamma\right)^{-1} & 0 \\ 0 & \left(E_\mu+E_\nu+\omega_\gamma\right)^{-1} \end{bmatrix}.\tag{58}\] The corresponding preconditioned matvec is \[\label{eq:gmres95matvec} \begin{bmatrix} X(\omega_\gamma) \\ Y(\omega_\gamma) \end{bmatrix} \mapsto \begin{bmatrix} \left(E_\mu+E_\nu-\omega_\gamma\right)^{-1} & 0 \\ 0 & \left(E_\mu+E_\nu+\omega_\gamma\right)^{-1} \end{bmatrix} \left(\begin{bmatrix} A & B \\ B & A \end{bmatrix} - \omega_\gamma \begin{bmatrix} I & 0 \\ 0 & -I \end{bmatrix}\right) \begin{bmatrix} X(\omega_\gamma) \\ Y(\omega_\gamma) \end{bmatrix},\tag{59}\] which can be obtained by performing the following steps:

1. Use the amplitudes \(X_{\mu\nu}(\omega_\gamma)\) and \(Y_{\mu\nu}(\omega_\gamma)\) to compute \(\delta H_{\mu\nu}^{20}(\omega_\gamma)\) and \(\delta H_{\mu\nu}^{02}(\omega_\gamma)\).

2. Add \(\left(E_\mu+E_\nu-\omega_\gamma\right)^{-1}\delta H_{\mu\nu}^{20}(\omega_\gamma)\) to \(X(\omega_\gamma)\) and add \(\left(E_\mu+E_\nu+\omega_\gamma\right)^{-1}\delta H_{\mu\nu}^{02}(\omega_\gamma)\) to \(Y(\omega_\gamma)\).

4.1.2.2 Lanczos.

In order to perform the matvec involving the matrix \(MK\), where \(M=A+B\) and \(K=A-B\), for the symmetric Lanczos iterations on the reduced eigenvalue problem, we need to repeat the steps of 57 . The crucial distinction in the Lanczos method is that we need to set \(\omega_\gamma=0\), so that we have a matvec with the matrix \[\begin{bmatrix} A & B \\ B & A \end{bmatrix}.\] Given a vector \(X\in\mathbb{R}^n\), the matvec with \(MK\) is computed as follows:

1. Set up a vector \(\begin{bmatrix} X \\ 0 \end{bmatrix}\) by appending \(X\) with a zero vector of size \(n\).

2. Apply the matvec to \(\begin{bmatrix} X \\ 0 \end{bmatrix}\) and obtain \(\begin{bmatrix} AX \\ BX \end{bmatrix}\).

3. Compute \(Z:=AX-BX\).

4. Apply the matvec to \(\begin{bmatrix} Z \\ 0 \end{bmatrix}\) and obtain \(\begin{bmatrix} AZ \\ BZ \end{bmatrix}\).

5. Set \(MKX:=AZ+BZ\).

We note that in addition to these two matvecs, we need one additional matvec to do \(K\)-inner product orthogonalization, another matvec for reorthogonalization, and a final one for normalization with respect to \(K\)-inner product. Thus, the Lanczos recurrence itself incurs five matvecs per iteration. In addition, after the \(m\)-step Lanczos iteration, we need to perform \(m\) matvecs to compute \(KZ_m\Theta_m^{-1/4}\) in order to derive the weights \(S_i(\hat{F})\) and \(\tilde{S}_i(\hat{F})\). Therefore, \(5m+m=6m\) matvecs are required in total.

4.2 Medium-Mass and Heavy Nuclei↩︎

The QRPA matrix dimensions are \(94,482\) for the medium-mass nucleus \(^{112}\)Sn and \(101,324\) for the heavier rare-earth nucleus \(^{150}\)Nd. We use converged frequency-by-frequency solvers as references for the smeared strength functions. We compute the Gamow–Teller strength by combining the intrinsic components of the \(J^\pi=1^+\) spin–isospin operator. In an axially-symmetric system, this amounts to performing separate calculations for the \(K=0\) and \(K=1\) components and reconstructing the total strength as the sum of the \(K=0\) contribution and twice the \(K=1\) contribution. The factor of two accounts for the degeneracy of the \(K=\pm1\) branches, which contribute equally to the total Gamow–Teller strength.

4.2.1 GMRES as a Frequency-by-Frequency Reference↩︎

We first compare GMRES with IFAM to identify a practical reference solution for the subsequent Lanczos comparison. Both methods solve the same linear-response problem at each prescribed frequency, but GMRES treats it directly as a Krylov linear solve rather than as a fixed-point iteration. Thus, when both methods are converged to the same tolerance, agreement between them is expected; the main distinction is computational efficiency. Figure 2 combines the comparison between GMRES and IFAM for both nuclei. The upper panels show that the two methods produce visually indistinguishable total Gamow–Teller strength functions for \(^{112}\)Sn and \(^{150}\)Nd, while the lower panel quantifies this agreement through the absolute pointwise difference, which remains below \(10^{-6}\) over the full frequency interval.

The iteration-count comparison in Figure 3 reveals a clear efficiency advantage for GMRES over the iterative FAM calculation. This improvement is consistent with the fact that GMRES typically requires fewer Krylov iterations to reach the same convergence tolerance. For the representative \(K=0\) calculations shown in Figure 3, GMRES requires \(2347\) total iterations for \(^{112}\)Sn compared with \(3781\) for IFAM, and \(3762\) for \(^{150}\)Nd compared with \(5336\) for IFAM over the range \(0\)–\(50\) MeV.

These results show that GMRES preserves the converged IFAM strength distribution while reducing the cost of the frequency-by-frequency calculation. For this reason, in the remainder of this section we use the converged GMRES result as the frequency-by-frequency reference against which the Lanczos approximation is compared.

4.2.2 Lanczos Approximation↩︎

We next compare the Lanczos approximation with the GMRES reference. Unlike GMRES and IFAM, the Lanczos method does not solve a separate linear system at each frequency. Instead, a single Krylov projection of the reduced eigenvalue problem provides the Ritz values and approximate transition amplitudes used to reconstruct the Lorentzian-smeared strength function over the entire interval. This subsection is organized around three questions: how the strength profile converges as the Krylov dimension increases, how close the resulting normalized distributions are to the GMRES reference, and how the runtime compares with the frequency-by-frequency methods.

4.2.2.1 \(^{112}\)Sn profile convergence.

We first consider the medium-mass nucleus \(^{112}\)Sn. To illustrate convergence with respect to the Krylov dimension, we compare Lanczos approximations with \(m=50\) and \(m=100\) iterations against the GMRES reference. Figure 4 first shows the full \(0\)–\(50\) MeV interval and then zooms in on the lower-energy windows \(0\)–\(10\) MeV and \(10\)–\(20\) MeV, where the dominant peaks are located. On the full interval, even \(m=50\) captures the broad distribution of strength. The zoomed views, however, show that \(m=50\) can miss some peak heights and local structure, whereas \(m=100\) gives visibly better agreement with the GMRES reference. The pointwise absolute difference in the lower-right panel confirms this convergence trend: increasing the Lanczos dimension from \(50\) to \(100\) reduces the discrepancy throughout the plotted interval.

4.2.2.2 \(^{150}\)Nd profile convergence.

For the heavy rare-earth nucleus \(^{150}\)Nd, we compare Lanczos approximations with \(m=100\) and \(m=200\) iterations against the GMRES reference. Figure 5 first shows the full \(0\)–\(50\) MeV interval and then zooms in on the lower-energy windows \(0\)–\(15\) MeV and \(15\)–\(25\) MeV, where the dominant peaks are located. On the full interval, \(m=100\) captures the broad distribution of strength, but the zoomed views show visible differences in local peak structure. The \(m=200\) reconstruction gives better agreement with GMRES in these zoomed regions. Once again, the pointwise absolute difference in the lower-right panel confirms this convergence trend: increasing the Lanczos dimension from \(100\) to \(200\) reduces the discrepancy throughout the plotted interval.

4.2.2.3 Distribution-level error and runtime.

As an additional distribution-level accuracy measure, we compute the Kullback–Leibler (KL) divergence between the GMRES strength distribution and each Lanczos reconstruction. Because KL divergence is defined for probability distributions, the strength values on the sampled energy grid are first normalized separately for each nucleus as \[\label{eq:KL95normalization} p_j = \frac{S_j^{\mathrm{GMRES}}}{\sum_{\ell} S_{\ell}^{\mathrm{GMRES}}}, \qquad q_j^{(m)} = \frac{S_j^{\mathrm{Lanczos},m}}{\sum_{\ell} S_{\ell}^{\mathrm{Lanczos},m}}.\tag{60}\] We then evaluate \[\label{eq:KL95divergence} D_{\mathrm{KL}}\left(p\middle\|q^{(m)}\right) = \sum_j p_j \log\left(\frac{p_j}{q_j^{(m)}}\right),\tag{61}\] where the GMRES distribution \(p\) is treated as the reference, or “true,” distribution for the corresponding nucleus. The KL divergence is nonnegative and equals zero only when the two normalized distributions agree exactly; smaller values therefore indicate better agreement with the GMRES reference. The values in Table 1 provide the same type of convergence evidence for both \(^{112}\)Sn and \(^{150}\)Nd: increasing the Lanczos dimension substantially reduces the distribution-level discrepancy. For \(^{112}\)Sn, the KL divergence decreases by more than an order of magnitude from \(m=50\) to \(m=100\). For \(^{150}\)Nd, the KL value at \(m=200\) is likewise much smaller than the value at \(m=100\), consistent with improved convergence as the Krylov dimension is increased.

Together, these comparisons isolate the convergence behavior of the Lanczos reconstruction for the two realistic test nuclei. The runtime comparison in Figure 6 shows the corresponding efficiency differences among IFAM, GMRES, and Lanczos. The improved runtime of GMRES relative to IFAM is consistent with its smaller iteration count: for the representative \(K=0\) calculations in Figure 3, GMRES requires \(2347\) total iterations for \(^{112}\)Sn compared with \(3781\) for IFAM, and \(3762\) for \(^{150}\)Nd compared with \(5336\) for IFAM over the range \(0\)–\(50\) MeV. The Lanczos timings can be interpreted through the matvec-equivalent cost model in Section 4.1.2. Including the additional products needed to form the transition amplitudes, an \(m\)-step Lanczos calculation requires \(6m\) matvecs. Thus, the Lanczos runs used for the timing comparison cost \(600\) matvecs for \(^{112}\)Sn with \(m=100\) and \(1200\) matvecs for \(^{150}\)Nd with \(m=200\). In the same bookkeeping, IFAM and GMRES require one matvec-equivalent operation per iteration, so the corresponding GMRES costs are \(2347\) and \(3762\) matvecs, respectively. These estimates explain why Lanczos remains faster than the frequency-by-frequency solvers in Figure 6, even though the precise wall-clock speedups also depend on implementation overhead and post-processing costs.

5 Conclusion↩︎

We have presented a symmetric Lanczos framework for approximating charge-changing QRPA strength functions in a matrix-free setting. Starting from the QRPA linear-response equation and passing through eigendecomposition-based and reduced-eigenproblem formulations, we obtained a representation of the strength function that is well suited to Krylov projection. The resulting method differs conceptually from conventional frequency-by-frequency FAM calculations in that a single Lanczos run captures spectral information across an entire energy interval.

The numerical experiments show that GMRES provides an efficient frequency-by-frequency reference for the medium-mass and heavy nuclei considered here, reproducing converged IFAM strength profiles while requiring fewer iterations. Against this reference, the Lanczos approximation reproduces the same overall Gamow–Teller strength profiles for \(^{112}\)Sn and \(^{150}\)Nd while requiring substantially fewer matvec-equivalent operations.

Taken together, these results indicate that the symmetric Lanczos method provides a practical alternative for QRPA strength-function calculations when broad spectral information is needed. The approach is especially attractive for deformed and heavy nuclei, where explicit QRPA matrices are too large to build and repeated frequency-by-frequency solves become expensive. Natural directions for future work include extending the method to additional external fields, refining the implementation of the reduced problem, and exploring larger systematic calculations in realistic nuclear-structure applications.

References↩︎

---

1. Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, United States↩︎

2. Department of Physics and Astronomy, University of North Carolina, Chapel Hill, NC 27599-3255, United States↩︎