QRAP: a numerical code for projected (Q)uasi-particle (RA)ndom (P)hase approximation

A computer code for quasiparticle random phase approximation-QRPA and projected quasiparticle random phase approximation-PQRPA models of nuclear structure is explained in details. An important application of the code consists in evaluating nuclear matrix elements involved in neutrino-nucleus reactions. As an example, cross section for 56Fe and 12C are calculated and the code output is explained. The application to other nuclei and the description of other nuclear and weak decay processes is also discussed.


I. INTRODUCTION
The new age of the physics beyond the standard model of electroweak interaction has as one of the most promising pathways the search of neutrino oscillations. Several experimental efforts are oriented to find the neutrino masses and the related oscillations involving atmospheric, solar, reactor and accelerator neutrinos [1][2][3][4][5]. Since neutrinos interact so weakly with matter, they bring information on the dynamics of supernova collapse and posterior explosion as well as on the synthesis of heavy nuclei [6,7]. * Electronic address: krmpotic@fisica.unlp.edu.ar, arturo˙samana@tamu-commerce.edu, carlos˙bertulani@tamu-commerce.edu The detection signal of neutrinos is measured trough the weak interaction of incoming neutrinos with the nuclei present in, e.g., a liquid scintillator detector, as well as with the surrounding blockhouse detector-shield. The flux-averaged ν-nucleus cross sections are the measured observables. Recently, Ref. [8] has studied the effect of neutrino oscillations on the expected supernova neutrino signal with the LVD detector, through their interactions with protons and carbon nuclei in a liquid scintillator and with iron nuclei in the support structure.
Charged and neutral ν e -nucleus cross sections on 12 C (liquid scintillator) as well as on 56 Fe (detector surrounding shield) were measured by the KARMEN Collaboration [9,10]. Other experiments such as LAMPF [11,12] and LSND [13,14] have also used 12 C to search for neutrino oscillations and to measure neutrino-nucleus cross sections. Furthermore, future experiments will use 12 C as liquid scintillator, such as in the spallation neutron source (SNS) at Oak Ridge National Laboratory (ORNL) [15], or in the LVD (Large Volume Detector) experiment [8].
On the other hand, the cross sections ν e (ν e )− 56 Fe are important to test the ability of nuclear models in explaining reactions on nuclei with masses around iron, which play an important role in supernova collapse [16]. The iron is used as material detector in experiments on neutrino oscillations such as MINOS [17], whereas future experiments, such as SNS at ORNL [15] plan to use the same material.
There have been great efforts on nuclear structure models to describe consistently semileptonic weak processes with 12 C such as RPA-like models. A brief summary on the different models employed for 12 C is sketched in Ref. [18].
The puzzle with the Random Phase Approximation (RPA) and the quasiparticle RPA (QRPA), when applied to the weak observables in the triad { 12 B, 12 C, 12 N}, is well known. That is, to get agreement with data for the ground state triplet T = 1 (β ± -decays, µ-capture, and the exclusive 12 C(ν e , e − ) 12 N reaction) the continuum RPA (CRPA) calculations of Kolbe, Langanke, and Krewald [19] needed to be rescaled by a reduction factor ∼ = 4. The reason for such a large discrepancy is very simple: within the RPA the transitions 12 C→ 12 N(1 + 1 ) and 12 C→ 12 B(1 + 1 ) are engendered mostly by the particlehole excitation p 3/2 → p 1/2 , what is physically incorrect. In fact, since late 1980's we know from several hadronic charge-exchange reaction measurements, and the consecutive Shell Model (SM) calculations, that the excitations p 3/2 → p 3/2 , p 1/2 → p 1/2 , and p 1/2 → p 3/2 participate quite significantly in these processes (see, for instance, [20, Table I]). It is the involvement of these configurations that brings about the necessary quenching of the Gamow-Teller (GT) resonances and β-decay rates. To make them come into play it is mandatory to open the p 3/2 shell by means of pairing correlations, which is done within both the SM and the QRPA. But, a new problem emerges in the application of the QRPA to 12 C, as first observed by Volpe et al. [21] who noted that within this approach the lowest state in 12 N irremediable turned out not to be the most collective one. As a consequence the QRPA also fails in accounting for the exclusive processes to the isospin triplet T = 1. Soon after it was shown [22][23][24] that the origin of this difficulty arises from the degeneracy among the p 1/2 and p 3/2 quasiparticle energies (both for protons and neutrons), which is inherent to the non-conservation of particle number. Therefore, for a physically sound description of the weak processes among the A = 12 iso-triplet it is imperative to use the SM or the number projected QRPA (PQRPA).
The QRAP code is based on Refs. [22][23][24], where a new formalism for neutrino-nucleus scattering has been developed, and the PQRPA is used as the nuclear model framework. The residual interaction was done with the simple δ-force, which has been used extensively in the literature to describe the single and double beta decays [25][26][27][28][29][30].
Before proceeding we address briefly on the genesis of the QRPA and PQRPA in a manner appropriate in the present context. Although this is not a topic of central interest for the application-oriented computer code, it belongs to the physics background. The neutron-proton QRPA was developed in 1967 by Hableib and Sorenson [31] in order to account for the hindrance of the allowed β-transitions. Almost 20 years later, when Vogel and Zirnbauer [32] and Cha [33] discovered the importance of the particle-particle force in the S = 1, T = 0 channel, the QRPA became to be the most frequently used nuclear structure method for evaluating double beta (ββ) rates. It was quickly realized, however, that a small change in the particle-particle interaction strength caused a large change in the lifetimes and eventually the breakdown (called a collapse) of the entire method. Later on several modifications of the QRPA were proposed to make it more reliable. One of these was the charge-exchange PQRPA, which has been formulated to evade the disadvantages inherent in the nonconservation of particle number, and was derived from the time-dependent variational principle [29]. But, the PQRPA did not yield substantially different result from the plain QRPA, and was unable to avoid the collapse in the study the two-neutrino ββ-decay in 76 Ge. As a matter of fact, the problem of the QRPA collapse has not yet been settled down, in spite of enormous effort invested for this purpose by many nuclear physicists (compare, for instance, Fig. 1 from Ref. [29] with Fig. 5 from a recent work of Yousef et al. [34]). However, the PQRPA turned out to be quite important for the description of relatively light nuclei such as 12 C. For example, the employment of PQRPA for the inclusive 12 C(ν e , e − ) 12 N cross section, instead of the continuum RPA (CRPA) used by the LSND collaboration in the analysis of ν µ → ν e oscillations of the 1993-1995 data sample, leads to an increased oscillation probability [24].
The PQRPA was recently also used to calculate the 56 Fe(ν e , e − ) 56 Co cross section [35]. A comparison between the QRPA and PQRPA for the same interaction and employing the same model space shows that the projection procedure could be important for medium mass nuclei. Moreover, several approximations such as: i) Hybrid Model (HM) [36], ii) QRPA with Skyrme interaction [37], iii) relativistic QRPA (RQRPA) [38], and iv) QRPA and PQRPA with the δ-force [35] yield different results for the neutrino cross section as a function of the neutrino energy. It is a hard task to find the origin for the differences, mainly because these models are not using the same interaction and/or the same single-particle configuration space, carrying different types of correlations in each case.
The cross sections for charged-and neutral-current neutrino-induced reactions on the iron isotopes 52−60 Fe were also evaluated within the HM for various supernova neutrino spectra [39]. Here, large-scale SM calculations were used for the GT-like contributions, while transitions for other multipoles are based on the RPA. More precisely, the authors scale the SM cross sections using the ratios obtained from the RPA calculations with and without this dependence of the multipole operator. The reason for such a procedure is twofold: i) the limitation of the SM to account for momentum-transfer dependence of the GT operator, and ii) the lack of pairing correlations in the RPA. It should be also mentioned that SM calculations of inelastic neutral-current neutrino-nucleus cross sections in medium-mass nuclei, present in supernova environment, have been constrained by the highly precise data on the magnetic dipole strength distributions for the nuclei 50 Ti, 52 Cr, and 54 Fe, which are dominated by spin-isospin flipping (GT-like) contributions [40]. In spite of the agreement between data and calculations it was necessary to consider also here the effects of finite momentum transfer what was done via the RPA. Briefly, the HM is neither fish nor fowl, and a comparison of the results from Refs. [39,40] with self-consistent calculations, such as the QRPA, PQRPA and RQRPA, could be enlightening.
This brief introduction shows: 1) the importance of neutrino-nucleus cross sections for astrophysical purposes and, 2) that these cross sections are strongly correlated with the nuclear structure model employed. The QRAP code, with a simple residual interaction, is able to access the sources of these problems and it can calculate several weak interaction processes mentioned above. Needless to stress that this code can be easily adapted for the evaluation of ββ-decays.
The write-up is organized as follows. In section II we make a short survey of the theoretical description of weak interaction processes, with emphasis on the formulation implemented in this numerical code. In sections III and IV we describe the QRPA, and PQRPA formalisms, making explicit the differences among them. In section V we show how the code is organized, how to make an input and how to understand the output. Section VI explains the role of each subroutine of the code. Finally, section VII proposes a few cases to practice with the code.

