Pion-nucleon correlations in finite nuclei in a relativistic framework: effects on the shell structure

The relativistic particle-vibration coupling (RPVC) model is extended by the inclusion of spin- and isospin-flip excitation modes into the phonon space, introducing a new mechanism of dynamical interaction between nucleons with different isospin in the nuclear medium. Protons and neutrons exchange by collective modes which are formed by isovector $\pi$ and $\rho$-mesons, in turn, softened considerably because of coupling to nucleons of the medium. These modes are investigated within the proton-neutron relativistic random phase approximation (pn-RRPA) and relativistic proton-neutron time blocking approximation (pn-RTBA). The appearance of isospin-flip states with sizable transition probabilities at low energies points out that they are likely to couple to the single-particle degrees of freedom and, in addition to isoscalar low-lying phonons, to modify their spectroscopic characteristics. Such a coupling is quantified for the shell structure of $^{100,132}$Sn and found significant for the location of the dominant single-particle states.


I. INTRODUCTION
A consistent treatment of pion degrees of freedom is still one of the challenges for nuclear structure theory, despite the fact that the meson-exchange nucleon-nucleon interaction is known since the work of Yukawa [1]. Being a mediator of the medium-range nucleon-nucleon interaction, the pion played the central role in approaches developed by the so-called 'Moscow school' [2][3][4][5][6][7][8][9][10][11] and 'Stony Brook school' [12][13][14][15][16][17][18][19]. Characteristics of the pion bound in nuclear medium were extensively discussed, in particular, in the context of possible formation of pion condensate in finite nuclei and various phases of nuclear matter [11,18,19]. Later on, chiral effective field theory (EFT) based on the approximate spontaneously broken chiral symmetry of QCD emerged from the pioneering work of Weinberg [20]. Since then, ab initio methods based on chiral EFT demonstrate an increasing success in the description of nuclear structure properties although their applicability is still limited to light nuclei only [21,22].
The covariant (relativistic) version of the nuclear field theory [48][49][50][51][52][53][54][55][56][57][58] is based on the quantum hadrodynamics (QHD) Lagrangian, which, in contrast to chiral EFT, includes heavy mesons explicitly, although their masses and coupling constants are adjusted to bulk properties of medium-mass nuclei (see review [59] and references therein). Numerous calculations within these models have shown that the effects of quasiparticlevibration coupling (QVC) modify nuclear shell structure and spectra of excitations considerably and bring the results in essentially better agreement to experimental data, in some cases approaching spectroscopic accuracy [50,52,54,60,61]. In particular, it turned out that the effects of QVC in open-shell nuclei are quantitatively stronger than those in closed-shell ones, which can be attributed to the existence of phonons with very low energies (∼1 MeV or even lower) in medium-mass open-shell nuclei, compared to ∼4 MeV in medium-mass doublymagic systems. Obviously, superfluidity plays a very important role providing better conditions for the QVC, which can be seen from the structure of the QVC selfenergy [50]. The pion field does not contribute on the mean-field level, if the parity conservation is imposed, however, its dynamical contribution is decisive for the description of spin-isospin transitions [56,57,62,63].
Recently, it has been realized that the spin-isospin interaction, especially the tensor interaction, is of particular and direct importance for single-particle properties of exotic nuclei [64,65]. This type of interaction is com-monly associated with one-pion exchange, although other mesons can also give contributions of the tensor type, which has been obtained in the relativistic Hartree-Fock [66] and relativistic RPA (RRPA) calculations for spinisospin-flip transitions [67]. However, the latter calculations are confined by the static approximation for the pion field, neglecting retardation effects, and include only the first-order contribution to the pion-exchange interaction. Nevertheless, RPA in the spin-isospin channel gives a sizable softening of the pion modes, where the softness is caused by medium polarization effects, as it has been recognized already in earlier calculations [11]. The appearance of such modes at low energies was interpreted as a signature for the onset of the pion condensation, although it was not clear which observables are sensitive enough to the presence of pion condensate and which of those can be detected experimentally with reasonable precision.
The main goal of the present work is to reveal dynamical contributions of the pion-exchange interaction to nuclear structure properties. The most important dynamical contributions are identified as terms in the nucleonic self-energy, which represent exchange by collective modes. Technically, besides the isoscalar modes, which are exchanged between nucleons of the same isospin and commonly considered in NFT, the isovector modes exchanged between nucleons with different isospin are added. First of all, the approach of Refs. [56,62], namely the proton-neutron relativistic random phase approximation (pn-RRPA) and proton-neutron relativistic time blocking approximation (pn-RTBA) is applied to calculations of the isospin-flip and spin-isospin-flip excitations. The obtained spectra are evaluated to identify the lowlying isovector modes of excitation, which give the most important contribution to the nucleonic self-energy, according to their degree of collectivity. Then, the selected most collective modes are coupled to the single-particle degrees of freedom by means of introducing corresponding terms into the single-nucleon self-energy. Thus, on one hand, the relativistic NFT is extended by isospin vibrations and, on the other hand, dynamical effects of the pion exchange are introduced into the theory beyond the Hartree-Fock approximation. The effects of these new terms on the single-particle characteristics, such as their energies and spectroscopic factors, are discussed. In the context of the discussion of pion condensate, it is shown that nuclear single-particle characteristics can represent observables which are sensitive to the presence of soft pion modes.

