
DOI: 10.22661/AAPPSBL.2020.30.1.67
A BeyondMeanField Description for Nuclear Excitation Spectra:
The Subtracted Second Randomphase Approximation
MARCELLA GRASSO^{1} ^{1}LABORATOIRE DE PHYSIQUE DES 2 INFINIS IRÈNE JOLIOT CURIE,
UNIVERSITÉ PARISSACLAY, 91406, ORSAY, FRANCE
communicated by Yongseok Oh
This is a review of recent results that were obtained by applying the subtracted second random phaseapproximation (SSRPA) method. This theoretical model is an extension of the randomphase approximation (RPA) where two particletwo hole configurations are also included, together with the one particleone hole RPA configurations. In addition, a subtraction procedure is applied to remove all the problems related to the overcounting of correlations, to various instabilities, as well as to the occurrence of an ultraviolet divergence generated by the zero range of the used effective interaction.
Recent applications are presented, namely a study of dipole excitations in the nucleus ^{48}Ca, a systematic analysis of giant quadrupole excitations in several mediummass and heavy nuclei, an estimate of beyondmeanfield effects on the effective mass of nuclear matter, and, finally, the study of soft monopole compression modes in the neutronrich nucleus ^{68}Ni. Such excitations are found to be neutrondriven and a link with a compressibility modulus defined for isospin asymmetric nuclear matter is discussed.
PACS numbers: 21.60.Jz, 21.10.Re, 27.20.+n, 27.40.+z
INTRODUCTION
Nuclear excitation spectra are described in a quite satisfactory way by employing the RPA model, based on individual degrees of freedom, where excited states are constructed as superpositions of singleparticle particlehole (ph) configurations. This leads in general to a fair description of centroid energies for giant resonances. However, by construction, the fragmentation and the spreading of such excitations around their centroids can be accounted for only through the socalled Landau damping, which is related to the structure of singleparticle spectra. In other words, only a meanfieldbased description can be provided and correlations beyond the meanfield approximation, which could generate an additional fragmentation in the excitation spectra, are totally neglected.
This is the main motivation for going beyond the meanfield approximation in the description of excited modes, both lowlying and collective giant resonances. Along this direction, we have performed such an extension by employing the second RPA (SRPA) model. The SRPA model represents an intuitive and natural extension of the RPA model, done by including also two particletwo hole (2ph) configurations in the definition of the excitation operators, in addition to the RPA ph configurations. Such a mixing between individual degrees of freedom and multiparticlemultihole configurations provides a way to account for a beyondmeanfield description of the spreading width and the fragmentation of excited states.
The formal derivation of the SRPA model was developed several decades ago using different methods such as the equationsofmotion method [1], a variational procedure [2], or the smallamplitude limit of the time dependent density matrix [3, 4]. Numerous applications have been published over the years. However, owing to the extremely heavy numerical effort which is required (in practice, the diagonalization of a very large matrix), strong cuts and approximations were employed such as using very low cutoff values for 2ph configurations, making a diagonal approximation in the block of the matrix where 2ph configurations are coupled among themselves, or simply relying on a second TammDancoff model.
In the last decade, first applications could be done in nuclear physics by two groups without such strong truncations. One of the groups used realistic microscopic interactions derived from Argone V18 [5, 6], whereas the second group (ours) conceived the first largescale SRPA calculations with phenomenological effective interactions of Skyrme and Gogny type [710]. After the very first applications, attention was also paid to the definition of the socalled rearrangement terms of the residual interaction in the context of SRPA, in cases where density dependent effective interactions are used. To address this problem, a method was proposed which is an extension of the previously mentioned variational procedure introduced by Providencia [8].
These first largescale applications could allow us to identify several formal problems and limitations which characterize the SRPA theoretical framework. Some of these problems are indeed quite general and are encountered in all the applications of the SRPA model, whatever is the interaction or the treated manybody system. Such problems are linked to the fact that the Thouless theorem cannot be extended in an automatic way to the SRPA framework. Thus, even if the HartreeFock (HF) state corresponds to the state which minimizes the HF Hamiltonian, there is no guarantee that the SRPA stability matrix is positive semidefinite. This may generate instabilities such as imaginary eigenvalues or eigenvalues with positive (negative) energy and negative (positive) norm. Related to such a limitation, another general problem of the SRPA model in its standard form exists, which is an anomalous (and unphysical) global shift of the energy spectrum towards lower energies, found compared to the RPA spectrum. In the nuclear case, such a shift may be of several MeV. These instabilities were analyzed in the last years for example by Tselayev [11] and Papakonstantinou [12].
In addition, the SRPA model suffers from other two limitations which occur when the model is formulated within the socalled energydensityfunctional (EDF) framework. First, a risk of double counting of correlations exists. In most cases, the parameters of the used density functionals are adjusted to reproduce specific observables by performing meanfield calculations. This implies that, to some extent, some correlations are contained in the numerical values of the adjusted parameters. Thus, if correlations are explicitly introduced when the meanfield approximation is overcome, a risk of double counting of correlations may exist. In addition, if the used functionals are derived from zerorange interactions, an ultraviolet divergence occurs in the beyondmeanfield SRPA model.
Following the suggestion of Tselayev [11], we have implemented the SRPA model using a subtraction procedure that was introduced by Tselayev in 2007 [13] to account for the problem of double counting of correlations within particlevibrationcoupling models. Tselyaev demonstrated in Ref. [11] that the use of such a prescription allows one to enlarge the Thouless theorem to be valid for extensions of RPA such as the SRPA model. We have then formulated a subtracted SRPA (SSRPA) model based on both Skyrme [14] and Gogny [15] effective interactions and seen that, using this procedure, also the ultraviolet divergences produced by contact interactions are canceled.
After a brief review of the formal model (Sec. II), I will resume in Sec. III some recent applications that we have performed, namely the study of dipole excitations and of the electric dipole polarizability in ^{48}Ca, a systematic analysis of quadrupole excitations in a series of nuclei, from mediummass systems to heavy nuclei, an estimation of a beyondmeanfield effect on the effective mass, and the study of soft monopole modes in neutron rich nuclei (link with the incompressibility modulus of asymmetric nuclear matter). Section IV illustrates the main conclusions of this manuscript.
A FEW WORDS ON THE EMPLOYED FORMAL MODEL
Details on the SRPA and the SSRPA equations may be found in several articles quoted in the Introduction. The compact form of the SRPA equations is
 (1) 