II. WEAK INTERACTING PROCESSES
In this section we give a brief summary of the main formulae developed in Ref. [18,23] for: • neutrino scattering (NS) • antineutrino scattering (AS) • muon capture (MC) rate where ℓ = e, µ. The comparison with other formalisms [41][42][43] can be found is in just mention works.
The weak Hamiltonian is expressed in the form where G = (3.04545 ± 0.00006)×10 −12 is the Fermi coupling constant (in natural units), is the hadronic current operator 1 , and is the plane wave approximation for the matrix element of the leptonic current in the case of neutrino reactions, with p ℓ ≡ {p, iE ℓ } and q ν ≡ {q, iE ν } being, respectively, the lepton and the neutrino momenta.
For the sake of convenience we will use spherical coordinates (m = −1, 0, +1) for the three-vectors, and the Walecka's notation [42], with the Euclidean metric, for four-vectors, i.e., x = {x, x 4 = ix ∅ }. The only difference is that we substitute Walecka's indices (0, 3) by our indices (∅, 0), i.e. we use the index ∅ for the temporal component and the index 0 for the third spherical component.
The quantity is the momentum transfer, where P i and P f are momenta of the initial and final nucleus, M is the nucleon mass, m ℓ is the mass of the charged lepton, and g V , g A , g M and g P are, respectively, the vector, axial-vector, weakmagnetism and pseudoscalar effective dimensionless coupling constants. Their numerical values are: In the numerical calculations we use an effective axialvector coupling g A = 1 [44]. The finite nuclear size (FNS) effect is incorporated via the dipole form factor with a cutoff Λ = 850 MeV, i.e., (2.6) To use (2.1) with the non-relativistic nuclear wave functions, the Foldy-Wouthuysen transformation has to be performed on the hadronic current (2.2). When the velocity dependent terms are included this yields [45]: (2.7) wherek = k/κ, κ ≡ |k|, and v ≡ −i∇/M is the velocity operator, acting on the nuclear wave functions. The following short notation has also been introduced.
In performing the multipole expansion of the nuclear operators (2.9) it is convenient: 1) to take the momentum k to be along the z axis, i.e., where ρ = κr, and 2) to introduce the operators O αJ , defined as Thus, where the geometrical factors are listed in Table I of Ref. [23]. Explicitly, from (2.7) The elementary operators are given by Here we make use of the conserved vector current (CVC). From (2.14), (2.15), and [46,Eq. (10.45) and are not, and it is convenient to put in evidence their real and imaginary parts, expressing them as with M R 1J , and M I 1J arising, respectively, from the terms in (2.16) with L = J ± 1, and L = J. Note that F ±1JJ = ∓1/ √ 2. It is also convenient to separate the elementary operators into: being a NP (UP) operator. This is a very important finding because it implies that O R αJ and O I αJ do not contribute simultaneously, and, therefore, one always can deal only with real operators.
In summary, natural and unnatural parity operators are, respectively: and

22)
A. Neutrino-nucleus cross section For the neutrino-nucleus reaction, the momentum transfer is k = p ℓ − q ν , and the corresponding cross section reads where F (Z ± 1, E ℓ ) is the Fermi function (Z + 1, for neutrino, and Z − 1, for antineutrino), θ ≡q ·p is the angle between the incident neutrino and ejected lepton, and the transition amplitude is (2.24) After expressing the spatial part of the lepton traces L αβ in spherical coordinates, and applying the Wigner-Eckart theorem, one can cast the transition amplitude in the compact form [23] T The explicit expressions for the traces L ∅ ≡ L ∅∅ , L m ≡ L mm , and L ∅0 are [23] L ∅∅ = 1 + |p| cos θ E ℓ , being the z-components of the neutrino and lepton momenta, and S 1 = ±1 for NS and AS, respectively.