II. THEORETICAL FRAMEWORK
A. Nucleonic self-energy and Dyson's equation The motion of protons and neutrons in nuclei is quantified by one-nucleon self-energy, or mass operator. In contrast to free nucleons, protons and neutrons in the nuclear medium are exposed to strong polarization effects and, hence, their motion is modified considerably. In the simplest mean field (Hartree or Hartree-Fock) approximation, this effect can be described by introducing an energy-independent effective mass m * instead of the bare nucleon mass, which helps to reproduce bulk nuclear properties reasonably well, but gives a very poor description of the single-particle spectra. The latter problem can only be solved by going beyond the mean-field approximation and by introducing an energy-dependent effective mass. The most important origin of the energy dependence is given by the coupling of the single particle motion to low-lying collective vibrations [24][25][26][27], see also the discussion in Ref. [48], where a relativistic particle-vibration coupling (RPVC) model was originally proposed. In both the RPVC model and relativistic quasiparticle-vibration coupling (RQVC) model for systems with pairing correlations [50], the energyindependent part of the nucleonic self-energy is given by the relativistic mean field (relativistic Hartree (RH) or Hartree-Bogoliubov (RHB)) approximation. The RH/RHB eigenstates {k, η} are used as a basis for calculating the dynamical part of the self-energy, matrix elements of which in this basis read: where the indices k, k 1 , k 2 represent the complete set of the single-particle quantum numbers, and the indices η, η 1 , η 2 indicate the upper and lower components of the matrix elements in the Nambu space. In the limit of no superfluidity η 1 = 1, if the state k 1 is above the Fermi energy, and η 1 = −1, if the state k 1 is below the Fermi energy. The index µ in Eq. (1) runs over the set of phonons taken into account. Their vertices γ µ and frequencies Ω µ can be calculated in the self-consistent relativistic quasiparticle random phase approximation, as described in Ref. [52], and E k are the energies of the Bogoliubov's quasiparticles obtained in the RHB approach, so that the product ηE k takes the correct limit of the mean-field energy in the case of no superfluidity. The index η µ labels forward and backward going diagrams in the self-energy (1), see Eq. (7) below. Using the same representation, the Dyson's equation for the single-quasiparticle Green's function can be formulated as follows: By solving this equation, as described in detail in Ref. [48], one can obtain the energies of the fragmented states E (ν) k and the corresponding spectroscopic factors S (ν) k , which characterize the probability of occupying the levels E (ν) k by quasiparticles. In the following we will focus on the isospin structure of the self-energy of Eq. (1). It represents the lowest-order phonon polarization correction to the singleparticle propagator and contains already an infinite sum of the diagrams dressing the fermionic line. If the summation is performed for the ring type of diagrams, the vertices γ µ and frequencies Ω µ are approximated by an R(Q)RPA phonon solution. Formally, the sum in Eq. (1) runs over all possible intermediate nucleonic states k. The common practice is to include only the states k of the same isospin as the one of the external states k 1 , k 2 . In the relativistic framework, this practice leads to a very good description of the single-quasiparticle spectra in open-shell nuclei, compared to data, and to a considerable improvement of the results in the closed shell nuclei, although in the latter case the description of the energies have not yet reached spectroscopic accuracy. In this context, the inclusion of components in the self-energy connecting states of different isospin can serve as the first step of the search for possible missing mechanisms of PVC.
The next subsections are devoted to the protonneutron phonons in nuclei and to establishing their characteristics.

