Implementation of the finite amplitude method for the relativistic quasiparticle random-phase approximation
Tamara Niksic, Nenad Kralj, Tea Tutis, Dario Vretenar, Peter Ring
aa r X i v : . [ nu c l - t h ] N ov Implementation of the finite amplitude method for the relativisticquasiparticle random-phase approximation
T. Nikˇsi´c , N. Kralj , T. Tutiˇs , D. Vretenar , P. Ring Physics Department, Faculty of Science,University of Zagreb, 10000 Zagreb, Croatia and Physik-Department der Technischen Universit¨at M¨unchen, D-85748 Garching, Germany (Dated: July 11, 2018)
Abstract
A new implementation of the finite amplitude method (FAM) for the solution of the relativis-tic quasiparticle random-phase approximation (RQRPA) is presented, based on the relativisticHartree-Bogoliubov (RHB) model for deformed nuclei. The numerical accuracy and stability ofthe FAM – RQRPA is tested in a calculation of the monopole response of O. As an illustrativeexample, the model is applied to a study of the evolution of monopole strength in the chain of Smisotopes, including the splitting of the giant monopole resonance in axially deformed systems.
PACS numbers: 21.60.Jz, 21.60.Ev . INTRODUCTION Vibrational modes and, more generally, collective degrees of freedom have been a re-curring theme in nuclear structure studies over many decades. An especially interestingtopic that has recently attracted considerable interest is the multipole response of nuclei farfrom stability and the possible occurrence of exotic modes of excitation [1, 2]. For theo-retical studies of collective vibrations in medium-heavy and heavy nuclei the tool of choiceis the random-phase approximation (RPA), or quasiparticle random-phase approximation(QRPA) for open-shell nuclei [3]. The (Q)RPA equations can be obtained by linearizing thetime-dependent Hartree-Fock (Bogoliubov) equations, and the standard method of solutionpresents a generalized eigenvalue problem for the (Q)RPA matrix. As the space of quasi-particle excitations can become very large in open-shell heavy nuclei, the standard matrixsolution of the QRPA equations is often computationally prohibitive, especially for deformednuclei. In addition to the fact that a very large number of matrix elements has to be com-puted, the calculation of each matrix elements is complicated by the fact that most modernimplementations of the (Q)RPA are fully self-consistent, that is, the residual interactionis obtained as a second functional derivative of the nuclear energy density functional withrespect to the nucleonic one-body density. Since energy density functionals can be quitecomplicated and include many terms [4, 5], this produces complex residual interactions andthe computation of huge (Q)RPA matrices becomes excessively time-consuming.Although several new implementations of the fully self-consistent matrix QRPA for axiallydeformed nuclei have been developed in recent years [6–10], and even applied to studies ofcollective modes in rather heavy deformed nuclei, the huge computational cost has so farprevented systematic studies of multipole response in deformed nuclei. An interesting andvery useful alternative solution of the (Q)RPA problem has recently been proposed, basedon the finite-amplitude method (FAM) [11] . In this approach one avoids the computationand diagonalization of the (Q)RPA matrix by calculating, instead, the fields induced bythe external one-body operator and iteratively solving the corresponding linear responseproblem. The FAM for the RPA has very successfully been employed in a self-consistentcalculation of nuclear photo absorption cross sections [12], and in a study of the emergenceof pygmy dipole resonances in nuclei far from stability [13]. More recently the FAM has beenextended to the quasiparticle RPA based on the Skyrme Hartree-Fock-Bogoliubov (HFB)2ramework [14–16]. The feasibility of the finite amplitude method for the relativistic RPAhas been investigated in Ref. [17].In this work we report a new implementation of the FAM for the relativistic quasipar-ticle random-phase approximation (RQRPA), based on the relativistic Hartree-Bogoliubov(RHB) model for mean-field studies of deformed open-shell nuclei [18, 19]. The standard ma-trix RQRPA for spherical nuclei was formulated in the canonical single-nucleon basis of theRHB model [20], extended to the description of charge-exchange excitations (pn-RQRPA) inRef. [21], and further extended to deformed systems with axial symmetry in Ref. [8]. Herewe develop a FAM method for the small-amplitude limit of the time-dependent Hartree-Bogoliubov framework based on relativistic energy density functionals and a pairing forceseparable in momentum space, and perform tests and illustrative calculations for the newmodel.The paper is organized as follows. In Sec. II we briefly recapitulate the small-amplitudelimit of the time-dependent RHB model and present a new implementation of the FAM forthis particular framework. Numerical details and test calculations are included in Sec. III,and in Sec. IV we apply the model to a study of the evolution of monopole strength in thechain of Sm isotopes. Section V summarizes the results and ends with an outlook for futureapplications. Details on the expansion of single-nucleon spinors in the axially symmetricharmonic oscillator basis, calculation of the matrix elements monopole operator, time-oddterms in the FAM equations, and RHB and FAM equations with time-reversal symmetry,are included in Appendix A-D.
II. SMALL AMPLITUDE LIMIT OF THE TIME-DEPENDENT RHB MODELAND THE FINITE AMPLITUDE METHOD
The relativistic Hartee-Bogoliubov (RHB) model [18, 19] provides a unified descriptionof nuclear particle-hole ( ph ) and particle-particle ( pp ) correlations on a mean-field level bycombining two average potentials: the self-consistent nuclear mean field that encloses allthe long range ph correlations, and a pairing field ˆ∆ which sums up the pp -correlations. Inthe RHB framework the nuclear single-reference state is described by a generalized Slaterdeterminant | Φ i that represents a vacuum with respect to independent quasiparticles. Thequasiparticle operators are defined by the unitary Bogoliubov transformation, and the cor-3esponding Hartree-Bogoliubov wave functions U and V are determined by the solution ofthe RHB equation: h D − m − λ ∆ − ∆ ∗ − h ∗ D + m + λ U k V k = E k U k V k . (1)In the relativistic case the self-consistent mean-field is included in the single-nucleon DiracHamiltonian ˆ h D , ∆ is the pairing field, and U and V denote Dirac spinors. In the formalismof supermatrices introduced by Valatin [22], the RHB functions are determined by the Bo-goliubov transformation which relates the original basis of particle creation and annihilationoperators c n , c † n (e.g. an oscillator basis) to the quasiparticle basis α µ , α † µ cc † = W αα † with W = U V ∗ V U ∗ . (2)In this notation a single-particle operator can be represented in the matrix form:ˆ F = 12 (cid:16) α † α (cid:17) F αα † + const . (3)with F = F F F − ( F ) ⊺ . (4)In particular, for the generalized density R : R = ρ κ − κ ∗ − ρ ∗ , (5)where the density matrix and pairing tensor read ρ = V ∗ V ⊺ , κ = V ∗ U ⊺ , (6)respectively, and the RHB Hamiltonian is given by a functional derivative of a given energydensity functional with respect to the generalized density: H = δE [ R ] δ R = h ∆ − ∆ ∗ − h ∗ . (7)4he evolution of the nucleonic density subject to a time-dependent external perturbationˆ F ( t ) is determined by the time-dependent relativistic Hartree Bogolyubov (TDRHB) equa-tion: i∂ t R ( t )= [ H ( R ( t )) + F ( t ) , R ( t )] . (8)For a weak harmonic external fieldˆ F ( t ) = η ( ˆ F ( ω ) e − iωt + ˆ F † ( ω ) e iωt ) , (9)characterized by the small real parameter η , the density undergoes small-amplitude oscil-lations around the equilibrium with the same frequency ω , that is, in the small-amplitudelimit of the TDRHB: R ( t )= R + η ( δ R ( ω ) e − iωt + δ R † ( ω ) e iωt ) , (10)and therefore H ( t )= H + η ( δ H ( ω ) e − iωt + δ H † ( ω ) e iωt ) . (11)The matrices δ R ( ω ) and δ H ( ω ) are not necessarily Hermitian. By linearizing the equationof motion (8) with respect to η , one obtains the linear-response equation in the frequencydomain: ω δ R = [ H , δ R ] + [ δ H ( ω ) , R ] + [ F ( ω ) , R ] . (12)In the stationary quasiparticle basis the matrices H and R are diagonal H = E − E , R = , (13)and, because the density matrix is a projector ( R = R ) at all times, only the two-quasiparticle matrix elements of the time-dependent matrix δ R do not vanish in this basis δ R = R R := XY . (14)This relation defines the QPRA amplitudes X µν and Y µν . In the quasiparticle basis Eq.(12) takes the form ( E µ + E ν − ω ) X µν + δH µν = − F µν (15)( E µ + E ν + ω ) Y µν + δH µν = − F µν . (16)5ince δ H ( ω ) depends on δ R ( ω ), that is, on the amplitudes X µν and Y µν , this is actuallya set of non-linear equations. The expansion of δH µν and δH µν in terms of X µν and Y µν up to linear order leads to the conventional QRPA equations. These equations containsecond derivatives of the density functional E [ R ] with respect to R as matrix elements. Fordeformed nuclei in particular, the number of two-quasiparticle configurations can becomevery large and the evaluation of matrix elements requires a considerable, and in many casesprohibitive, numerical effort. In many cases this has prevented systematic applications ofthe conventional QRPA method to studies of the multipole response of medium-heavy andheavy deformed nuclei.In the finite amplitude method for the QRPA [14, 15], the amplitude X µν and Y µν areformally expressed X µν = − F µν + δH µν E µ + E ν − ω , (17) Y µν = − F µν + δH µν E µ + E ν + ω , (18)and δ H ( ω ) is calculated by numerical differentiation δ H ( ω ) = lim η → η ( H ( R + ηδ R ( ω )) − H ( R )) , (19)using a stationary RHB code for the evaluation of H ( R ). We start from Eq. (14) with δ R ( ω ) in the stationary quasiparticle basis. To use it in the stationary code it has to betransformed back to the original single-particle basis δ R ( ω ) = δρ δκ − δ ¯ κ ∗ − δρ ∗ = W XY W † , (20)and one finds δρ = U XV ⊺ + V ∗ Y U † , (21) δκ = U XU ⊺ + V ∗ Y V † , (22) δ ¯ κ ∗ = − U ∗ Y U † − V XV ⊺ . (23)In this basis we derive the matrix elements of δ H ( ω ) in Eq. (19) δh = lim η → η ( h ( ρ + δρ ) − h ( ρ )) , (24) δ ∆ = lim η → η (∆( κ + δκ ) − ∆( κ )) , (25) δ ¯∆ = lim η → η (∆( κ + δ ¯ κ ) − ∆( κ )) , (26)6nd δH ( ω ) and δH ( ω ) are obtained by transforming back to the quasiparticle basis δ ¯ H ( ω ) = U † V † V ⊺ U ⊺ δh δ ∆ − δ ¯∆ ∗ − δh ⊺ U V ∗ V U ∗ . (27)The explicit expressions for δH and δH read δH ( ω ) = U † δhV ∗ − V † δh ⊺ U ∗ + U † δ ∆ U ∗ − V † δ ¯∆ ∗ V ∗ (28) δH ( ω ) = V ⊺ δhU − U ⊺ δh ⊺ U + V ⊺ δ ∆ V − U ⊺ δ ¯∆ ∗ U . (29)Eqs. (17) and (18) are solved iteratively using the Broyden method [15], and the transitiondensity for each particular frequency ω reads δρ tr ( r ) = − π Im δρ ( r ) . (30)The transition strength is calculated from S ( f, ω ) = − π Im Tr[ f ( U XV ⊺ + V ∗ Y U † )] , (31)and in the present study we only consider isoscalar monopole transitions induced by thesingle-particle operator f = A X i =1 r i . (32) III. NUMERICAL IMPLEMENTATION AND TEST CALCULATIONS
The FAM for the relativistic QRPA is implemented using the stationary RHB code inwhich the single-nucleon Hartree-Bogoliubov equation (1) is solved by expanding the Diracspinors in terms of eigenfunctions of an axially symmetric harmonic oscillator potential (cf.Appendix A). The expressions for the matrix elements of the monopole operator in this basisare given in Appendix B.In the present illustrative study we employ the relativistic functional DD-PC1 [23]. Start-ing from microscopic nucleon self-energies in nuclear matter, and empirical global propertiesof the nuclear matter equation of state, the coupling parameters of DD-PC1 were fine-tunedto the experimental masses of a set of 64 deformed nuclei in the mass regions A ≈ − A ≈ − h k | V S | k ′ i = − Gp ( k ) p ( k ′ )will be here used in the pp channel. By assuming a simple Gaussian ansatz p ( k ) = e − a k ,the two parameters G and a were adjusted to reproduce the density dependence of the gapat the Fermi surface in nuclear matter, calculated with the pairing part of the Gogny inter-action. When transformed from momentum to coordinate space, the interaction takes theform: V ( r , r , r ′ , r ′ ) = Gδ ( R − R ′ ) P ( r ) P ( r ′ ) 12 (1 − P σ ) , (33)where R = ( r + r ) and r = r − r denote the center-of-mass and the relative coordi-nates, respectively, and P ( r ) is the Fourier transform of p ( k ): P ( r ) = 1 / (4 πa ) / e − r / a .The actual implementation of the FAM does not, of course, depend on the choice of therelativistic density functional or the pairing functional.To avoid the occurrence of singularities in the right-hand side of Eqs. (17) and (18), thefrequency ω is replaced by ω + iγ with a small parameter γ , related to the Lorentzian smearingΓ = 2 γ in RQRPA calculations. Eqs. (17) and (18) are solved iteratively. The solution isreached when the maximal difference between collective amplitudes corresponding to twosuccessive iterations decreases below a chosen threshold ( ǫ = 10 − ). The stability and rapidconvergence of the FAM iteration procedure is ensured by adopting the modified Broyden’sprocedure [25, 26], which is also implemented in the calculation of the RHB equilibriumsolution. Compared to ground state calculations, the use of Broyden’s method in the FAMfor QRPA requires an increase of the number of vectors retained in Broyden’s history ( M =20 for the FAM, compared to M = 7 for the RHB). With this modification FAM solutionshave been achieved with less than 40 iterations for all examples considered in the presentillustrative calculations. The FAM for QRPA necessitates the inclusion of time-odd terms(currents) in the calculation of induced fields (cf. Appendix C). The FAM equations for thecase of time-reversal, reflection and axial symmetries are detailed in Appendix D.To verify the numerical implementation and accuracy of our FAM model, a simple testcalculation has been performed for the light spherical nucleus O. In this case we coulddirectly compare the FAM results to those obtained using the standard computer code forthe RQRPA matrix [20]. This comparison presents an excellent test of both codes becausethe FAM formalism employs only numerical derivatives of the single-particle Hamiltonian8 S t r eng t h ( f m M e V - ) RQRPARQFAM 0 5 10 15 20 25 30 35 40Energy (MeV)050100150200250300 O (a) (b) FIG. 1: (Color online) Strength functions of the isoscalar monopole operator for O. The solidcurves denote the RQRPA response, FAM results are indicated by (red) symbols. The two panelscorrespond to a calculation without dynamical pairing (panel (a)), and to a fully self-consistentcalculation with pairing included in the RQRPA residual interaction and FAM induced fields (panel(b)). The single-nucleon wave functions are expanded in a basis of 10 oscillator shells, and theresponse is smeared with a Lorentzian of Γ = 2 γ = 0 . and the pairing field, whereas the QRPA codes uses explicit expressions for the matrixelements of the residual interaction. In Fig. 1 we display the isoscalar strength functions ofthe monopole operator P Ai =1 r i for O. The panel (a) corresponds to a calculation withoutdynamical pairing, that is, pairing is only included in the calculation of the RHB groundstate but not in the residual interaction (QRPA) or induced fields (FAM). The strengthfunctions in the panel (b) are calculated fully self-consistently with dynamical pairing. Inboth panels the solid curves denote the RQRPA response, whereas symbols correspond tothe FAM results. Firstly we note that in both cases the RQRPA and FAM results coincide9xactly at all excitation energies. In the calculation without dynamical pairing, that is,by including pairing correlations only in the RHB ground state, one notices the occurrenceof a strong spurious response below 10 MeV. This Nambu-Goldstone mode is driven toapproximately zero excitation energy (in this particular calculation it is located below 0.2MeV) when pairing correlations are consistently included in the QRPA residual interactionand FAM induced fields.
10 20 30 40 50 60 70 80Energy (MeV)10 -6 -5 -4 -3 -2 -1 ∆ S ( η ) / S ( η ) η =10 -7 η =10 -6 η =10 -5 η =10 -4 η =10 -3 η =10 -2 O FIG. 2: (Color online) The relative accuracy of the strength function for the isoscalar monopoleoperator in O (see Eq. (34)). The curves are plotted for several values of the parameter η andspan a broad interval of excitation energies. Fig. 2 shows the stability of the current implementation of the FAM method for a broadrange of values of the parameter η that is used to calculate the numerical derivatives inEqs. (24) – (26). The relative accuracy of the strength function is defined as∆ S ( ω, η ) S ( ω, η ) = 1 S ( ω, η ) | S ( ω, η ) − S ( ω, η ) | . (34)10n practice the accuracy can only be improved by reducing η down to 10 − . A furtherdecrease of this parameter introduces numerical noise which deteriorates the accuracy of theFAM method, and thus η = 10 − has been used throughout this study. IV. ILLUSTRATIVE CALCULATIONS: SAMARIUM ISOTOPES
Collective nucleonic oscillations along different axes in deformed nuclei and mixing ofdifferent modes lead to a broadening and splitting of giant resonance structures [27]. Thegiant dipole resonance (GDR), for instance, displays a two-component structure in deformednuclei and the origin of this splitting are the different frequencies of oscillations along themajor and minor axes. In axially deformed nuclei the isoscalar giant quadrupole (ISGQR)resonance displays three components with K π = 0 + , + , + [28, 29], where K denotes the pro-jection of the total angular momentum I = 2 + on the intrinsic symmetry axis. The isoscalargiant monopole resonance (ISGMR) in deformed nuclei mixes with the K π = 0 + componentof the ISGQR and a two-peak structure of the monopole resonance is observed [30, 31]. In arecent study of the roles of deformation and neutron excess on the giant monopole resonancein neutron-rich deformed Zr isotopes [32], based on the deformed Skyrme – HFB + QRPAmodel, the evolution of the two-peak structure of the ISGMR has been investigated. Thetheoretical analysis has shown that the lower peak is associated with the mixing betweenthe ISGMR and the K π = 0 + component of the ISGQR, and the transition strength of thelower peak increases with neutron excess. Here we apply the FAM method for the rela-tivistic QRPA to a calculation of the isoscalar K π = 0 + strength functions in the chain ofeven-even Sm isotopes, starting from the neutron-deficient Sm isotope and extending tothe neutron-rich
Sm isotope. The calculations have been performed in the harmonic os-cillator basis with N max = 18 oscillator shells for the upper component and N max = 19 shellsfor the lower component of the Dirac spinors [33]. It has been demonstrated in Ref. [15]that by using N max = 18 oscillator shell basis one obtains convergent results even for thesuperdeformed states.Fig. 3 displays the energy curves of Sm isotopes calculated with the constraint on theaxial quadrupole moment, as functions of the axial deformation parameter β . Energies arenormalized with respect to the binding energy of the absolute minimum for each isotope. Forthe isotopes with a prolate equilibrium deformation ( − Sm and − Sm), an additional11 B i nd i ng ene r g y ( M e V ) β β β Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm(a) (b) (c)(d) (e) (f)(g) (h) (j)(k) (l) (m)(n) (o) (p)
FIG. 3: (Color online) Self-consistent RHB binding energy curves of the even-even − Smisotopes as functions of the axial deformation parameter β . Energies are normalized with respectto the binding energy of the absolute minimum for each isotope. minimum is predicted on the oblate side and the two minima are separated by a potentialbarrier. In neutron-deficient isotopes both the oblate minimum and the potential barrierare considerably lower compared to the neutron-rich nuclei. Both Sm and
Sm, that is,nuclides at the borders of the region of weakly deformed and/or spherical systems around theneutron shell-closure at N = 82, exhibit soft potentials with wide minima on the prolate side. − Sm display two weakly deformed and almost degenerated minima, and the isotopes , Sm are spherical.For each isotope in the chain − Sm the calculated K π = 0 + response is shown inFig. 4. The principal result is the splitting of the K π = 0 + strength into two peaks for thedeformed isotopes. The arrows indicate the positions of the mean energies m /m , that is,the ratio of the energy-weighted sum (EWS) and the non-energy-weighted sum, calculated12
46 246 246246 246 24682468 S t r eng t h ( f m M e V - ) Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm Sm(a) (b) (c)(d) (e) (f)(g) (h) (i)(j) (k) (l)(m) (n) (o)
FIG. 4: (Color online) Evolution of the K π = 0 + strength functions in − Sm. The arrowsindicate the positions of the mean energies m /m calculated in the energy intervals 10 < E < . . < E <
20 MeV. in the energy intervals 10 < E < . . < E < Sm, whereasthe LE peak appears in the energy region where the giant quadrupole resonance in
Sm islocated ( E ISGQR = 14 MeV). With increasing deformation (cf. Fig. 3) the HE peak is shiftedto higher energy because of the coupling with the K π = 0 + component of the ISGQR, andthe LE peak is simultaneously lowered in energy. It should be noted that the K π = 0 + components of other resonances also contribute to the LE and HE peaks, but to a muchlesser extent.In the panel (b) of Fig. 5 we display the mean energies of the HE (squares) and LE (circles)peaks as functions of the equilibrium deformation parameter β . The calculated energies are13 .0 0.1 0.2 0.3 0.4 0.5 β E . A / ( M e V ) Sm, LE
Sm, HE
Sm, LE
Sm, HE
Sm, J π =2 +144 Sm, J π =0 + S , f r a c t i on o f t he E W S R (a)(b) FIG. 5: (Color online) Panel (a): fraction of the EWSR for the HE and LE components of the K π = 0 + strength in the deformed nuclei − Sm and − Sm, calculated in the energyintervals 14 . < E <
20 MeV and 10 < E < . m /m of the HEand LE peaks, denoted by squares and circles, respectively, as functions of the equilibrium valueof β . A / to account for the empirical mass dependence of the ISGMRexcitation energy E ∼ A − / . With increasing equilibrium deformation the splitting betweenthe LE and HE components becomes larger, although the trend is not quite the same for theisotopes with A <
144 and
A > < E < . . < E <
20 MeV for the high-energy (HE) region. Generally thefraction of the EWSR in the LE mode increases with deformation, but again the trend isslightly different for for the isotopes with
A <
144 and
A > K π = 0 + strength in deformed systems can be studiedin more detail by performing a deformation-constrained calculation for a single isotope. InFig. 6 we show the K π = 0 + strength distributions in Sm isotope for eight different valuesof the axial quadrupole constraint, from β = − . β = 0 .
4. The blue dashed and reddot-dashed curves correspond to the monopole and quadrupole strength distributions for theequilibrium spherical configuration of
Sm, respectively. Both for the prolate and oblateconstrained configurations the splitting between the LE and HE components of the K π = 0 + strength increases with deformation. An interesting result is that the HE component of the K π = 0 + strength distribution is more pronounced for prolate configurations, whereas foroblate configurations the LE component becomes dominant.Fig. 7 compares the mean energies, that is, the ratio of the energy-weighted sum (EWS)and the non-energy-weighted sum m /m of the HE and LE components of the K π = 0 + strength distribution in Sm, as functions of the constrained quadrupole deformation β .For the prolate configurations the moments of the strength distribution are calculated inthe energy intervals 10 − . . −
20 MeV (HE region). Thecorresponding intervals for oblate configurations are 10 − . . − . K π = 0 + distribution is consistent with the result shown in Fig. 1 of Ref. [37]. We note15 S t r eng t h ( f m M e V - )
12 14 16 18 20 22 2424681012 14 16 18 20 22 24246810 12 14 16 18 20 22 2424681012 14 16 18 20 22 24
E (MeV)
12 14 16 18 20 22 24
E (MeV) β =0.1 β =0.1 β =0.2 β =0.3 β =0.4 β =-0.1 β =-0.2 β =-0.3 β =-0.4 Sm (a) (e)(b) (f)(g)(c)(d) (h) FIG. 6: (Color online) The K π = 0 + strength distributions in Sm for eight different values ofthe axial quadrupole constraint. The blue dashed and red dot-dashed curves denote the J = 0 + and J = 2 + strengths for the Sm equilibrium spherical configuration, respectively. that the monopole ( I = 0 + ) and quadrupole ( I = 2 + ) strength distributions, calculated forthe spherical equilibrium configuration of Sm, are somewhat fragmented and this leads tothe broadening of the K π = 0 + strength distribution for deformed configurations since eachmonopole state couples to each quadrupole state.Individual modes of collective excitations can be studied qualitatively by analyzing thecorresponding transition densities. For K π = 0 + the intrinsic transition densities are axiallysymmetric: δρ tr ( r ) = δρ tr ( r ⊥ , z ) . (35)By projecting the two dimensional intrinsic transition densities δρ tr ( r ⊥ , z ) onto good angularmomentum, one obtains the transition densities in the laboratory frame of reference. For a16 β E ne r g y ( M e V ) HE regionLE region Sm FIG. 7: (Color online) The mean energies m /m of the HE and LE components of the K π = 0 + strength distribution in Sm, as functions of the constrained quadrupole deformation β . particular value of the angular momentum J ≥ K , the projected transition density reads δρ Jtr ( r ) = δρ Jtr ( r ) Y JK (Ω) , (36)with the radial part of the projected transition density δρ Jtr ( r ) = Z d Ω δρ tr ( r ⊥ , z ) Y JK (Ω) . (37)Although the last equation is not exact, it yields accurate results for large deformations.As an example we chose two axially-constrained deformed configurations of Sm: theoblate configuration at deformation β = − .
3, and the prolate configuration at β = 0 . J = 0 and J = 2 angular-momentum-projected transition densities, and the intrinsictransition densities for the LE and HE peaks of the ISGMR strength distributions are shownin Figs. 8 (prolate deformed configuration, β = 0 .
3) and 9 (oblate deformed configuration,17 = − . J = 0 and J = 2 components of the transition density is constructive at thepoles and destructive in the equatorial plane of the density ellipsoid of Sm. Figure 10compares the radial parts of the angular momentum projected transition densities δρ J =0 tr ( r )and δρ J =2 tr ( r ) that correspond to the LE and HE peaks in the Sm isotope: the prolateconfiguration at β = 0 . β = − . δρ J =0 tr ( r ) component displays the characteristic radial dependenceof the monopole (compression) transition strength with a single node in the surface region.Both the volume and surface contributions are more pronounced for the HE componentat prolate deformation, whereas for the oblate deformed configuration the LE componentdominates. This is consistent with the strength distributions displayed in Fig. 7. In all cases δρ J =2 tr ( r ) has a radial dependence characteristic for quadrupole oscillations. We also noticethat for the LE component the surface contributions of the δρ J =0 tr ( r ) and δρ J =2 tr ( r ) transitiondensities are in phase when the nucleus has prolate deformation, and out of phase when thedeformation is oblate. The opposite is found for the HE energy component. V. SUMMARY AND OUTLOOK
Realistic QRPA calculations for deformed nuclei still present a considerable computa-tional challenge, particularly if one considers heavy nuclei. The dimension of the configu-ration space increases rapidly in heavier open-shell nuclei, and thus it becomes increasinglydifficult to compute and store huge QRPA matrices. Although several relativistic QRPAstudies have been performed for axially deformed nuclei, computationally this task is simplyto complex for systematic large scale calculations. One possible solution is to employ thefinite amplitude method in the solution of the corresponding linear response problem. In thiswork we have implemented a recently proposed efficient method for the iterative solutionof the FAM – QRPA equations in the framework of relativistic energy density functionals.Several numerical tests have been performed to verify the stability of the FAM-RQRPAiterative solution and its consistency with the solution of the matrix RQRPA equations.As an illustrative example, the FAM-RQRPA model has been applied to an analysis of the18 x −10−8−6−4−20246810 z (a) Sm, LE mode −0.30−0.24−0.18−0.12−0.060.000.060.120.180.24 δ ρ J = (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (d) Sm, HE mode −0.30−0.24−0.18−0.12−0.060.000.060.120.180.240.300.360.420.48 δ ρ J = (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (b) −0.30−0.24−0.18−0.12−0.060.000.060.120.180.24 δ ρ J = (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (e) −0.30−0.24−0.18−0.12−0.060.000.060.120.180.240.300.360.420.48 δ ρ J = (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (c) −0.30−0.24−0.18−0.12−0.060.000.060.120.180.24 δ ρ t o t (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (f) −0.30−0.24−0.18−0.12−0.060.000.060.120.180.240.300.360.420.48 δ ρ t o t (cid:0) f m − (cid:1) FIG. 8: (Color online) The J = 0 and J = 2 angular-momentum-projected transition densities,and the intrinsic transition densities for the LE (left column) and HE (right column) peaks ofthe ISGMR strength distribution in Sm. The stationary density corresponds to the prolateconfiguration with the constraint deformation β = 0 . x −10−8−6−4−20246810 z (a) Sm, LE mode −0.30−0.24−0.18−0.12−0.060.000.060.120.180.240.300.360.420.48 δ ρ J = (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (d) Sm, HE mode −0.30−0.24−0.18−0.12−0.060.000.060.120.180.240.300.36 δ ρ J = (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (b) −0.30−0.24−0.18−0.12−0.060.000.060.120.180.240.300.360.420.48 δ ρ J = (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (e) −0.30−0.24−0.18−0.12−0.060.000.060.120.180.240.300.36 δ ρ J = (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (c) −0.30−0.24−0.18−0.12−0.060.000.060.120.180.240.300.360.420.48 δ ρ t o t (cid:0) f m − (cid:1) x −10−8−6−4−20246810 z (f) −0.30−0.24−0.18−0.12−0.060.000.060.120.180.240.300.36 δ ρ t o t (cid:0) f m − (cid:1) FIG. 9: (Color online) Same as in the caption to Fig. 9 but for the oblate configuration with theconstraint deformation β = − . Sm. δ ρ t r J ( f m - ) δ ρ t r J ( f m - ) β =0.3, LE β =0.3, HE β =-0.3, LE β =-0.3, HE Sm (a) (b)(c) (d) FIG. 10: (Color online) Radial parts of the angular-momentum-projected transition densities thatcorrespond to the LE and HE ISGMR peaks in the
Sm isotope: the prolate configuration at β = 0 . β = − . splitting of the giant monopole resonance in deformed nuclei. In particular, we have inves-tigated the evolution of the K π = 0 + strength in the chain of samarium isotopes, from theproton-rich Sm, to systems with considerable neutron excess close to
Sm. To studythe splitting of the monopole strength in more detail, we have performed a deformation-constrained FAM-RQRPA calculation for the nucleus
Sm. A significant mixing of themonopole and quadrupole modes has been found for both the low-energy and high-energycomponents of the K π = 0 + strength, consistent with the standard interpretation of thesplitting of monopole strength in deformed nuclei.The advantage of developing and employing the FAM-RQRPA formalism would, of course,be limited unless it is extended to higher multipoles. This development is already in progress.Another important extension of the FAM-RQRPA model is the one to the charge-exchange21hannel. The FAM for the charge-exchange RQRPA will enable the description and modelingof a variety of astrophysically relevant weak-interaction processes, in particular beta-decay,electron capture, and neutrino reactions in deformed nuclei. Since many isotopes that arecrucial for the process of nucleosynthesis display considerable deformation, a reliable mod-eling of elemental abundances necessitates a microscopic and self-consistent description ofunderlying transitions and weak-interaction processes, and this can be attained using thecharge-exchange extension of the framework introduced in this work. Acknowledgments
Appendix A: The single-nucleon basis
The Dirac single-nucleon spinors are expanded in the basis of eigenfunctions of an axiallysymmetric harmonic oscillator in cylindrical coordinates:Φ α ( r ⊥ , z, φ, s ) = φ n z ( z ) φ Λ n ⊥ ( r ⊥ ) 1 √ π e i Λ φ χ m s , α ≡ { Λ n ⊥ n z m s } , (A1)where φ n z ( z ) = N n z √ b z H n z ( ξ ) e − ξ / , ξ = z/b z , N n z = 1 p √ π n z n z ! , (A2)with the Hermite polynomial H n z ( z ), and for the r ⊥ coordinate φ Λ n ⊥ ( r ⊥ ) = N Λ n ⊥ b ⊥ √ η Λ / L Λ n ⊥ ( η ) e − η/ , η = r ⊥ b ⊥ , N Λ n ⊥ = s n ⊥ !( n ⊥ + Λ)! , (A3)with the Laguerre polynomial L Λ n ⊥ ( η ). The time-reversed state readsΦ ¯ α ( r ⊥ , z, φ, s ) = ˆ T Φ α ( r ⊥ , z, φ, s ) = φ n z ( z ) φ Λ n ⊥ ( r ⊥ ) 1 √ π e − i Λ φ ( − / − m s χ − m s . (A4)22 ppendix B: Matrix elements of the monopole operator The matrix elements of the monopole operator can be calculated analytically in theharmonic oscillator basis. After separating the spin and angular parts of the matrix element,the following expression is obtained: f αα ′ = δ m s ,m ′ s δ ΛΛ ′ Z ∞−∞ dz Z ∞ dr ⊥ r ⊥ φ Λ n ⊥ ( r ⊥ ) φ n z ( z )( r ⊥ + z ) φ Λ ′ n ′⊥ ( r ⊥ ) φ n ′ z ( z ) . (B1)Using the orthogonality of the eigenfunctions in the z and r ⊥ coordinates, the matrix elementcan be written as f αα ′ = δ m s ,m ′ s δ ΛΛ ′ × (cid:20) δ n z n ′ z Z ∞ dr ⊥ r ⊥ φ Λ n ⊥ ( r ⊥ ) φ Λ ′ n ′⊥ ( r ⊥ ) + δ n ⊥ n ′⊥ Z ∞−∞ dzz φ n z ( z ) φ n ′ z ( z ) (cid:21) . (B2)For the integral over z direction one uses the recursive relation ξ H n ′ z ( ξ ) = 14 H n ′ z +2 ( ξ ) + 12 (2 n ′ z + 1) H n ′ z ( ξ ) + n ′ z ( n ′ z − H n ′ z − ( ξ ) , (B3)which can be expressed in terms of harmonic oscillator eigenfunctions ξ φ n ′ z ( ξ ) = 12 p ( n ′ z + 2)( n ′ z + 1) φ n ′ z +2 ( ξ ) + 12 (2 n ′ z + 1) φ n ′ z ( ξ ) + 12 p n ′ z ( n ′ z − φ n ′ z − ( ξ ) . (B4)This yields I z = Z ∞−∞ φ n z ( z ) z φ n ′ z ( z ) dz = p n ′ z ( n ′ z − b z , n z = n ′ z − (cid:0) n z + (cid:1) b z , n z = n ′ z p n z ( n z − b z , n z = n ′ z + 20 , otherwise . (B5)To calculate the integral in r ⊥ the recursive relation ηL Λ n ′ ( η ) = (2 n ′ + Λ + 1) L Λ n ′ ( η ) − ( n ′ + 1) L Λ n ′ +1 ( η ) − ( n ′ + Λ) L Λ n ′ − ( η ) , (B6)can expressed in the form ηφ Λ n ′⊥ ( r ⊥ ) = (2 n ′⊥ + Λ + 1) φ Λ n ′⊥ ( r ⊥ ) − q n ′⊥ ( n ′⊥ + Λ) φ Λ n ′⊥ − ( r ⊥ ) (B7) − q ( n ′⊥ + 1)( n ′⊥ + Λ + 1) φ Λ n ′⊥ +1 ( r ⊥ ) . I ⊥ = Z ∞ φ Λ n ⊥ ( η ) ηφ Λ n ′⊥ ( η ) dη = (2 n ⊥ + Λ + 1) b ⊥ , n ′⊥ = n ⊥ − p n ′⊥ ( n ′⊥ + Λ) b ⊥ , n ′⊥ = n ⊥ + 1 − p n ⊥ ( n ⊥ + Λ) b ⊥ , n ′⊥ = n ⊥ −
10 otherwise (B8)The monopole operator does not mix states from different K π blocks, and the matrix ele-ments are real and symmetric. Appendix C: Time-odd terms
The time-odd current reads j ( r ) = X α ˜ α h ρ α ˜ α Φ † α σ Φ ˜ α + ρ ˜ αα Φ † ˜ α σ Φ α i , (C1)where α ( ˜ α ) denotes the harmonic oscillator quantum numbers for the large (small) com-ponent of the single-nucleon Dirac spinor. The σ matrix can be expressed in cylindricalcoordinates σ = e − iφ σ + e ⊥ + e iφ σ − e ⊥ − ie − iφ σ + e φ + ie iφ σ − e φ + σ z e z . (C2)The following expressions can easily be evaluatedΦ † α σ + e − iφ Φ β = 12 π δ m αs , / δ m βs , − / φ n αz ( z ) φ n βz ( z ) φ Λ α n α ⊥ ( r ⊥ ) φ Λ β n β ⊥ ( r ⊥ ) e i (Λ β − Λ α − φ , (C3)Φ † α σ − e iφ Φ β = 12 π δ m αs , − / δ m βs , / φ n αz ( z ) φ n βz ( z ) φ Λ α n α ⊥ ( r ⊥ ) φ Λ β n β ⊥ ( r ⊥ ) e i (Λ β − Λ α +1) φ . (C4)Next, we use the condition for monopole excitations Ω α = Ω β , that is, Λ α + m αs = Λ β + m βs ,Φ † α (cid:0) σ + e − iφ + σ − e iφ (cid:1) Φ β = 12 π δ m αs , − m βs φ n αz ( z ) φ n βz ( z ) φ Λ α n α ⊥ ( r ⊥ ) φ Λ β n β ⊥ ( r ⊥ ) , (C5)Φ † α (cid:0) σ − e iφ − σ + e − iφ (cid:1) Φ β = 12 π ( − / − m βs δ m αs , − m βs φ n αz ( z ) φ n βz ( z ) φ Λ α n α ⊥ ( r ⊥ ) φ Λ β n β ⊥ ( r ⊥ ) . (C6)and, finally, calculate the contribution from the z componentΦ † α σ z Φ β = 12 π ( − / − m βs δ m αs ,m βs φ n αz ( z ) φ n βz ( z ) φ Λ α n α ⊥ ( r ⊥ ) φ Λ β n β ⊥ ( r ⊥ ) . (C7)The following relations are validΦ † α (cid:0) σ + e − iφ + σ − e iφ (cid:1) Φ β = Φ † β (cid:0) σ + e − iφ + σ − e iφ (cid:1) Φ α , (C8)Φ † α (cid:0) σ − e iφ − σ + e − iφ (cid:1) Φ β = − Φ † β (cid:0) σ − e iφ − σ + e − iφ (cid:1) Φ α , (C9)Φ † α σ z Φ β = Φ † β σ z Φ α , (C10)24nd also Φ † ¯ α (cid:0) σ + e − iφ + σ − e iφ (cid:1) Φ ¯ β = − Φ † α (cid:0) σ + e − iφ + σ − e iφ (cid:1) Φ β , (C11)Φ † ¯ α (cid:0) σ − e iφ − σ + e − iφ (cid:1) Φ ¯ β = Φ † α (cid:0) σ − e iφ − σ + e − iφ (cid:1) Φ β , (C12)Φ † ¯ α σ z Φ ¯ β = − Φ † α σ z Φ β . (C13)The corresponding elements of the Hamiltonian matrix read h α | σ · V | β i = h α | ( e − iφ σ + + e iφ σ − ) V ⊥ + i ( e iφ σ − − e − iφ σ + ) V φ + σ z V z | β i . (C14) Appendix D: FAM equations for time-reversal symmetry
We consider systems with time-reversal, reflection and axial symmetries. The single-quasiparticle states can be ordered so that we first list states with Ω >
0, and then stateswith Ω <
0. The HFB matrices U and V read U = u u ∗ , V = − v ∗ v . (D1)This generates a density matrix and pairing tensor with block-diagonal structure ρ = V ∗ V T = v ∗ − v ∗ v ∗ − v T = vv † v ∗ v T = ρ ρ , (D2) κ = V ∗ U T = v ∗ − v ∗ u u ∗ T = vu † − v ∗ u T = κ κ . (D3)The FAM amplitudes X and Y are antisymmetric matrices X = x − x T , Y = y − y T , (D4)where x and y are symmetric complex matrices. The explicit expressions for the densitymatrix and pairing tensor read ρ = ( v − ηux ) ( v − ηuy ∗ ) † , ρ = (cid:0) v − ηux † (cid:1) ∗ (cid:0) v − ηuy T (cid:1) T ,κ = (cid:0) v − ηux † (cid:1) ∗ (cid:0) u + ηvy T (cid:1) T , κ = − ( v − ηux ) ( u + ηvy ∗ ) † , ¯ κ = (cid:0) v − ηuy T (cid:1) ∗ (cid:0) u + ηvx † (cid:1) T , ¯ κ = − ( v − ηuy ∗ ) ( u + ηvx ) † . (D5)25t should be noted that since the x and y matrices are complex, the relations ρ = ρ ∗ and κ † i = κ i are no longer fulfilled. The matrices H ( ω ) and H ( ω ) read δH ( ω ) = δh − [ δh ] T , δH ( ω ) = δh − [ δh ] T (D6)with δh ( ω ) = − u † δh v + u † δ ∆ u + v † δ ¯∆ v − v † δh T u, (D7) δh ( ω ) = v T δh u ∗ − u T δ ¯∆ ∗ u ∗ − v T δ ∆ v ∗ + u T δh T v ∗ . (D8)The matrices F and F of the external operator are decomposed in an analogous way.Time-reversal symmetry reduces by half the dimension of the equations of motion( E µ + E ν − ω ) x µν + δh µν + f µν = 0 , (D9)( E µ + E ν + ω ) y µν + δh µν + f µν = 0 . (D10) [1] N. Paar, D. Vretenar, E. Khan, and G. Colo, Rep. Prog. Phys. 70, 691 (2007).[2] D. Savran, T. Aumann, and A. Zilges, Prog. Part. Nucl. Phys. 70, 210 (2013).[3] P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer-Verlag, Heidelberg, 1980).[4] M. Bender, P.-H. Heenen, P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).[5] Extended Density Functionals in Nuclear Structure Physics, Lecture Notes in Physics 641,edited by G. A. Lalazissis, P. Ring, and D. Vretenar (Springer, Heidelberg, Germany, 2004).[6] S. P´eru and H. Goutte, Phys. Rev. C 77, 044313 (2008).[7] K. Yoshida and N. Van Giai, Phys. Rev. C 78, 064316 (2008).[8] D. Pe˜na Arteaga, E. Khan, and P. Ring, Phys. Rev. C 79, 034311 (2009).[9] J. Terasaki and J. Engel, Phys. Rev. C 82, 034326 (2010).[10] C. Losa, A. Pastore, T. Dossing, E. Vigezzi, and R. A. Broglia, Phys. Rev. C 81, 064307(2010).[11] T. Nakatsukasa, T. Inakura, and K. Yabana, Phys. Rev. C 76, 024318 (2007).[12] T. Inakura, T. Nakatsukasa, and K. Yabana, Phys. Rev. C 80, 044301 (2009).[13] T. Inakura, T. Nakatsukasa, and K. Yabana, Phys. Rev. C 84, 021302 (2011).
14] P. Avogadro, and T. Nakatsukasa, Phys. Rev. C 84, 014314 (2011).[15] M. Stoitsov, M. Kortelainen, T. Nakatsukasa, C. Losa, and W. Nazarewicz, Phys. Rev. C 84,041305 (R) (2011).[16] N. Hinohara, M. Kortelainen, and W. Nazarewicz, Phys. Rev. C 87, 064309 (2013).[17] H. Liang, T. Nakatsukasa, Z. Niu, J. Meng, Phys. Rev. C 87, 054310 (2013).[18] D. Vretenar, A. V. Afansjev, G. A. Lalazissis, P. Ring, Phys. Rep. 409, 101 (2005).[19] J. Meng, H. Toki, S. G. Zhou, S. Q. Zhang, W. H. Long, and L. S. Geng, Prog. Part. Nucl.Phys. 57, 470 (2006).[20] N. Paar, P. Ring, T. Nikˇsi´c and D. Vretenar, Phys. Rev. C 67, 034312 (2003).[21] N. Paar, T. Nikˇsi´c, D. Vretenar and P. Ring, Phys. Rev. C 69, 054303 (2004).[22] J. G. Valatin, Phys. Rev. 122, 1012 (1961)[23] T. Nikˇsi´c, D. Vretenar, P. Ring, Phys. Rev. C 78, 034318 (2008).[24] Y. Tian, Z.Y. Ma, P. Ring, Phys. Lett. B 676, 44 (2009).[25] D.D. Johnson, Phys. Rev. B 38, 12807 (1988).[26] A. Baran, A. Bulgac, M. McNeil Forbes, G. Hagen, W. Nazarewicz, N. Schunck, M.V. Stoitsov,Phys. Rev. C 78, 014318 (2008).[27] M. N. Harakeh and A. van der Wounde, Giant Resonances: Fundamental High-Energy Modesof Nuclear Excitation (Oxford University Press, Oxford, 2001).[28] T. Kishimoto, J. M. Moss, D. H. Youngblood, J. D. Bronson, C. M. Rozsa, D. R. Brown, andA. D. Bacher, Phys. Rev. Lett. 35, 552 (1975).[29] H. Miura, and Y. Torizuka, Phys. Rev. C 16, 1688 (1977).[30] U. Garg, P. Bogucki, J. D. Bronson, Y.-W. Lui, C. M. Rozsa, and D. H. Youngblood, Phys.Rev. Lett. 45, 1670 (1980).[31] H. P. Morsch, M. Rogge, P. Turek, C. Mayer-Boricke, and P. Decowski, Phys. Rev. C 25, 2939(1982).[32] K. Yoshida, Phys. Rev. C 82, 034324 (2010).[33] Y. K. Gambhir, P. Ring, and A. Thimet, Ann. Phys. (N.Y.) 198, 132 (1990).[34] C.E. Price and G.E. Walker, Phys. Lett. B 155, 17 (1985).[35] J.R. Shepard, E. Rost and J.A. McNeil, Phys. Rev. C 40, 2320 (1989).[36] J.A. McNeil, R.J. Furnstahl, E. Rost and J.R. Shepard, Phys. Rev. C 40, 399 (1989).[37] S. Nishizaki, K. Ando, Prog. Theor. Phys. 73, 889 (1985).14] P. Avogadro, and T. Nakatsukasa, Phys. Rev. C 84, 014314 (2011).[15] M. Stoitsov, M. Kortelainen, T. Nakatsukasa, C. Losa, and W. Nazarewicz, Phys. Rev. C 84,041305 (R) (2011).[16] N. Hinohara, M. Kortelainen, and W. Nazarewicz, Phys. Rev. C 87, 064309 (2013).[17] H. Liang, T. Nakatsukasa, Z. Niu, J. Meng, Phys. Rev. C 87, 054310 (2013).[18] D. Vretenar, A. V. Afansjev, G. A. Lalazissis, P. Ring, Phys. Rep. 409, 101 (2005).[19] J. Meng, H. Toki, S. G. Zhou, S. Q. Zhang, W. H. Long, and L. S. Geng, Prog. Part. Nucl.Phys. 57, 470 (2006).[20] N. Paar, P. Ring, T. Nikˇsi´c and D. Vretenar, Phys. Rev. C 67, 034312 (2003).[21] N. Paar, T. Nikˇsi´c, D. Vretenar and P. Ring, Phys. Rev. C 69, 054303 (2004).[22] J. G. Valatin, Phys. Rev. 122, 1012 (1961)[23] T. Nikˇsi´c, D. Vretenar, P. Ring, Phys. Rev. C 78, 034318 (2008).[24] Y. Tian, Z.Y. Ma, P. Ring, Phys. Lett. B 676, 44 (2009).[25] D.D. Johnson, Phys. Rev. B 38, 12807 (1988).[26] A. Baran, A. Bulgac, M. McNeil Forbes, G. Hagen, W. Nazarewicz, N. Schunck, M.V. Stoitsov,Phys. Rev. C 78, 014318 (2008).[27] M. N. Harakeh and A. van der Wounde, Giant Resonances: Fundamental High-Energy Modesof Nuclear Excitation (Oxford University Press, Oxford, 2001).[28] T. Kishimoto, J. M. Moss, D. H. Youngblood, J. D. Bronson, C. M. Rozsa, D. R. Brown, andA. D. Bacher, Phys. Rev. Lett. 35, 552 (1975).[29] H. Miura, and Y. Torizuka, Phys. Rev. C 16, 1688 (1977).[30] U. Garg, P. Bogucki, J. D. Bronson, Y.-W. Lui, C. M. Rozsa, and D. H. Youngblood, Phys.Rev. Lett. 45, 1670 (1980).[31] H. P. Morsch, M. Rogge, P. Turek, C. Mayer-Boricke, and P. Decowski, Phys. Rev. C 25, 2939(1982).[32] K. Yoshida, Phys. Rev. C 82, 034324 (2010).[33] Y. K. Gambhir, P. Ring, and A. Thimet, Ann. Phys. (N.Y.) 198, 132 (1990).[34] C.E. Price and G.E. Walker, Phys. Lett. B 155, 17 (1985).[35] J.R. Shepard, E. Rost and J.A. McNeil, Phys. Rev. C 40, 2320 (1989).[36] J.A. McNeil, R.J. Furnstahl, E. Rost and J.R. Shepard, Phys. Rev. C 40, 399 (1989).[37] S. Nishizaki, K. Ando, Prog. Theor. Phys. 73, 889 (1985).