where, as in the RPA case, the solutions are the excitation energies ω and the amplitudes (X, Y) which characterize the wave functions of the corresponding excited states. The structure of the matrix is obviously much more complicated than in RPA. For instance, the matrix A has the following form:
where the indices '1' and '2' indicate ph and 2ph configurations, respectively. With this notation, A_{11'} is the standard RPA matrix. It is possible to write the SRPA equations as RPAlike equations with energydependent matrix elements. In this way, the RPAlike A matrix has the following form:
 (2) 
where the first term is the RPA matrix A and the following terms represent the energydependent selfenergy insertion. Within the subtraction procedure, the energy dependent matrix element is replaced by a subtracted one having the form
 (3) 
where E_{11'}(ω) is the energydependent part (selfenergy) of A_{11'}(ω) in Eq. (2).
This subtraction is done for the following reasons. Within a densityfunctionaltheorylike spirit, one assumes that the RPA response corresponds to the exact static response. To avoid any double counting of correlations when the meanfield approximation is overcome, one has to require that the static (zeroenergy) response is not modified with respect to the RPA response. For this, one imposes that the static polarizability is conserved (equal to the RPA one). One way to guarantee this is to perform the subtraction which is described in Eq. (3).
SOME RECENT APPLICATIONS
In all the recent applications that I will present in this Section, the calculations were performed in a fully self consistent way, in the sense that the residual interaction contains all the terms corresponding to the Skyrme interaction. Several Skyrme parametrizations were used.
The dipole strength, both lowlying and in the giant resonance region, was studied in Ref. [16] for the nucleus ^{48}Ca. It was first shown, by comparing with the corresponding experimental results [17], that the subtraction procedure considerably improves the lowlying response (between 5 and 10 MeV), both in the distribution of the main energy peaks and in the prediction of the transition probability, compared to previous SRPA calculations, where the subtraction procedure was not yet employed [9]. One has to mention that the RPA model does not provide any strength between 5 and 10 MeV for this nucleus. The prediction for the giant dipole resonance was compared with more recent experimental results [18]. It was shown that the centroid energy was slightly underestimated compared to the experimental distribution (for instance by 1.5 MeV with the Skyrme parametrization SGII). However, the important result of this work was that the width of the excitation mode was indeed very well reproduced.
The corresponding electric dipole polarizability was also measured and Fig. 1 shows the comparison with several theoretical predictions. The dipole polarizability is calculated using the expression
 (4) 
where B(E1) is the dipole transition probability. RPA and SSRPA results were obtained with the Skyrme parametrization SGII. The RPA curve indicates that the centroid is quite well reproduced. However, the increase of the polarizability is much sharper than the experimental result in the region of the giant resonance, indicating that the width of such an excitation is largely underestimated. The abinitio band is a collection of coupledcluster results. The figure shows that abinitio predictions overestimate the centroid of the giant resonance. In addition, they provide a quite sharp increase of the polarizability, implying that the width is too small compared to the experimental value. The SSRPA prediction, on the other side, provides a slightly underestimated centroid (1.5 MeV with the used parametrization SGII). The most important finding is that the slope of the curve is in very good agreement with the experimental one. This may be more easily seen in the figure by looking at the green squares, which represent SSRPA results after a global shift of 1.5 MeV.
Fig. 1: Electric dipole polarizability versus the excitation energy. The grey band corresponds to abinitio results and the region between the two solid red lines describes the experimental results reported in Ref. [18]. RPA (blue dotted line) and SSRPA (black dashed line) results are shown. The green squares and dashed line represent SSRPA results shifted upwards by 1.5 MeV. The Skyrme interaction SGII is used. Adapted from Ref [16].
The first conclusion from this study was that the SSRPA model produces a very good agreement with the experimental measurement in the description of the lowlying spectrum of ^{48}Ca. At the same time, such a model is able to describe very well the fragmentation and the width of the giant dipole resonance for the same nucleus.
We then moved to the description of isoscalar giant quadrupole resonances and performed a systematic analysis of their centroids and widths, exploring the region of mediummass and heavy nuclei [19]. Figure 2 shows for example our results for the centroids. Experimental results are displayed for comparison. They were extracted from Ref. [20]. RPA and SSRPA predictions can be compared with the experimental results for those nuclei which are indicated in red (vertical dotted red lines). One can clearly observe that the SSRPA description is in overall better agreement with the experimental centroids, compared to the RPA case. For the widths, we have found that, systematically, the SSRPA model provides larger widths than the RPA model, as obviously expected.
To give an example, I show in Fig. 3 the quadrupole response predicted in RPA (a) and in SSRPA (c) for the nucleus ^{120}Sn. The used interaction is SLy4. Panel (c) shows both the discrete SSRPA spectrum and the spectrum obtained by folding with a very narrow Lorentzian of 40 keV width, which corresponds to the experimental accuracy. It is easy to see how the SSRPA naturally describes the width of excited states, so that the folding does not add any artificial spreading effect. The corresponding highresolution experimental measurement [21] is shown in panel (b). One observes the strong improvement that the SSRPA model produces with respect to RPA in the description of the fragmentation. The RPA spectrum is composed in practice by a unique energy main peak.
Fig. 2: Centroid energies of the isoscalar GQR. The experimental data, shown as black circles, are extracted from Ref. [20]. SSRPA (RPA) predictions are plotted as green diamonds (violet circles). At A = 48, there are two experimental measurements, corresponding to ^{48}Ti and for ^{48}Ca. The point for ^{48}Ca is the highest one. Theoretical calculations are carried out for ^{30}Si, ^{34}Si, ^{36}S, ^{40}Ca, ^{48}Ca, ^{56}Ni, ^{68}Ni, ^{90}Zr, ^{114}Sn, ^{116}Sn, ^{120}Sn, ^{132}Sn, and ^{208}Pb with the Skyrme parametrization SLy4. Adapted from Ref. [19].
Fig. 3: Isoscalar giant quadrupole response function for the nucleus ^{120}Sn. RPA and SSRPA results are shown on panels (a) and (c), respectively. Panel (b) displays the experimental measurement of Ref. [21]. For RPA and SSRPA discrete spectra, units are e^{2} fm^{4}. For the SSRPA folded case, units are e^{2} fm^{4} MeV^{1}. Adapted from Ref [19].
Such a study on isoscalar giant quadrupole resonances led also to another interesting analysis. A very well known and used relation exists between the effective mass m^{*} of infinite matter and the centroid energy of isoscalar giant quadrupole resonances in nuclei. Such a relation is very general for a manybody system and may be derived within the Landau theory of Fermi liquids using the localdensity approximation [22, 23]. In the nuclear case, the centroid energy of isoscalar giant quadrupole resonances is proportional to the square root of m / m^{*}, where the effective mass is not the effective mass in nuclei but in nuclear matter. The definition of the effective mass of a particle with energy E and momentum k is
 (5) 
where
 (6) 
One recognizes in the above equation the selfenergy Σ_{k}+Σ_{k,E}, where the first term is the meanfield contribution Σ_{k} and the second term describes the beyondmean field energydependent contribution Σ_{k,E}. The mean field contribution is different from 1 only if the effective interaction has an explicit momentum dependence (gradient terms). An energy dependence may be acquired only beyond the meanfield approximation and, for this reason, the beyondmeanfield contribution to the effective mass is equal to 1 if genuine meanfield calculations are carried out. One can derive the following expression:
 (7) 
where the socalled Emass m^{*}_{E} /m and kmass m^{*}_{k} /m have been introduced.
Using the known relation between the centroid energies of isoscalar giant quadrupole resonances and the effective mass of infinite matter, we proposed a methodology for extracting and estimating a beyondmeanfield effect on the effective mass of matter [24]. Figure 4 shows RPA centroids obtained for the nuclei ^{48}Ca and ^{90}Zr with the four Skyrme parametrizations SkP, SGII, SLy4, and Ska, which have quite different values for the effective mass, ranging from 0.6 to 1 (in the meanfield approximation). The RPA centroids are then associated with the corresponding values of m^{*} / m calculated with the same Skyrme interaction for infinite matter in the meanfield scheme. The above mentioned linear relation is clearly visible. A linear fit is performed (blue dotted lines). On these blue dotted lines, we have reported the SSRPA centroids obtained for the two nuclei and the two interactions SGII and SLy4. One can in this way extract, for each of these four points, an estimate for a beyondmeanfield effective mass (which is also nucleus dependent, and not only interaction dependent as it was for the meanfield case).
To make a more global analysis, I show in Fig. 5, panel (a), the meanfield effective masses for infinite matter associated with three interactions SGII, SLy4, and Ska. The colored bands describe the associated theoretical uncertainty (orange area for the dispersion between SGII and SLy4 and orange plus grey area for the dispersion between SGII and Ska). It is known that typical mean field values for the effective masses are located in the window between 0.6 and 0.8, which implies that an uncertainty band of 23% represents a typical theoretical dispersion for the meanfield case. Let us now go to the beyondmeanfield case, illustrated in panel (b). First of all, one may observe the global shift upwards of the colored bands, indicating the enhancement of the effective mass induced beyond the mean field. The extraction procedure of the beyondmeanfield effective mass was applied also to the nucleus ^{120}Sn, in addition to the nuclei ^{48}Ca and ^{90}Zr. The values which were obtained for the three Skyrme parametrizations are reported. One may observe that it is possible to define also in this case uncertainty bands associated with the three used interactions. It is interesting to notice that the uncertainty bands remain comparable to the meanfield case, which means that the additional nucleus dependence which is introduced by such a specific methodology to extract the effective mass does not generate an increase of the theoretical uncertainty band. We have obtained that the E mass got an increase (from the meanfield value, which is 1) between 6 and 16%, according to the each specific case (nucleus and interaction).
Fig. 4: Centroids for the isoscalar giant quadrupole resonances for ^{48}Ca and ^{90}Zr as a function of _{}. Four RPA centroids (black squares) are shown, obtained with four Skyrme parametrizations. The corresponding meanfield effective masses in nuclear matter are associated with them. The blue dotted lines represent the linear fits that were performed on such points. The SSRPASLy4 and SSRPASGII centroids are shown on the blue dotted lines (green triangles and magenta squares, respectively). The grey bands illustrate the experimental values for the centroids in the two nuclei. Adapted from Ref. [24].
Fig. 5: (a) Meanfield effective masses for the interactions SGII, SLy4, and Ska; (b) Beyondmeanfield effective masses for the same interactions as in (a) and for the nuclei for ^{48}Ca, ^{90}Zr, and ^{120}Sn. Adapted from Ref. [24].
Fig. 6: Monopole isoscalar response for ^{68}Ni. Adapted from Ref. [25]
Fig. 7: Neutron and proton transition densities multiplied by r^{2} (in units of fm^{1}) for the peaks located at 11.2, 12.53, and 12.92 MeV in the monopole response of ^{68}Ni. Adapted from Ref. [25].
More recently, we have analyzed lowenergy or soft monopole compression modes in several neutronrich nuclei, among them ^{68}Ni [25]. The monopole spectrum, covering both the lowenergy part and the giant monopole resonance, is shown in Fig. 6. The used interaction is SGII. We have found several lowenergy peaks and the transition densities associated with three of them are plotted on Fig. 7. Such transition densities reveal the nature of these soft excitations. It can be clearly seen that these lowenergy compression modes are not breathing collective modes where neutrons and protons participate coherently. The excitations are mainly driven by neutrons in the whole volume of the nucleus, both in its internal region and at the surface.
Now, it is known that giant monopole compression modes may be linked to the compressibility modulus of symmetric infinite nuclear matter. The compressibility for symmetric nuclear matter (having an associated equation of state E^{sym}/A) is defined as
 (8) 
where ρ_{0} is the equilibrium density of symmetric matter, equal to 0.16 fm^{3} (accepted empirical value for the saturation density of symmetric matter). If the nucleus is described as a liquid drop with a compressibility modulus equal to K, the following relation between the centroid energy of a giant monopole resonance and the compressibility modulus may be derived:
 (9) 
In the above expression η_{0} is the root mean square radius. If for instance η_{0} = r_{0 }A^{1/3}, with r_{0} ~ 1 fm, one can write
 (10) 
Such a very simple and general relation generates quite reasonable estimations for the centroid energies of isoscalar giant monopole resonances, knowing the compressibility modulus produced by each specific effective interaction.
Now, we have seen that soft monopole modes are neutrondriven excitations. We have for this reason proposed a method to relate the excitation energies of these modes to a compressibility modulus defined for asymmetric nuclear matter. Let us introduce an isospin asymmetry parameter X, which will describe the isospin asymmetry of nuclear matter as well as the isospin asymmetry of the neutrondriven excitation mode that we want to describe: X =(ρ_{n}  ρ_{p}) / ρ, where ρ_{n}, ρ_{p}, and ρ are neutron, proton and total densities, respectively. From the SSRPA solutions, computing neutron
ρ_{ν}^{n}
and proton ρ_{ν}^{p} transition densities for a given energy peak ν, one can define the following quantities
 (11) 
and
 (12) 
From this, one can estimate the isospin asymmetry for a given excitation as
 (13) 
This will provide a microscopic prediction for the value of X associated with a given excited mode.
We can also define a compressibility modulus for asymmetric matter with asymmetry X,
 (14) 
which is defined exactly as the compressibility for symmetric matter, where now the equation of state describes asymmetric matter and ρ_{eq} is the equilibrium density associated with this equation of state. Such equations of state are plotted on Fig. 8 for the interaction SGII, where the green circles describe the equilibrium point of each equation of state. In the legend of the figure one can find the compressibility modulus associated with each equation of state.
Fig. 8: Equations of state computed with the parametrization SGII, from symmetric to pure neutron matter. Dotted lines represents the equations of state of asymmetric matter with asymmetry values X (from the bottom to the top) equal to 0.1 (K_{X} = 211.59 MeV), 0.2 (K_{X} = 203.54 MeV), 0.3 (K_{X} = 187.33 MeV), 0.4 (K_{X} = 166.40 MeV), 0.5 (K_{X} = 139.90 MeV), 0.6 (K_{X} = 108.30 MeV), 0.7 (K_{X} = 72.35 MeV), 0.8 (K_{X} = 33.74 MeV), 0.85 (K_{X} = 13.32 MeV). The compressibility value associated with symmetric matter for the parametrization SGII is also indicated in the figure. The equilibrium points (at which the compressibility values K_{X} are calculated) are represented by green circles for each equation of state. Adapted from Ref. [25].
With similar arguments as those used for the case of the giant monopole resonance and the compressibility of symmetric matter, one can introduce the following relation,
 (15) 
where now the energy refers to a soft neutrondriven mode and the compressibility to asymmetric nuclear matter. This equation is plotted for three values of A, 34, 60, and 68 on Fig. 9. Using these curves, one can extract in a very simple way an estimation of the isospin asymmetry of a soft mode located at a given energy. For instance, one can extract the value of X associated with the excitation located at 11.02 MeV in ^{68}Ni: a value of X between 0.7 and 0.8 is obtained in this way. If now one computes X microscopically by using Eq. (13), one finds the value of 0.78 for the excitation energy of 11.02 MeV, in agreement with the previous qualitative estimation extracted by using the compressibility of asymmetric matter.
Fig. 9: Trends provided by Eq. (15) for three values of A, 34, 60, and 68.
CONCLUSIONS
This breaf review contains and illustrates several recent results that have been obtained within the SSRPA model. After a short introduction of this theoretical model, several applications are presented. Dipole excitations and the electric dipole polarisability are analyzed for the nucleus ^{48}Ca. A systematic analysis of isoscalar giant quadrupole resonances, carried out for several mediummass and heavy nuclei, is illustrated. Related to giant quadrupole energies, an estimation of beyondmeanfield effects on the effective mass of infinite nuclear matter is discussed. Finally, soft neutron driven compression modes are studied and the link with a compressibility modulus defined for isospin asymmetric matter is explored.
Acknowledgments: I acknowledge the warm hospitality of the Korean colleagues who invited me at the APCTP Focus Program 'Nuclear ManyBody Theories: Beyond the meanfield approaches', which took place in July 2019 in Pohang.
I acknowledge the work done together with my collaborators. In particular, for the applications illustrated here, I acknowledge the work done in collaboration with Danilo Gambacurta, Olivier Sorlin, and Olivier Vasseur.
References
[1] C. Yannouleas, Phys. Rev. C 35, 1159 (1987).
[3] M. Tohyama, Gong, Z. Phys.A 332, 269 (1989).
[2] J. da Providencia, Nucl. Phys. 61, 87 (1965).
[4] D. Lacroix, S. Ayik, Ph. Chomaz, Prog. Part. Nucl. Phys. 52, 497 (2004).
[5] P. Papakonstantinou and R. Roth, Phys. Lett. B 671, 356 (2009).
[6] P. Papakonstantinou and R. Roth, Phys. Rev. C 81, 024317 (2010).
[7] D. Gambacurta, M. Grasso, and F. Catara, Phys. Rev. C 81, 054312 (2010).
[8] D. Gambacurta, M. Grasso, and F. Catara, J. Phys. G 38, 035103 (2011).
[9] D. Gambacurta, M. Grasso, and F. Catara, Phys. Rev. C 84, 034301 (2011).
[10] D. Gambacurta, M. Grasso, V. De Donno, G. Co', and F. Catara, Phys. Rev. C 86, 021304(R) (2012).
[11] V.I. Tselyaev, Phys. Rev. C 88, 054301 (2013).
[12] P. Papakonstantinou, Phys. Rev. C 90, 024305 (2014).
[13] V.I. Tselyaev, Phys. Rev. C 75, 024306 (2007).
[14] D. Gambacurta, M. Grasso, and J. Engel, PRC 92, 034303 (2015).
[15] D. Gambacurta and M. Grasso, Eur. Phys. J. A 52, 198 (2016).
[16] D. Gambacurta, M. Grasso, and O. Vasseur, Phys. Lett. B 777, 163 (2018).
[17] T. Hartmann et al., Phys. Rev. Lett. 93, 192501 (2004).
[18] J. Birkhan et al., Phys. Rev. Lett. 118, 252501 (2017).
[19] O. Vasseur, D. Gambacurta, and M. Grasso, Phys. Rev. C 98, 044313 (2018).
[20] G. Scamps and D. Lacroix, Phys. Rev. C 88, 044310 (2013).
[21] A. Shevchenko et al, Phys. Rev. Lett. 93, 122501 (2004).
[22] J.P. Blaizot, Phys. Rep. 64, 171 (1980).
[23] BaoAn Li et al., Prog. Part. Nucl. Phys. 99, 29 (2018).
[24] M. Grasso, D. Gambacurta, and O. Vasseur, Phys. Rev. C 98, 051303 (R) (2018).
[25] D. Gambacurta, M. Grasso, and O. Sorlin, Phys. Rev. C 100, 014317 (2019).