B. Isospin-flip modes in 100,132 Sn
Isospin-flip modes of excitation in nuclei are the most sensitive probes to reveal the properties of the isospin part of nucleon-nucleon interaction. Besides that, precise knowledge about various strength distributions involving isospin transfer, especially in exotic nuclei, is of great importance for astrophysics as well as for numerous applications.
Another aspect of isospin dynamics is its direct connection to pionic degrees of freedom, that can be revealed in calculations which include pion-nucleon interaction explicitly. Already in earlier works (see Ref. [11] and references therein) it has been realized that the explicit separation of pion-exchange processes lead to a considerable improvement of the description of isospin transfer modes, beta-decay half-lives and magnetic momenta.
Lately, self-consistent calculations in the relativistic framework have become possible on the (Q)RPA [62,63,67] and beyond QRPA [56,57] levels, using the free-space pseudovector coupling strength for pionnucleon interaction. In particular, a very reasonable description of the position of the Gamow-Teller resonance (GTR) and beta-decay half-lives in medium-mass nuclei has been achieved already on the QRPA level without a non-relativistic reduction of the pion-exchange interaction [62,63]. Besides that, it has been shown that the explicit inclusion of the meson exchange (Fock) term in both the mean-field and effective interaction gives a relatively good description of the GTR's centroid on the level of self-consistent RPA [67]. Calculations beyond R(Q)RPA demonstrate an important role of the phonon coupling for spreading of the GTR [57] and spin-dipole resonance (SDR) [56]. Analogously to the situation for giant resonances in the neutral channel, the inclusion of the phonon coupling mechanism shows a considerable improvement of the description for both GTR and SDR.
SDR occurs as a nuclear response to the spin-isospin flip operator with an additional one unit angular momentum transfer, i.e. to the operator: with L = 1, λ = 0, 1, 2 and parity π = −1. Components of different values of λ can be disentangled in chargeexchange reaction experiments on stable nuclei and show a hierarchy of the λ = 0, 1, 2 centroids from higher to lower energy: [69]. Although in the calculations based on the relativistic Hartree approach (without Fock term in the mean field) this hierarchy is not always reproduced giving a higher value for the E(1−), the overall SDR distribution agrees with the data very reasonably because of the dominance of the 2 − component [56]. Recently, the latter model has been applied to calculations of the SDR in neutron-rich nuclei, such as 132 Sn and 78 Ni [58]. It has been found that in these nuclei in the 2 − component of the SDR a low-energy structure is developed, resembling the well-known low-lying collective quadrupole excitation in the isoscalar channel. Here we consider the case of 132 Sn in more detail and in comparison to a neutron-deficient 100 Sn. Figure 1 displays both the S − and S + branches of the SDR in the latter nuclei. The 0 − , 1 − , 2 − components calculated within the pn-RTBA approach of Ref. [56] are represented by the violet, orange and green curves, respectively. The black solid curve shows the sum of all three components and the dashed curve shows the total SDR strength calculated in pn-RRPA, where the particle-vibration coupling mechanism is absent. Thus, the difference between the solid and dashed black curves shows the effect of the SDR fragmentation due to the PVC. Another general observation is that in 100 Sn both branches are well developed and, although not exactly symmetric because of the presence of the Coulomb interaction, exhibit similar strength distributions, while for 132 Sn very little strength is found in the S + branch, as compared to S − because of the large isospin asymmetry in this nucleus.
Although the PVC effect is quite noticeable for the overall picture, we are interested here in the properties of the low-lying isospin-flip modes. One can see that in both 100 Sn and 132 Sn strong 2 − states form on the lowenergy shoulders of the S − strength distributions. They can bee seen already in the pn-RRPA approach and they are not affected very much by the PVC. In 100 Sn such a state can be seen at about 25 MeV in the S − and also around 0 MeV in the S + branch. In 132 Sn such a state in S − is located at -1.7 MeV and it is visibly separated from the giant resonance. Notice that all the energies are given with respect to the ground state of the 100,132 Sn nuclei and, of course, become positive with respect to the daughter nuclei after the corresponding shifts. Here and in the following we keep the energy scale relative to the ground states of the tin nuclei, in particular, in the self-energies considered below in Eqs. (4), (5).
Thus, there appears another aspect of the isospin transfer modes which has been barely discussed in the literature. As low-energy collective excitations can appear in the isovector channel, similarly to the ones in the isoscalar channel, they are likely to couple to singleparticle motion. Thus, such modes should be included in the nucleonic self-energy and their contributions should be investigated numerically [35]. The next sections are devoted to the implementation of this task.