B. µ-capture rates
The muon capture transition amplitude T M C (J f ) can be derived from the result (2.25) for the neutrino-nucleus reaction amplitude, by keeping in mind that: i) the roles of p and q are interchanged within the matrix elements of the leptonic current, which makes that in (2.26) S 1 → −1, ii) the momentum transfer turns out to be k = q − p, and therefore the signs on the right-hand sides of (q 0 , p 0 ) have to be changed, and iii) the threshold values (p → 0 : q → k, k ∅ → E ν − m ℓ ) must be used for the lepton traces. All this yields q 0 = E ν , p 0 = 0, and (2.28) Instead of summing over the initial lepton spins s ℓ , as done in (2.24), one has now to average over the same quantum number. We get where φ 1S is the muonic bound state wave function evaluated at the origin, and is the binding energy of the muon in the 1S orbit. Thus from (2.25) and (2.28) In the case of MC it is convenient to rewrite the effective coupling constants (2.8) as Thus, natural and unnatural parity operators are now, respectively: The PQRPA for charge-exchange excitations was derived from the time-dependent variational principle in Ref. [29]. In the same reference is also described in details the projected Barden-Cooper-Schiffer (PBCS) approximation. Basically one employs the number projection operatorsP N on the |BCS state. That is:P 0 =P ZPN for a ground state with Z protons and N neutrons, and P µ =P Z+µPN −µ , with µ = ±1, for excited states in nuclei with Z + µ protons and N − µ neutrons. In this section we give a brief description of both the PBCS and PQRPA approximations.
The PBCS gap equations are are the pairing gaps, and are the dressed single-particle energies, where are the PBCS number projection integrals. The PBCS correction term ∆e k can be found in Ref. [29], G(kkk ′ k ′ ; 0), and F(kkk ′ k ′ ; 0) stand for the usual proton or neutron particle-particle (pp), and particle-hole (ph) matrix elements of the residual interaction V , i.e., G(klmn; J) = kl; J|V |mn; J . F (klmn; J) = kl −1 ; J|V |mn −1 ; J . (3.5) Note that these relations are valid for both identical and non-identical particles. The forward-going (X µ ), and backward-going (Y µ ) PQRPA amplitudes are obtained by solving the RPA equations with the PQRPA matrices defined as: where are the unperturbed energies, are the norms, are the projected quasiparticle energies, and the quantities R K are defined as [29] (3.11) Both positive and the negative solutions are physically meaningful. For µ = ±1 the positive solutions describe excitations in the (Z ±1, N ∓1) nuclei, while the negative energy solutions represent the positive energy excitations in the (Z ∓1, N ±1). Thus only one RPA equation has to be solved, either for µ = +1, or for µ = −1, to describe the excitations to the Z ± 1, N ∓ 1 nuclei. This is well known feature of the charge-exchange modes [47][48][49][50].
Let us be more specific, and take advantage of the index f to label different final states |J π f with same spin and parity. Evidently, f will run from 1 up the total number f max of two-quasiparticle configurations |pnJ π . Moreover, the eigenvalue problem (3.6) has 2f max solutions, and we will use the index F to label them. Thus, if µ = +1 one have: Finally, to store the eigenvalues and eigenfunctions it is convenient to define the index F as The usual gap equations are obtained from Eqs. (3.6)-(3.7) by: 1. Making the replacement e k → e k − λ k , with λ k being the chemical potential, and taking the limit I K → 1. That is, the Eq. (3.1) remains as it is, but instead of (3.2) and (3.3) one has now and (3. 16) 2. Impose the subsidiary conditions as the number of particles is not any more a good quantum number.
In this way the usual BCS gap equations read This equation, together with the normalization condition u 2 k + v 2 k = 1, has as solution the occupation probabilities (for example, from the Chapter I of Rowe [47]) which depend on the quasiparticle energies 20) and the pairing gaps The QRPA equations are recovered from (3.6) by i) dropping the index µ, ii) taking the limit I K → 1, and iii) substituting the unperturbed PBCS energies by the BCS energies relative to the Fermi level, defined by equation where E k are the usual BCS quasiparticle energies defined in (3.20). In this way the unperturbed energies in (3.7) are replaced by These energies, however, are not used in the QRPA eigenvalue problem. Namely, the coefficients X(pnJ f ) and Y (pnJ f ), and the eigenvalues ω(J f ) are obtained from where In the pn-QRPA the eigenvalues occur in pairs ±ω(J f ), but the negative energies don't have a direct physical meaning. The perturbed energies for daughter (Z + µ, N − µ) nuclei are defined as There is, however, only one set of eigenfunctions (X(pnJ f ), Y (pnJ f )) for both µ = 1, and µ = −1. This is a very important difference in relation to the PQRPA case, which is crucial for the distribution of the transition strengths.

C. Nuclear Matrix Elements
When the excited states |J f in the final (Z ± 1, N ∓ 1) nuclei are described within the PQRPA, the transition amplitudes for the multipole charge-exchange operators (2.21) and (2.22) read with the one-body matrix elements given by 28) and J = J.
In the QRPA case, using the limit I K → 1 in (3.27), the nuclear matrix elements for the multipole chargeexchange operators O J are with the same one-body matrix elements (3.28).
The unperturbed and perturbed transition strengths are defined, respectively, as One might be particularly interested in the Gamow-Teller (GT) and Fermi (F) β-decay strengths (B-values), in which case (3.30) and (3.31) are evaluated for the operators O J , which don't contain the radial form factors j J (ρ). That is, O 0 = 1, and O 1 = σ, for F and GT operators, respectively. We denote the B-values as S µ (pnJ) and S µ (J f ). Occasionally one also might want to calculate the the energy distribution of the last one, i.e., where one usually takes η = 1 MeV.

IV. COMPUTER PROGRAM AND USER'S MANUAL
The QRAP code evaluates the electron neutrinonucleus interaction described by equation (2.1) (IREAC=1 for NS, IREAC=2 for AS) and (2.2) (IREAC=0 for MC). The processes, from the ground state of the even-even father nucleus (Z, N ) to the excited states with spin (ISPIN) and parity (IPARI) in the odd-odd daughter nucleus (Z ± 1, N ∓ 1), are calculated by using the QRPA model (IQP=0) or the PQRPA (IQP=1). These options must be setup in the input data file, qrapin.dat, which is supplemented with two included files: (a) sp.inc, containing the dimensions of single-particle quantum numbers, occupation probabilities, quasiparticle quantities and strength amplitudes for allowed transitions; (b) conf.inc, which has the dimensions for the quasiparticle state configurations, the hamiltonian matrices (A, B) or (A, B) which are diagonalized, the forward and backward amplitudes, and the eigenvalues.
There are two input files: 1) qrapout.dat, where is listed the output file that shows the neutrino/antineutrino (ν/ν) cross section, as a function of the incident neutrino, or the muon capture rate, for each state of a given nuclear spin in the daughter nucleus, and 2) the above mentioned qrapin.dat, which contains: a) the quantum numbers of all single-particle state (sps), and the corresponding single particle energies (s.p.e.), b) the mass and the proton number of the parent nucleus, c) the neutron and proton pairing strengths for the BCS approximation, d) the particle-particle, and particle-hole strengths of the residual interaction, e) the position of the Fermi level, and the experimental gap for neutrons and protons, and f) the Q-value for the ν/ν scattering.
There are three default output files. Two of them, AUXI.OUT and OUT.OUT, contain the results of the nuclear structure model, whereas the results for the weak processes appear in the file created by qrapout.dat. For example, if one is interesting in the multipole J π f = 1 + with a single-particle space of six levels in 12 C ("set 1"), we can introduce in qrapout.dat the file names QNC.out (PNC.out) for neutrino capture, QAC.out (PAC.out) for antineutrino capture, QMC.out (PMC.out) for muon capture, using the QRPA (PQRPA) model. The auxiliary output files AUXI.OUT and OUT.OUT are relabeled to (QAUXI.OUT, QOUT.OUT) and (PAUXI.OUT, POUT.OUT) for QRPA and PQRPA respectively.
All just mentioned outputs are included as examples. The following units are employed: i) 10 −42 cm 2 , for neutrino or antineutrino-nucleus cross sections, ii) 10 4 s −1 , for muon capture rates, and iii) MeV, for energies.