C. Nucleonic self-energy with pion exchange
As in Refs. [48,50] and many other applications, here we use the diagonal approximation for the Eqs. (1), (2). In a spherical non-superfluid system, the proper neutron and proton self-energies Σ where γ ηµ;ηnηn′ and where the round brackets indicate the full sets of quantum numbers with excluded magnetic ones (reduced matrix elements), and the analogous convention for the matrix elements with interchanged proton and neutron indices applies. In Eqs. (4), (5) we distinguish between the isoscalar phonons with the vertices γ µ and frequencies Ω µ and isovector (proton-neutron (pn) and neutron-proton (np)) ones with the vertices ζ λ and frequencies ω λ . These quantities are to be extracted from the response function of various multipolarities J µ , J λ and parities. For the isoscalar phonons this procedure is described in detail in Refs. [48,52] and consists of the RRPA calculations for the particle-hole (ph) and hole-particle (hp) components of the vertices γ µ and subsequent determination of the particle-particle (pp) and hole-hole (hh) components entering the first sum of the Eqs. (4), (5). The characteristics of the isovector phonons are calculated analogously. First, the proton-neutron RRPA (pn-RRPA) equations are solved for the ph-components of the transition densities R: where the upper indices indicate the location of the proton and neutron states with respect to the Fermi energy (FE): 'p' (particle) is above the FE (η = 1) and 'h' (hole) is below the FE (η = -1). The quantities ε p,n are the mean-field energies for protons and neutrons, respectively, and the interaction in the isovector channelṼ (±) is represented by the exchange of positively/negatively charged π and ρ mesons, respectively: After the solution of the Eq. (8), the isovector phonon vertices are obtained as follows: and then their reduced matrix elements and frequencies obtained from Eq. (8) are used in the self-energies (4), (5). Their diagrammatic representation is given in the Fig. 2. For both proton and neutron self-energies, the first two terms on the right-hand side represent the contributions from isoscalar vibrations (incoming, outgoing and intermediate single-particle states are of the same isospin). The last two terms are the contributions from isospin-flip vibrations (the intermediate states have different isospin).