A. Reading the data
There are three sets of input data in qrapin.dat separated in modules labeled as: *Data set 1 for a singleparticle space of six levels in 12 C (0, 1, and 2 ω oscillator shells), *Data set 2 for a single-particle space of ten levels in 12 C (0, 1, 2, and 3 ω oscillator shells), and *Data set 3 for single-particle space of 12 levels in 56 Fe (2, 3, and 4 ω oscillator shells).
For each one of these input data, the number of sps represents the available space where one wants to solve the BCS (or PBCS) problem given by equations (3.17) and (3.18) ((3.1)-(3.2)). It contains the necessary number of harmonic oscillator shells leading to a smooth smearing of the Fermi's surface. The Fermi level with the neighboring levels constitute the active shell for the mentioned smearing. For example, in 12 C (ground state with J = 0 + ) the active shell is composed by the 1p 3/2 and 1p 1/2 levels. According to the single-particle shell model the sps filled up the 1p 3/2 orbital, and the nucleons can be promoted to the 1p 1/2 , creating a particle state in 1p 1/2 and a hole state in 1p 3/2 . This scheme describes the first particlehole (ph) excitation on 12 C in order to obtain the 12 N or 12 B ground state with J = 1 + , by promoting a proton or a neutron, respectively.
Let us show, as example, a data input of six sps for 12 C: *Data set 1. The rows starting with a symbol "*" are not read as input and just serve to remind the user on the meaning of the physical quantities. Taking out the comments "*" in the first lines of this file, we have First line: Nuclear spin (ISPIN=1) of the daughter nucleus, parity (IPARI=+1), coupling strength of particle-particle channel (t = 0.0) (see definition below), index of neutrino reaction (IREAC=1) and the index (IQP=1) to solve the PQRPA problem.
Third and Fourth lines: quantum numbers and the s.p.e. for each neutron and proton sps, respectively. They are represented in the same way as in the shell model scheme, with their respective quantum numbers (n ℓ (j + 1 2 )). For example, 101 → 1s 1/2 where 1 is principal quantum number, corresponding to the first harmonic oscillator level (n), 0 corresponds to the orbital angular momentum ℓ ≡ s and the last number 1 ≡ j + 1 2 = 1 2 + 1 2 . Table I shows the notation and the corresponding quantum numbers, as well as the PBCS quasiparticle energies e N j , and e Z j , defined by (3.10).   1s 1/2 1p 3/2 1p 1/2 2s 1/2 1d 5/2 1d 3/2 and sixth lines respectively. Eighth line: Q-value minus the lepton mass for ν/ν scattering (EGS=17.338 for 12 N [51]). It can be fixed as being the energy of the ground state in the daughter nucleus. The lepton mass must be added to EGS to obtain the Q-value for the reaction.