III. SINGLE-PARTICLE STATES IN 100,132 SN: DETAILS OF CALCULATIONS, RESULTS AND DISCUSSION
The calculations are performed for 100,132 Sn and the procedure is divided into the three steps. (i) RRPA [70] and pn-RRPA (8) calculations with the parameter set NL3* [71], with the interaction of Eq. (9) and g ′ = 0.6 are done for determining the frequencies and the reduced transition probabilities of the phonon modes. The modes below 15 MeV for 132 Sn and below 20 MeV for 100 Sn (with respect to the daughter nuclei in the isospin-flip case) with J π µ = 2 + , 3 − , 4 + , 5 − , 6 + in the isoscalar channel (only natural parities are included in the phonon space because spin-flip phonons are known to contribute very little [34,35]) and with J π λ = 0 ± , 1 ± , 2 ± , 3 ± , 4 ± , 5 ± , 6 ± in the isovector channel are included in the model space. After that the phonons are selected according to the usual truncation scheme, keeping the modes with the transition probabilities larger than 5% of the one with the maximal probability. The transition probabilities of the spin-isospin-flip phonons (unnatural parity) are evaluated for the operators of Eq.
(3) with L = λ − 1 for λ > 0 and with L = λ + 1 for λ = 0. For the modes with natural parity the usual multipole operators are assumed. (ii) The complementary vertex matrix elements were calculated according to Eq. (10) for the selected isospin-flip phonons, and by a similar procedure [48] for the isoscalar ones. (iii) The obtained frequencies and vertices are included in the neutron and proton self-energies (4), (5) and the Dyson equation (2) is solved in the diagonal approximation. It is transformed to a diagonalization problem, by the method described in Ref. [25]. The first illustrative calculations presented below neglect the terms in the self-energies (4), (5) with η p ′ = η n , η n ′ = η p , which correspond to the backward-going diagrams and contain effects of the ground states correlations (GSC). These GSC associated with the isovector phonons are represented by the last terms in the diagrammatic form of the self-energies Σ (e) pp ′ and Σ (e) nn ′ in Fig. 2. The GSC-terms associated with the isoscalar phonons and given by the second graphs on the r.h.s. of Fig. 2 are included in the calculations.
The results of the calculations for the single-particle states in 100,132 Sn are displayed in Fig. 3. The left columns (NL3*) in each panel correspond to the meanfield case, the second left columns (NL3*+PVC) show the RPVC calculations including only isoscalar phonons, the second right columns (NL3*+PVC+π-dyn) contain the results obtained with the additional contribution of the isovector phonons, and 'experimental' single-particle energies extrapolated from data [72,73] are presented in the right columns (Exp). Notice that the two middle columns show only the dominant single-particle states, i.e. those with the maximal spectroscopic amplitudes.
The influence of the isoscalar phonons on the singleparticle spectra was investigated systematically within  the RPVC model in Ref. [49] and a significant overall improvement of the description of the dominant single-particle states was demonstrated. The columns (NL3*+PVC) correspond to this approach and are to be compared to the next columns (NL3*+PVC+π-dyn) in order to reveal the dynamical contribution of pions to the positions of the dominant single-particle states. The contributions of the isospin-flip phonons are associated with mainly pionic processes, because the analysis of the matrix elements of the interaction (Ṽ ρ +Ṽ π +Ṽ δπ ) in Ref. [74] shows that the contribution from the ρ-meson exchange is relatively small. Thus, as can be seen from Fig. 3, the contribution of the pion dynamics provides additional shifts of the dominant single-particle states, which ranges from a few hundred keV to 1 MeV. The spectroscopic factors for the majority of the considered states change only little compared to those of Ref. [49]. However, for the states far from the Fermi surface like, for instance, the proton state 1f 7/2 in 132 Sn, the effect can look really dramatic, because of a redistribution of the strength and a change of the dominant fragment.
It follows from Eqs. (4,5) that the effect of various contributions to the nucleonic self-energies is determined by the matrix elements γ ηnη n ′ (µ;nn ′ ) , γ ηpη p ′ (µ;pp ′ ) and ζ ηnη p ′ (λ;np ′ ) , ζ ηpη n ′ (λ;pn ′ ) given by Eq. (7), and by the energy denominators in Eqs. (4,5), which include combinations of single-particle mean-field energies ε p ′ , ε n ′ and phonon frequencies: Ω µ for the isoscalar phonons and ω λ for the isovector ones. Although the self-energies of Eqs. (4,5) do not directly provide the splittings and shifts of the single-particle states and the Dyson equation (2) has to be solved for this, in the second-order perturbation theory the matrix elements of Σ (e) k can be estimated by setting ε = ε k in Eqs. (4,5), where the index k stands for both neutron and proton single-particle states. Then it is clear, that the largest contributions are given by the terms with the largest phonon vertices and the lowest phonon frequencies, that determines the above mentioned truncation.
The typical values of the γ µ are between 0.1 MeV and 1 MeV, in rare cases achieving few MeV. For the vertices ζ λ the pn-RRPA gives somewhat smaller (factor 2-3) values in average, while in both isoscalar and isovector cases the largest matrix elements are obtained for the lowest multipoles. The existence of the isoscalar collective phonons with frequencies below 10 MeV is the most obvious reason for the strong impact of PVC on the single-particle spectra (see numerous papers cited in the introduction). In Fig. 3 this effect on the dominant states can be seen by comparison of the first left and the second left columns of each panel. The typical frequencies of isovector spinflip phonons of the lowest multipoles are seen in Fig.  1 which shows that they can be quite low (in particular, 2 − states shown by green curves). The second-order self-energy denominators corresponding to isoscalar and isovector phonons have values of the same order of magnitude. In the latter case the phonon frequencies ω λ add up to proton-neutron single-particle energy differences ∆ pn ′ = ε p − ε n ′ , ∆ np ′ = ε n − ε p ′ , producing the energy denominators D pn ′ λ = ∆ pn ′ ∓ ω λ , D np ′ λ = ∆ np ′ ∓ ω λ . Among many terms with large denominators of several tens MeV, there are usually a few of them with the energies below 10 MeV, which give the leading contributions for both IS and IV phonons. Figures 4, 5, 6 illustrate a quantitative assessment of the contributions of isovector phonons with various multipolarities and frequencies to the single-particle characteristics. Those figures show single-particle strength distributions for the two neutron states of 132 Sn: 2p 1/2 (Figs. 4, 5) and 3s 1/2 (Fig. 6). The latter represents a so-called "good single-particle state" because of its closeness to the Fermi edge, and the former is an example of a "deep hole state", which is far from it.  included in addition to the IS ones, in comparison to the initial mean-field state (dashed blue). One can see that the 2p 1/2 state is already fragmented considerably by coupling to IS phonons, and two major fragments can be distinguished at - 15   with J π λ = 0 ± , 1 ± (b), J π λ = 0 ± , 1 ± , 2 ± , 3 ± (c) and J π λ = 0 ± − 6 ± (d) are included up to 15 MeV, compared to the initial mean-field state (dashed blue). These figures illustrate the role of various multipoles of the IV phonons in the single-particle strength distribution. Similarly to the case of IS phonons, known from literature, low-J modes give the major contribution to the strength redistribution while the phonons with J π λ = 4 ± , 5 ± , 6 ± produce a weaker effect, that justifies the truncation of the phonon basis at J λ = 6.
The change of the 3s 1/2 state because of the IV-type PVC, which is shown in Fig. 6, reduces to a gradual move toward the Fermi surface with every extension of the phonon space and a very little fragmentation. The effect also becomes smaller when the phonon energy increases while no clear saturation is seen for this state for the IV phonons up to 20 MeV. This points out that a proper regularization may be needed for the considered mass operator, for instance, a subtraction technique analogous to the one employed in the response theory [52] can be adopted.
Another consequence coming from the analysis above is that there is no clear dominance of the lowest IV phonons. This emphasizes the distinction between the energies of the most important IS modes and those of the of IV modes, due to the fact that the relevant energy differences between neutron or proton states are typically smaller than those between proton-neutron states.
In conclusion, except for a few cases, the shifts of the dominant states in Fig. 3 are toward the Fermi surface and give an additional improvement of the description, compared to available data. Proton particle states in 132 Sn seem to be very sensitive to the PVC effects, which tend to overestimate the shifts, however, the obtained picture can change further after the inclusion of ground state correlations associated with the IV phonons, which will be addressed by future efforts.