B. Running the code
As first step the QRAP solves the BCS problem. In this case, one needs to adjust the pairing strength to reproduce the experimental pairing gap.
Next, one can solve the PBCS problem or directly the QRPA if the option IQP=0 was selected. If IQP=1 then the PQRPA equations are solved. It means that QRAP firstly calculates the nuclear matrix elements in the QRPA or PQRPA by selecting the option IQP=0 or 1, appropriately. The option for which type of weak interaction process one wants to evaluate is adopted with IREAC in the input data. We recommend first to adjust the pairing strength as it is explained below. After this it is convenient to fit the parameters of the residual interaction using the option IREAC=3 for the muon capture rate because this calculation is fast. Physically you can check quickly how good is your choice of parameters be-cause the values for inclusive muon-capture rate, and GT B-values are available in the literature (see for example Refs. [52,53]).
For the residual interaction the code assumes a delta force, which has been used extensively in the literature [26][27][28] to describe single and double beta decays. Next, we explain how the parameters of the interaction are adjusted using for example the input data for six levels in 12 C. The results are presented in output file OUT.OUT. • where B(N, Z) is the binding energy of the even-even nucleus (Z, N ). This is the most common fit (Fit1) used in several works with standard QRPA [28][29][30]. In this case, the ∆ N (Z) must be equal or approximately equal to the energy ∆ N (Z) j k =F L of the corresponding to the Fermi level (FL). To solve the set of PBCS coupled equations (3.1)-(3.4) for u k and v k it is recommended to obtain first the solutions for the BCS problem, as these probability occupations are use as input for the PBCS case. The PBCS coupled nonlinear equations are solved consistently with Powell Hybrid method using subroutine HYBRD [55].
The results of the BCS or PBCS problem are shown as tables in the first lines of OUT.OUT for neutrons and protons, respectively. The quantities defined by (3.10) and (3.11) are presented there. In particular, the projected quasiparticle energy defined in (3.10) are which means that E(+) corresponds to a particle state, and E(−) to a hole state. The values of ∆ N (Z) j k are shown in the ninth column of the table labeled as CONFIGU-RATION SPACE. This Fit1 comes from the fact that the experimental energy difference between the states that lie just above (p state) and just above (h state) the FL is approximately twice the experimental gap, i.e., There is another fitting procedure for the pairing gap that is called by Fit2. In Fit2, all the s.p.e. e N (Z) j from Table I are varied with a χ 2 search to account for the experimental spectra E j : p ≡ E(+), for a particle state, In Fit 2, the Eq. (4.4) is automatically satisfied. This procedure was employed to obtain the e j spectra shown in [23, Table III], whereas the e j for the reduced space of six levels in the present example are shown in Table I. These s.p.e. are used in input data Data set 1.
To make the calculations as simple as possible the Fit1 procedure is the usual choice, with the e j spectra obtained either from a harmonic oscillator or from a Wood-Saxon potential, and by varying the coupling v pairN s and v pairP s to satisfy the condition ∆ For 56 Fe the input data is called *Data set 3. The s.p.e of the active 3 ω shell were taken from the experimental energies of 56 Ni, and the 2 ω and 4 ω shell energies were taken from the harmonic oscillator energies with ω/MeV = 45 A 1/3 − 25 A 2/3 . Fit1 was employed to adjust the experimental ∆ N (Z) for 56 Fe.
• Adjusting the particle-hole couplings r and p In the particle-hole matrix element F , defined in Eq. (3.5), the couplings v s and v t appear as linear combinations v s + v t and 3v t − v s . Therefore, it is convenient (4.7) • Adjusting the particle-particle couplings s and t Here also it is convenient to normalize to v pair s the coupling constants v pp s and v pp t that appear in the pp matrix elements G in Eq. ∼ v pp s [25,30]. However, in Ref [22] it was shown that this parametrization might not be suitable for N = Z. In fact, the best agreement with data in 12 C was obtained when the pp-channel is totally switched off, i.e., v pp s ≡ v pp t = 0, and three different set of values for the ph-coupling strengths were used. These conditions are related with s = t for 12 C (N = Z), and with s = 1 and t variable in nuclei with N > Z. For 56 Fe were adopted the values s = 1 and t = 0. In the code QRAP the following conditions are standard: (i) s = t with t as a variable parameter for N = Z; and (ii) s = 1, and t as a variable parameter for N > Z, i.e., the residual interaction is defined as a function of two adjustable parameters v ph t and t.
Several experimental data are available in the literature that can be used for fixing the residual interaction coupling constants, such as: ground state energies of daughter nuclei, B(GT )-values for the β + or β − decay, and partial muon capture rate [51,58,59].
One can use the reduced space of six levels to identify in the output file, the quantities shown in Figure 1. The results for three values of t are shown in Table III. The values of ω µ (1 + f ), and S µ (1 + f ) in 12 N and 12 B can be found in the output file AUXI.OUT. In the present case the largest value of index f is f max =16. Both set of states, with µ = +1, and µ = −1, are ordered from highest to lowest energies. In the PQRPA, the most collective ones are that of the corresponding ground states: |1 + F =16 in 12 N (and |1 + F =17 in 12 B) although there also are significant strengths in the states F = 7, 11, and 14. In QRPA, the ground state is in |1 + F =16 for both 12 N and 12 B. These wave functions are presented below. For the PQRPA case, we also show the unperturbed energies ω 0 µ (pn1 + ) (which are not ordered), and the corresponding singleparticle GT strengths S 0 µ (pn1 + ), given respectively, by (3.8), and (3.30) for the GT operator σ. The largest ones are S 0 +1 (1p π 1/2 , 1p ν 3/2 ; 1 + ), and S 0 −1 (1p π 3/2 , 1p ν 1/2 ; 1 + ), which in the particle-hole limit correspond to excitations 1p ν 3/2 → 1p π 1/2 , and 1p π 3/2 → 1p ν 1/2 . For spins and parities J π f = 1 + , or 0 + , in AUXI.OUT are shown the energies ω µ , but not the strengths S µ . The results for the eigenvalue problem are displayed in the output file OUT.OUT. For the option MAPR= 1 are printed out the matrix elements (A, B) Eq. (3.25) for the QRPA, or (A µ=1 , B) eq. (3.7) for PQRPA. The nuclear wave functions (X(pn; J F ), Y (pn; J F )) are grouped to four, with the index F , defined in (3.14), going from 1 to f max in the QRPA case, and from 1 to 2f max in the PQRPA case. To make easy reading together with each set of wave functions are also printed: the value of f , the two quasiparticle configurations (p and n), and the unperturbed and perturbed energies.
• Output for the ν-nucleus processes The output of the results for the weak processes is selected according to the value of IREAC: IREAC=0 prints the results for the muon capture rate in the file (QMC.out or PMC.out). For J π = 0 + or J π = 1 + are shown in in this output file the folded strengths S µ (J f , E) (SˆTILDE) defined by (3.32), where 'ENERGY' represents E. The partial capture rate for each state f , the perturbed energy ω µ=−1 , and the strength S µ=−1 (J f ) (if J π = 0 + , or J π = 1 + ) are shown in the table labeled CAPTURE RATE. The total capture for the evaluated spin J π is presented in the last line.
IREAC=1 or IREAC=2 prints the results for the neutrino or antineutrino cross sections in the files (QNC.out/PNC.out) or (QAC.out/PAC.out). We repeat in this output the folded strengths S µ (J f , E) (3.32) for J π = 0 + or J π = 1 + . The cross sections (SIGMA(Enu)) are calculated as a function of the neutrino energy (Enu) for each nuclear spin from f = 1 to f = f max . The perturbed ω µ=±1 energies for the daughter nucleus are also shown according the process related. The absolute value of maximum (cos θ = −1) and minimum (cos θ = 1) nuclear momentum transfer (|k|) in units of MeV/c for each energy are also printed.
Note: The cross sections are printed up to a maximum energy of 250 MeV. Depending upon the single particle space employed, the cross sections, as a function of the neutrino energy, should be restricted to lower energies. This issue will be discussed and explained in details in a next work [60]. Anyway, the PQRPA cross sections obtained within the single particle space provided as examples are well behaved up to E ν/ν < 100 MeV on averaged according to *Data set 1, *Data set 2 and *Data set 3. This interval of energies is important for supernova neutrinos and low-energy decay-at-rest neutrinos [10] .
For the RSMPE dependent on the tensor product of spherical harmonic times the nucleon velocity we have We use here the angular coupling |( 1 2 , l)j ,Ĵ ≡ √ 2J + 1 and (l p L|l n ∓ 1) is the short notation for the Clebsh-Gordon coefficient (l p 0L0|(l n ∓ 1) 0).
For the scalar product of spin times nucleon velocity, we have p, (l p u np,lp (r)u nn,ln (r) j J (κr) r 2 dr.
Finally, the RSMPE dependent of the tensor product of spherical harmonic times the spin operator reads with the radial part R 0 L (pn; κ) given by (6.6).

B. Fermi function and effective momentum approximation (EMA)
To account for the Coulomb interaction between the charged lepton and the residual nucleus, the QRAP code is setup to use by default the Fermi function [45,46]. This correction was employed in several works for reactions on 12 C with neutrinos from the DAR of µ + . As pointed out in Ref. [21], the quantity p ℓ R A is of the order of 0.5, where p ℓ is the lepton momentum, and R A is the radius of the nucleus. Thus, the correction is well described by a Fermi function. Yet, for high energy neutrinos, e.g. neutrinos from the DIF of π + , the outgoing muons have p ℓ R A > 0.5. For these relativistic leptons, the effective momentum approximation (EMA) [66] should take care of the Coulomb field of the daughter nucleus, instead of the Fermi function. This prescription for the Coulomb correction is considered in the code within the subroutine SECCION. More precisely, with EMA=0 the Fermi function is employed, while with EMA=1 the EMA prescription is used. In the EMA procedure, the lepton energy and momentum are modified by a constant electrostatic potential within the nucleus with V ef f = 4V C (0)/5 = −6Z f α/5R A [38,67]. These two approximations for the Coulomb correction were tested in the calculation of the inclusive cross section for neutrino scattering on 208 Pb [38]. As shown in Ref. [38], the Fermi function correction overestimates the cross sections at higher neutrino energies where the EMA provides a more reliable approach. Thus, we recommend to use the Fermi function correction in the range of neutrino energies for which the cross section is below the corresponding EMA value, whereas the EMA could be employed at higher energies, as shown in previous studies [21,38].
As a final comment, the QRPA code could be easily extended to calculate ν µ -induced processes. This was done in Refs. [22,23] to calculate ν µ − 12 C cross sections using the EMA prescription for the DIF regime of the LSND experiment. The nuclear structure calculations remain the same, while the kinematics changes by changing the electron mass to the muon mass in the variable RMLEP.