IV. SUMMARY AND OUTLOOK
The dynamical contribution of the pion is included into the single-particle self-energy of a finite nuclear system in a self-consistent relativistic framework. It has been shown that in medium-mass nuclei low-lying isospin-flip states associated with soft pionic modes can occur with sizable transition probabilities and, therefore, their coupling to single-particle and other collective modes should be considered. The strength of this coupling and its effect on the single-particle states is evaluated for 100,132 Sn and the shell structure of these nuclei is found to be sensitive to pion dynamics. Thereby, (i) pion-nucleon correlations beyond the Hartree-Fock approximation are included in the theory based on the QHD Lagrangian; (ii) the latter can provide a link between the relativistic nuclear field theory and chiral effective field theories; (iii) the relativistic particle-vibration coupling model is extended self-consistently by incorporating isospin-flip phonons.
As the nucleonic self-energy of the presented approach contains a full summation of the ring diagrams with pion exchange, the pion contribution is included non-perturbatively, although its soft modes entering the nucleonic self-energy are obtained by calculations on the RPA level. The next natural step will be an inclusion of the ground state correlations associated with the exchange of the isospin-flip phonons, which will provide a more accurate description of the effects considered in the present work. Pairing correlations of the superfluid type should be taken into account for the isovector phonons, which will allow for more systematic studies of the effects of pion dynamics in open-shell nuclei. The considered effects of pion dynamics should have also an impact on nuclear response in various channels, therefore, a corresponding extension of the response theory can be another exciting future development.

V. ACKNOWLEDGEMENT
Enlightening discussions with E.E. Kolomeitsev, P. Ring, V. Tselyaev, T. Otsuka, and V. Zelevinsky are gratefully acknowledged. The author is very thankful to T. Marketin for providing a part of the code for pn-RRPA matrix elements. This work was supported by US-NSF grants PHY-1204486 and PHY-1404343.