{'fraction_non_alphanumeric': 0.035774911554265224, 'fraction_numerical': 0.01139085875652436, 'mean_word_length': 4.4893320560829055, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Continuum Random Phase Approximation for Relativistic Point Coupling Models

Relativistic Continuum Random Phase Approximation (CRPA) is used to investigate collective excitation phenomena in several spherical nuclei along the periodic table. We start from relativistic mean field calculations based on a covariant density functional with density dependent zero range forces. From the same functional an effective interaction is obtained as the second derivative with respect to the density. This interaction is used in relativistic continuum-RPA calculations for the investigation of isoscalar monopole, isovector dipole and isoscalar quadrupole resonances of spherical nuclei. In particular we study the low-lying E1 strength in the vicinity of the neutron evaporation threshold. The properties of the resonances, such as centroid energies and strengths distributions are compared with results of discrete RPA calculations for the same model as well as with experimental data.


Introduction
Density functional theory (DFT) provides a very successful description of nuclei all over the periodic table. Based on relatively simple functionals, which are adjusted in a phenomenological way to the properties of infinite nuclear matter and a few finite nuclei, this theory allows a highly accurate reproduction of many nuclear structure data, such as binding energies, radii, deformation parameters of finite nuclei and their dependence on mass number and isospin. In addition to these static properties, one can use the nuclear response to external multipole fields to investigate the dynamics of such systems. In the framework of time dependent density functional theory, this response can be calculated from the linearized Bethe Salpeter equation using an effective interaction derived from the same functional.
A very successful scheme of this type is covariant density functional theory (CDFT). It is based on Lorentz invariance, connecting in a consistent way the spin and spatial degrees of freedom of the nucleus. Therefore, it needs only a relatively small number of parameters which are adjusted to reproduce a set of bulk properties of spherical closed-shell nuclei. Numerous works have shown that observations involving both ground state and excited state phenomena, can be nicely interpreted in a relativistic framework.
The most popular applications of this type are based on the Walecka model [1], where the nucleus is described as a system of Dirac nucleons interacting with each other via the exchange of virtual mesons with finite mass and the electromagnetic field through an effective Lagrangian. In the mean field approximation this yields to various contributions to the nuclear self energy depending on the quantum numbers of these mesons. Early investigations have shown that this simple ansatz is not able to describe the incompressibility of infinite nuclear matter nor the surface properties of finite nuclei such as nuclear deformations. For that reason, a medium dependence has been introduced by including nonlinear meson self-interaction terms in the Lagrangian [2].
Several very successful phenomenological RMF interactions of this type have been adopted, as for instance the popular set NL3 [3]. Closer to the concept of density functional theory are models with an explicit density dependence for the meson nucleon couplings. This density dependence can be calculated from first principles in a microscopic Dirac-Brueckner scheme [4] or it can be adjusted in a completely phenomenological way to properties of finite nuclei [5,6].
One of the advantages of density functional theory is the fact that with a proper choice of the parameters the success of RMF for nuclear ground states ensures also a good basis upon which one can apply time-dependent density functional theory to study nuclear excitations. In order to investigate the dynamic behavior of the nuclear system, one considers oscillations around the self-consistent static solution. This can be done by solving the time dependent relativistic mean field equations (TDRMF) [7] or, in the limit of small amplitudes, by using the relativistic random phase approximation (RRPA) [8]. The corresponding eigen modes can be determined either by diagonalizing the RRPA equation in an appropriate basis or by solving the linear response equations in a time-dependent external field. This requires a matrix inversion for given frequency ω.
These two methods lead in principle to exactly identical results. There are, however, cases where one of them is clearly preferable. The proper treatment of the coupling to the continuum is such a case, which can be solved in a very elegant way, by the solution of the Bethe Salpeter equation within the response formalism.
We recall that the spectrum of the Dirac equations has a discrete and a continuous part. For the ground state properties of the nucleus, one needs only the single particle wave functions of the occupied orbitals in the Fermi sea. They are either determined by solving the corresponding differential equations in r-space or by expansion in an appropriate basis, given for for instance by a finite number of eigenfunctions of a harmonic oscillator [9] or of a Saxon-Woods potential in a finite box [10]. For the bound states both methods yield the same solutions with high accuracy. However, this is no longer true for the states in the continuum. Here we have, in the first case scattering solutions in r-space for each energy with proper boundary conditions, while in the second case, a finite number of discrete eigenstates which depend strongly on the dimension of the expansion. They provide only a basis and have little to do with physics.
These discrete eigenstates lead to a finite number of ph-configurations for the solution of response equations. with a discrete spectrum. They provide us with the so called spectral representation of the response function in contrast to the continuum representation, where the exact scattering states with the proper boundary conditions are used at each energy.
Self-consistent relativistic RPA (RRPA) calculations have a long history. The early investigations in the eighties [11,13,14,15] were based on the Walecka model with linear meson-nucleons couplings. They were able to describe the low-lying negative-parity excitations in 16 O by the method of matrix diagonalization [11], isoscalar giant resonances in light and medium nuclei [13] by the solution of the linear response equation in the spectral representation, and the longitudinal response for quasi-elastic electron scattering with a proper treatment of the continuum.
The first RRPA calculations based on non-linear models were carried out in the spectral representation including only normal particle-hole (ph) pairs with particles above the Fermi energy and holes in the Fermi sea. This seemed to be a reasonable approximation, since the configurations formed by particles in the Dirac sea and holes in the Fermi sea (ah-pairs) are more than 1.2 GeV away from the normal ph-pairs. Indeed, a proper coupling to the Dirac sea and current conservation was neglected in these investigations. They showed considerable deviations from the results obtained form time-dependent RMF-calculations with the same Lagrangian, particularly for isoscalar excitations [16]. A fully self-consistent treatment with current conservation requires the inclusion of a very large number of ah-pairs connected with a considerable numerical effort. Most of the very successful applications of RRPA theory based on non-linear meson-nucleon coupling models in the last ten years have been carried out in this way [17,18,19,20,21].
There are also relativistic continuum RPA calculations based on the nonspectral representation of the response function using the single particle Green's function in the continuum with proper boundary conditions [14]. These calculations are done for meson exchange forces with finite range. The early investigations were based on linear models. Later on the method was generalized to include non-linear coupling terms between the mesons [22]. This leads to a a more sophisticated density dependence which is crucial for a realistic description of giant resonances in nuclei [22,23].
Of course, because of the finite range of the effective force these models are relatively complicated not only for static applications to triaxially deformed or rotating nuclei, but also for investigations of nuclear dynamics, such as the solution of the relativistic RPA or linear response equations for the description of excited states. In particular one needs simpler forces for applications going beyond the mean field approach such as Particle Vibrational Coupling (PVC) [24] or configuration mixing calculations in the framework of the Generator Coordinate Method (GCM) [25]. Therefore over the years several attempts have been made to develop relativistic point coupling (PC) models with forces of zero range [26], in analogy to non-relativistic Skyrme-functionals. but only recently parameter sets have been found, which are comparable in quality to the density dependent meson-exchange models [27,28].
PC models do not contain mesonic degrees of freedom and are therefore closer to the philosophy of the density functional theory. Their essential advantage is of course the fact that the zero range of the effective interaction reduces considerably the numerical effort in practical applications. Because of their simplicity they are nowadays much used in many complex calculations going beyond the mean field approach [24,25]. However, so far they have not been used much for the dynamic investigations and it is only quite recently that a code has been developed to diagonalize the RPA equations for relativistic Point Coupling models [21] and it has been shown that this latter approach reproduces excitation and collective phenomena, in particular Giant Multipole Resonances, with a quality comparable to that of standard finite-range forces.
This manuscript is devoted to an investigation of relativistic point coupling models with an exact treatment of the coupling to the continuum. The relativistic response equations are solved both in the continuum and in the spectral representation and the corresponding results are compared. We use the Lagrangian PC-F1 [27], which is capable of reproducing a wide range of experimental data.
The paper is organized in the following way: In Sec. 2 we present the main characteristics of the point coupling RMF theory, while the relativistic RPA equations are derived in Sec. 3. The proper treatment of the continuum in connection with point coupling models is discussed in Sec. 4 and in Sec. 5 we finally present applications of this method for the spectra of in spherical nuclei. In particular we calculate the strength function of Isoscalar and Isovector Giant Resonances as well as their contributions to their respective energy weighted sum rules. The results are summarized in Sec. 6.
2 Relativistic mean field theory of zero range.
As in all the relativistic models, the nucleons are described as point like Dirac particles. In contrast to the Walecka model, however, where these particles interact by the exchange of effective mesons with finite mass, point coupling models [26] neglect mesonic degrees of freedom and consider only interactions with zero range. In principle, these models are similar to the Nambu Jona-Lasinio model [30] used extensively in hadron physics. There is, however, an important difference: in order to obtain a satisfactory description of the nuclear surface properties one needs gradient terms in the Lagrangian simulating a finite range of the interaction.
A general point-coupling effective Lagrangian is constructed to be consistent with the underlying symmetries of QCD (e.g., Lorentz covariance, gauge invariance, and chiral symmetry). It should in principle contain every possible term, allowed by these symmetries, but at the same time should also be described by the least possible number of parameters in order to give a quantitative solution.
In this work we use the point coupling Lagrangian introduced by Buervenich et al. in Ref. [27]. It presents an expansion in powers of the nucleon scalar, vector and isovector-vector densities. The Lagrangian consists of the term for free nucleons: the term for normal four-fermion interactions the term for higher order terms leading in mean field approximation to a density dependence the term containing derivative terms which simulate in a simple way the finite range of the forces: and finally the electro-magnetic part In these equations, ψ represents the nucleon spinors. The subscripts S and V are attributed to scalar and vector fields, while the subscript T is attributed to isovector fields. As usual, vectors in isospin space are denoted by arrows, where symbols in bold indicate vectors in ordinary three-dimensional coordinate space. From this Lagrangian and the corresponding energy momentum tensor we can derive a relativistic energy density functional. It has the form: where the energy density consists of a kinetic part an interaction part The interaction part depends on the local densities: and currents As in all relativistic mean field models, the no-sea approximation is used in the calculations of the nuclear densities by summing only over the singleparticle states with energies in the Fermi sea. Vacuum polarization effects are not taken into account explicitly but only in a global way by the correct choice of the Lagrangian parameters. All interactions in the Lagrangian (1) are then expressed in terms of the corresponding local densities Many effects, which go beyond mean field, seem to be neglected on the classical level, such as Fock-terms, vacuum polarization, short range Brueckner correlations etc. However, the coupling constants of the method are adjusted to experimental data, which, of course, contain all these effects and many more. Therefore these effects are not neglected. On the contrary, they are taken into account in an effective way. This concept of RMF methods is therefore equivalent to that of density functional theory.
The time-dependent variational principle allows us to derive from the energy density functional E[ρ] an equation of motion for the time-dependent relativistic single particle density: which has the form i∂ The self energy, i.e. the single particle hamiltonianĥ(ρ(t)) is obtained as the functional derivative of the energy density functional with respect to the relativistic density matrix:ĥ This yields the Dirac hamiltonian: with the self-consistent scalar and vector potentials The nucleon isoscalar-scalar, isovector-scalar, isoscalar-vector and isovectorvector self-energies are density dependent and defined by the following relations: Here we have neglected retardation effects, i.e. second derivatives with respect to the time for the various densities.
In the static limit we have thus the static densityρ 0 is obtained from the solution of the self-consistent Dirac equations upon all the nucleons with eigenvalues ε k and eigenfunctions ψ k (r):ĥ |ψ k (r) = ε k |ψ k (r) .
For spherical symmetry the spinors have the form: The subscripts n, κ and m are principal and angular momentum quantum numbers; κ = ∓(j + 1 2 ) for j = l ± 1 2 , where j and l are the total and the orbital angular momenta of the nucleon. As usual, m is the z component Coupling const. of the total angular momentum. The spherical spinors Y κ m (Ω) are given in terms of spherical harmonics Y lm l (Ω) and Pauli spinors χ ms as: while the functions f i (r) and g i (r) satisfy the static radial Dirac equations: The point coupling Lagrangian used in this work contains eleven coupling constants. Based on an extensive multi parameter χ 2 minimization procedure, Bürvenich et al. [27] have adjusted the parameter set P C-F 1 to reproduce ground state properties of infinite nuclear matter and spherical doubly closed shell nuclei. This set is listed in Table 1 and it has been tested in the calculation of many ground state properties of spherical and deformed nuclei all over the periodic table. The results are very well comparable with reasonable effective meson-exchange interactions.
The nuclear ground state is defined as the equilibrium point of the functional (7), thus, is associated with the density which minimizes E RMF [ρ]. Furthermore, small oscillations around this equilibrium point correspond to the vibrational nuclear states. They are usually described within the harmonic approximation, that is, using linear response theory. In nuclear physics, this is the so called Random Phase Approximation (RPA) which has been already mentioned in our discussion and will be described in more detail in the next section.

Relativistic RPA formalism
Under the influence of an external field F (ω) oscillating with the frequency ω the nucleus is excited. The cross section of this process is proportional to the strength function: where F αβ is the operator inducing the reaction and R αβγδ (ω) is the response function which, in an arbitrary representation indicated by the Greek indices α, β, . . . (e.g. the (r, s)-representation) is defined as: The imaginary part iη is infinitesimal and is introduced in order to fulfill the proper boundary conditions and to prevent R(ω) from diverging at ω = E ν − E 0 . We use here the response derived from the retarded Green's functions as defined in Ref. [53] In the independent particle model, |0 is the Slater determinant of the ground state, formed by the self-consistent solutions of the Dirac equation (30) and |ν = a + p a h |0 are ph-states, while E 0 and E ν are the corresponding energies. In the basis |k , where the single particle hamiltonian (22) is diagonal we obtain the free response function: with the occupation factors: The full response of Eq. (35) contains the transition densities: They can be deduced from the time-dependent density matrix in Eq. (19), which is derived from the variational principle in Eq. (18).
In the small amplitude limit one uses the linear response approximation to obtain the full response R(ω) of Eq. (35) as the solution of the linearized Bethe-Salpeter equation: The relativistic residual interaction is found as the second derivative of the energy density functional (7) with respect to the density matrix Once again, we have neglected retardation and this effective interaction has to be calculated at the static density. In a short hand notation the response equation (39) has the formal solution or introducing the inverse of R 0 we have The evaluation of the strength function (34) requires therefore three steps. The starting point is the calculation of the free response function R 0 (ω). In the next step one determines the interaction V ph and finally one solves the response equation by the inversion (41). In details there are several methods to proceed. In particular one can choose various basis sets to solve these equations.
a) As we have seen in Eq. (36) the free response has a particularly simple form in the basis of Dirac spinors (Dirac basis) diagonalizing the selfconsistent mean field equation (30). This is in particular simple for cases where the Dirac equation is solved in a discrete basis, as for instance the oscillator basis [9] or in a Saxon Woods basis [10] determined by the solution of the Dirac equation in a box with finite size. However, the simplicity in the calculation of R 0 (ω) is compensated by the computational effort required in the next steps. First we have to calculate a large number of matrix elements for the interaction (40) in the basis of the corresponding ph-states and in a second step the matrix (1 − R 0 (ω)V ph ) has to be inverted for each value of the frequency ω. In general the number of single particle states is rather large and this leads to a huge number of ph-states, requiring considerable computational sources, not only in memory but also in computer time. This is in particular a problem in the case of deformed nuclei. By this reason this method can only be used successfully for light spherical nuclei, where the number of ph-states is limited.
b) The inversion is particular simple in the RPA-basis. Inserting expression (36) into Eq. (42) we find that the response function is equivalent to the resolvent of the RPA matrix where Of course, the calculation of this matrix requires the same numerical effort as the evaluation of V ph in the Dirac basis discussed above. However there exist standard routines for the diagonalization of the RPA-matrix and this diagonalization has to be carried out only once, whereas the inversion of the response equation has to be done for each value of the frequency ω. In the RPA-basis given by the eigenvectors |µ the reduced response function defined in Eq. (52) has a particular simple form we find for R F F (ω) and for the strength function in Eq. (34) Here ∆ is a smearing parameter, which introduces a folding with a Lorentzian and is introduced by numerical reasons. c) In many cases the effective interaction V ph αβα β can formally be written as a sum of separable terms.
where Q c are single particle operators characterized by the channel index c.
As discussed in Appendix A, this is particular the case for the effective interaction of the relativistic point coupling model PC-F1 used in the present investigation. Working in the channels given by these operators Q c the numerical effort can be simplified considerably. We insert the effective interaction (51) into the Bethe-Salpeter equation (39) and introducing the reduced response function: equation (39) turns into the reduced Bethe Salpeter equation which has the same formal solution as given in Eq. (41). In all cases, where one has a continuous channel index c, as for instance the radial coordinate r, this is an integral equation. In Eq. (53) the interaction V ph c is diagonal with respect to the channel index c. This is not always the case. However, as we shall see in Appendix A, the relativistic interaction PC-F1 can be expressed to a large extent in this way. We have to allow only in specific cases also for non-diagonal interactions V ph cc , as for instance in the case of the Coulomb force or in the case of derivative terms. This is a rather simple extension of the present method and therefore, for the sake of simplicity, we will restrict ourselves in the following to an interaction diagonal in the cannel index c. If the external operator F in Eq. (34) can be expressed by the operators Q c as we finally obtain the strength function as If F cannot be expressed in terms of the operators Q c we obtain R F F from the Bethe-Salpeter equation (39) as 4 Treatment of the continuum.
As we have briefly discussed earlier, a proper treatment of the continuum is not possible by using a discrete basis, because one needs a tremendously large number of ph-states to fill up the continuum with. Instead, it can only be properly taken into account if one makes use of the more flexible linear response formalism in an appropriate channel space. Starting from Eq. (52) for the reduced response function and using Eq. (47) we derive the following expression for the reduced free response, which depends only on the energy ω and the channel indices c, c : where h stands for occupied (hole) and p for unoccupied (particle) states. It is easy to show that the sum over p can be safely extended to run over the full space, since terms of the form hh vanish due to the cancellation of forward and backward going parts. Using completeness we obtain: Here,ĥ is the Dirac hamiltonian (22) and G(E) = 1/(E −ĥ) is the corresponding single particle Green's function.
In this work we use relativistic zero range forces, thus it is appropriate to work in coordinate space. The method described in the following is a relativistic generalization of the method introduced by Bertsch et al [31] for non-relativistic zero range forces. In this case we solve the response equation in r-space, which is considerably simpler than the method introduced in Refs. [14] for finite range forces.
In coordinate representation the indices α,β, . . . in Eq. (35) are abbreviations for the "coordinates" 1 = (r 1 , d 1 , s 1 , t 1 ), where s is the spin, t the isospin coordinate, and d = 1, 2 labels large and small components. Starting from the energy density functional (7) we find the effective interaction in Eq. (40) to be of the form (51): with the local channel operators Q c (r) defined by where we distinguish the "coordinates" abbreviated by the upper index (1) and the channel index (r, c) used in Eq. (51). Due to this r-dependance, the dimension of the matrix R 0 cc (r, r ; ω) in the numerical applications will be the number of r-mesh points times eight, which represents the number of the covariant channels c, given in Table 6 of the Appendix A. This implies that all scalar, longitudinal, and transverse modes (isoscalar and isovector) are fully included and mixed by the matrix inversion of Eq. (41).
This channel index has now a continuous part given by the radial coordinate r and a discrete part characterized by the quantum numbers c = (D, S, L, T ) where the Dirac index D runs over three 2×2 matrices γ D = γ 0 , 1, γ 5 defined in Eq. (80), S = 0, 1 is the spin, L the orbital angular momentum and T = 0, 1 the isospin. Further details are given in Appendix A.
Inserting the channel operators (60) into Eqs. (52) and (58) we obtain the reduced free response function: The sum runs over all the occupied states (hole) states h with the 2-dimensional radial Dirac spinor h(r)| = (f * h (r) g * h (r)) in Eq. (33) and over all the quantum numbers κ = (lj) compatible with the selection rules in the reduced angular and isospin matrix elements where e Tc = 1 in the isoscalar channel (T c = 0) and e Tc = ±1 (for protons or neutrons) in the isovector channel (T c = 1). The reduced matrix elements of the operator σ Sc Y Lc J contain integrations over the orientation angles Ω and sums over the spin indices. The matrix elements of the form h|γ D G(E)γ D |h depend on r and r and are obtained by summing over the Dirac indices d = 1, 2 for large and small components. The Green's function G κ (r, r , E) describes the propagation of a particle with the energy E and the quantum numbers κ from r to r . It can either be calculated by spectral or non-spectral methods. In the spectral representation [32] it is obtained as a discrete sum over a complete set of eigenstates |n(r) of the radial Dirac equation (33) with the quantum number κ using box boundary conditions (or an oscillator expansion). In this case the continuum is discretized, in correspondence to the bound states inside the potential. In principle, the radial quantum number n runs over the whole single particle basis characterized by the angular quantum number κ, but one can show that this is identical to summing only over the unoccupied states, since the hole-hole pairs in Eq. (62) are not contributing, due to the cancellation between forward and backward going part. Furthermore, because of the no-sea approximation the states in the Dirac sea are empty and therefore the sum over n in Eq. (63) has also to be extended over the negative energy states. This corresponds to the sum over the ah-components discussed in the introduction. In practical applications one has to restrict this infinite set by a finite sum introducing an upper limit p − h < E ph cut in energy for the particle states p above the Fermi surface and a lower limit a − h > −E ah cut for the negative energy solutions a is introduced in order to make the -otherwise infinite -sum, tractable. This leads to a discretized spectrum.
In the spectral representation the response function R 0 (ω) has poles at the ph-energies ω = ±(ε p − ε h ) and the full response function R(ω) has poles at the eigenenergies Ω µ of the RPA-equation (46) in the same restricted space. For real frequencies ω it is purely real, and therefore the strength function vanishes everywhere apart from these poles. For complex energies ω + i∆/2, however, these poles are shifted from the real axis and one obtains a continuous spectrum, with the phenomenological width ∆. This procedure yields identical results as the diagonalization of the RPA-matrix in (46) along with a subsequent folding with a Lorentzian as discussed in Eq. (50).
In the non-spectral or continuum approach [31] the single particle Green's function is constructed at each energy from two linearly independent solutions of the Schroedinger equation with different boundary conditions at r = 0 and at r → ∞. In the relativistic case the Dirac-equation in r-space depending on the quantum number κ is a two-dimensional equation and therefore the corresponding single particle Green's function is a 2×2 matrix. Using the bracket notation of Dirac for the 2-dimensional spinors we can write [33]: where u(r) and w(r) are two independent Dirac spinors [33]: normalized in such a way that the Wronskian which is independent of r, is normalized to unity. The solution u κ (r) is regular at the origin and the solution w κ (r) fulfills outgoing wave boundary conditions [34]. Further details are given in Appendix B.
Provided that the free response function R 0 c,c (r, r ; ω) has been properly derived, we are able to solve the reduced Bethe-Salpeter equation (53) R c,c (r, r ; ω) = R 0 c,c (r, r ; ω) where the index c runs over the various discrete channels given in Table 6. Finally the strength function is obtained as: The sum rules are defined as moments of the strength function S(ω): They are helpful to characterize the spectral distribution of the oscillator strength. In particular they allow us to define the centroid energy by the ratio This quantity can be compared directly with experimental values. Of course, in most experiments only a restricted energy range is accessible and therefore one also has to restrict the integration in Eq. (69) to the same energy window.
Other important quantities are transition densities in various channels c with respect to the operator F as for instance the neutron and proton transition densities:

Applications
In the previous section we briefly described how conventional RPA methods treat the continuum part of the spectrum through the introduction of a potential "wall" far from the nucleus. In the credit side of this approach, general properties of collective excitations can be very well reproduced, either by using finite range or point coupling interactions (Nikšić et.al. [21]). Since CRPA can treat the coupling to the continuum exactly, it is of interest to see how well this model does in reproducing the properties of excited state in finite nuclei, in particular the giant resonances.
The most prominent resonances are the Isoscalar Monopole Resonance (ISGMR), which is a breathing of the nucleus as a whole, the Isovector Dipole Resonance (IVGDR) which corresponds to a collective excitation of the proton against the neutron density, and Isoscalar Quadrupole Resonance (ISGQR). In addition we have the Isoscalar Dipole Resonance (ISGDR) revealing the spurious state corresponding to a translational motion of the nucleus. These modes show up in an energy range of 10 − 30 MeV and they exhaust a major portion of the corresponding sum rules. In the next sections we investigate the ISGMR, the IVGDR and the ISGDR in more detail.

Numerical details
In the following, we perform several calculations using the relativistic continuum RPA approach in r-space with Point Coupling forces [27]. We select the doubly magic nuclei 16 O, 40 Ca, 132 Sn and 208 Pb to investigate how the collective excitation phenomena depend on an exact coupling to the continuum.
In a first step, the ground state of the nucleus is determined by solving the self-consistent RMF equations (33) for the parameter set PC-F1 given in Table 1. The method we are using is a fourth order Runge-Kutta in r-space (Dirac-mesh) where nucleons move in a spherical box with radius R D = 15 fm and with a mesh size d D = 0.05 fm.
Using the single particle wave functions and the corresponding energies of this static solution, we determine the free response R 0 of Eq. (62) in the same box radius but using a wider mesh in r-space (response-mesh). The size d R of this mesh depends on the excitation mode; for the monopole modes we use d R = 0.15 fm, while for the dipole a larger interval d R = 0.30 fm is sufficient. Then we solve the Bethe salpeter equation (68) to get the strength distribution S(ω).
At the same time, we perform similar calculations using the discrete RPA approach, where the continuum is not treated exactly, aiming of course to a more precise comparison with the CRPA results. For those calculations, an energy cut-off is necessary, so that a feasible diagonalization is achieved. In particular, we have used an energy cut-off | p − h | < E ph cut = 300 MeV for the configurations with particles above the Fermi sea and | a − h | < E ah cut = 1500 MeV for configurations with anti-particles in the Dirac sea.

Isoscalar Giant Monopole Resonances
Results for the isoscalar monopole strength distribution are attainable, once the corresponding external field is used. In this case, the classical energy weighted sum rule m 1 (E0) becomes: The doubly magic spherical nucleus 208 Pb is a particularly good example in perform our calculations, since it has been used in the literature to test numerous nuclear structure models in the past, in particular applications of the random phase approximation [35,22,23,36].
In Fig. 1 we show the ISGMR strength distribution obtained by continuum RPA (full red line) and compare it with the discrete B(E0) values (blue) obtained by the spectral representation of the response function for the same parameter set PC-F1 [27].
Using the CRPA approach, we find for the calculated centroid energy defined in Eq. (70) that m 1 /m 0 = 14.40 MeV, which is rather close to the result m 1 /m 0 = 14.17 MeV deduced from discrete RPA calculations as well as to the experimental value m 1 /m 0 = 13.96 ± 0.2 MeV [37].
In those two methods, no additional smearing ∆ = 0 has been used. This means that the observed width of the continuum RPA strength corresponds entirely to the escape width which in the Pb region is very small, due to the relatively high Coulomb and centrifugal barriers in this heavy nucleus. In contrast, discrete RPA provides no width at all. Otherwise, the agreement of these two methods in this nucleus is excellent. In the panel (b) of Fig. 1, we give the neutron and proton transition densities at the peak energy, as it is calculated in Eq. (72). They emphasize the collective character of the isoscalar breathing mode extended over the entire interior of the nucleus with neutrons and protons always in phase.
In addition, the energy weighted sum rule obtained in CRPA using Eq. ]. This shows that the results obtained in the literature by relativistic RPA calculations using the spectral method are very reliable for such heavy nuclei [17,8].
In Fig. 2 we show the E0 strength distributions for the lighter doubly magic nuclei 16 O, 40 Ca, and 132 Sn. As in Fig. 1, the smearing parameter ∆ is zero, but now the escape width is considerably larger for these nuclei. Fig. 3 summarizes the results for the isoscalar monopole strength distributions as a function of the mass number A. In panel (a), we plot the centroid energies of both continuum RPA (red dots) and discrete RPA (blue dots), together with the experimental centroid energies taken from Ref. [37]. We also show the phenomenological A-dependenceĒ 1 − ≈ 31.2 A −1/3 + 20.6 A −1/6 by the dashed line. It becomes clear that CRPA can successfully reproduce collective excitations over the known range of nuclei. In panel (b) of Fig. 3 we show the escape width Γ ↑ of E0 resonances. The red values correspond to the full width half maximum (FWHM) of the peak, using continuum RPA , while the experimental values are indicated in black. The evident disagreement is not surprising, if we consider that only 1p1h-configurations are taken into account, i.e. the major part of the width resulting from the coupling to more complicated configurations such as 2p2h etc. is not described well in this simple RPA approach. It has been shown in recent investigations of the coupling to complex configurations within the framework of the relativistic time-blocking approximation (RTBA) or the relativistic quasiparticle-time-blocking approximation (RQTBA) [24] that such couplings can be taken into account successfully in a fully consistent way starting from one density functional E[ρ]. So far, relativistic investigations of this type have been carried out with discrete methods. At present, investigations in this direction including the continuum properly go beyond the scope of this paper.

Isovector Giant Dipole Resonances
Isovector Giant Dipole resonance is the most well studied collective excitation and the first to be observed experimentally [38]. An external electromagnetic field of the form: causes protons and neutrons to oscillate in opposite phases to each other and this leads to a pronounced peak in the photoabsorption cross section. This mode has been well studied in many nuclei [39].
With the increasing number of experiments in systems far from stability and systems with large neutron excess, one has been able to observe also lowlying E1 strength in the area of the neutron emission threshold. It is called Pygmy Dipole Resonance PDR and can be interpreted as a collective mode with dipole character where the neutron skin oscillates against an isospin saturated proton-neutron core. This mode has first been predicted in phenomenological models [40] exhausting several percent of the electric dipole sum rule. In recent years, it has been intensively investigated both on the experimental side by the Darmstadt group [41,42] as well as on the theoretical side, using discrete relativistic RPA calculations based on NL3 [43].
In Fig. 4  In addition to the giant dipole resonance a smaller peak appears at the energy region of the neutron emission threshold around E ∼ 7.5 MeV, that corresponds to the pygmy resonance.
In panel (b) of Fig. 4 we give the transition densities associated the lowlying peak at E = 7.66 MeV and the GDR peak at E = 12.9 MeV. The higher peak has clearly an isovector character, since the neutrons are oscillating against the protons over a large radial range centered at the surface. The lower peak shows an isoscalar core, where neutrons and protons oscillate in phase and a pure neutron skin moving against the T = 0 core. This is the typical behavior of the pygmy mode.
Closer investigation of pygmy resonances have shown that this mode is in the neighborhood of the neutron separation threshold, slightly below for small and slightly above for large neutron excess (see for instance Ref. [20]). It is therefore of particular importance to study this mode with a proper treatment of the continuum, since in most of the previous investigations this has not been possible. We show in Fig. 5    In the same figure and below the threshold we have in both cases discrete lines. The solid blue ones are again the eigen-solutions of the DRPA-equation (46). The solutions of the CRPA equations lead in this region also to discrete poles. We show them by dashed red lines at the pole of the full response function. Numerically, the only way to determine the B(E1) values of these poles in CRPA is by using very small imaginary parts ∆ → 0 in the frequency ω + i 1 2 ∆ and then determining the B(E1) values by simple integration over a small interval around this pole.
By doing that, we finally observe that there are differences in the details between the continuum and the discrete RPA calculations close to the neutron separation threshold. In Table 2   In addition, we find that the escape width in this nucleus is considerably smaller in the E1 channel as compared to the E0 channel in Fig. 2. This has the following explanation: The selection rules for ph-excitations with E0 character is ∆j = 0 and no change in parity. It turns out that most of the ph-excitations contributing to the strong peak in the resonance region have rather small values for the particle configurations and therefore a very low or no centrifugal barrier. This is different for the E1 resonance, where one has a change in parity and in addition changes of ∆j = 0, ±1. In such a case, a large part of the contributing ph-pairs have particles with larger -values i.e. a strong centrifugal barrier and hence the width becomes smaller.     In Fig. 7 we show the region of the PDR in the doubly magic nucleus 132 Sn. As already found in Ref. [20], the theoretical neutron emission threshold at E = 7.13 MeV lies much below the area of interest. As before, we calculate the B(E1) values of the prominent peaks, for both discrete and continuum calculations with the total strength to be in good agreement. In Table 3 we show in what extent each level contributes to the total pygmy collective state. Finally, the energy weighted strength m 1 in this area is 13.24 [e 2 fm 2 ] (i.e. 2.35 % of the total sum rule) for CRPA and 20.45 [e 2 fm 2 ] (i.e. 3.46 % of the total sum rule) for DRPA.

CRPA
In Fig. 8 we show the electric dipole strength distribution of the lighter nuclei 16 O and 40 Ca. The strength obtained in CRPA calculations (red curves) are compared with the B(E1)-values resulting from discrete DRPA calculations (blue lines). The position of the corresponding peaks and poles with large strength are in rather good agreement, as explained in Table 4. We find, however, that in the continuum calculations a much larger escape width emerges, in particular for the nucleus 16 O.

Isoscalar Giant Dipole Resonances
Besides the distribution of the isovector dipole strength which is dominated by the IVGDR in many experimental spectra, in recent years there has also been considerable interest in measuring the isoscalar dipole strength distribution [47,48]. In a similar way, one expects to find the ISGDR, which corresponds to a compression wave going through the nucleus along a definite direction and to learn from such experiments more about the nuclear incompressibility. Relativistic calculations based on discrete RPA [17,23] have shown that the resonance energy of this mode is indeed closely connected to the incompressibility of nuclear matter.
Along with this ISGDR resonance built on 3 ω-excitations above 20 MeV, calculations based on both relativistic [17] and non-relativistic [49] RPA approaches have revealed a low-lying isoscalar dipole strength in the region below and around 10 MeV. Experimental investigations with inelastic scattering of α-particles at small angles [50,48] have also found isoscalar dipole strength in this region. This strength has been attributed in Ref. [19] to an exotic mode of a toroidal motion predicted already in early theoretical investigations on multipole expansions of systems with currents [51] and investigated also by semiclassical methods [52] On the theoretical point of view, there is further interest in the isoscalar dipole mode, characterized by the quantum numbers (J π = 1 − , T = 0), because it contains the Goldstone mode connected with the violation of translational symmetry in the mean field solutions. This mode corresponds to the center of mass motion of the entire nucleus. Because of the missing restoring force, this mode has vanishing excitation energy. It is one of the essential advantages of the RPA approximation, that it preserves translational symmetry and therefore it has an eigenvalue at zero energy with the eigenfunction given by the ph-matrix elements of the linear momentum operator.
Since the ISGDR is expected to be a 3 ω-excitation it is usually associated with the external field derived in Ref. [54] where the factor η = 5 3 r 2 is used to extract the spurious center of mass motion.
In the upper part of Fig. 9 we display the distribution of the isoscalar dipole strength in 208 Pb, calculated with the operator (77) for η = 0, that is, we take no action for the spurious state. We therefore observe a huge peak close to zero energy, which dominates the spectrum and corresponds to the spurious translational mode.
It turns out that the position of this spurious state is an extremely sensitive object which strongly depends on the numerics of the model. Of course the optimal would be to calculate the spurious state at exactly zero energy. Therefore this excitation mode presents an ideal benchmark for numerical shown that the exact separation of the spurious state requires a fully selfconsistent solution [22]; a fact which was not given in most of the older applications with Skyrme or Gogny forces. In many cases, only few of the different terms in the residual interaction had been taken into account in RPA calculations. In addition, the configuration space must be full. Indeed, the discussed drawback of the conventional spectral representation in a truncated phconfiguration space affects the position of the spurious state. Therefore, the convergence to zero eigenvalue of the spurious translational mode occurs very slowly and only in extremely large configuration space. In relativistic applications this is translated to including also large spectrum in the Dirac sea [15,8]. As a consequence, in the spectral representation, one has to take into account many configuration with particles in the Dirac and holes in the inevitably decreases the efficiency of the method. Fortunately, using the continuum RPA approach, one is free from such constraints and limitations, since the entire configuration space is automatically included. The results in Fig. 9 obtained with the operator (77) for η = 0 show clearly the spurious state dominating the entire spectrum (see the scale). Its position is not precisely at zero energy, rather it depends on the mesh size used for the solution of the continuum response equation (the response mesh). In panel (a) of Fig. 9 we present two calculations with different mesh-sizes, where in panel (b) we show how the spurious state moves to zero energy as we use a finer radial interval. For the ideal case of an infinitesimal mesh, the strength connected with the spurious state would be completely separated from the rest of the spectrum.
In Fig. 10 we show results obtained with the full operator (77), i.e. with η = 5 3 r 2 , in a scale increased by three orders of magnitude. Obviously this procedure removes the spurious state with high precision. We also did not observe any influence of the isoscalar mode in the isovector channel due to isospin mixing. In this context we have to remember, that the isospin mixing introduced on the mean field level is corrected on the RPA level to a large extend [12].
The main part of the remaining isoscalar dipole spectrum in Fig. 10 Table 5: Self-consistent (relativistic and non-relativistic) RPA calculations performed for the ISGDR in 208 Pb, compared with the most recent experimental data. The two columns refer to the centroid energies of both the lowand high-energy sides of the ISGDR mode.
located at E ≈ 25 MeV. This "exotic" mode is best described as a "hydrodynamical density oscillation", in which the volume of the nucleus remains constant and the state can be visualized as a compression wave oscillating back and forth through the nucleus [19].
Moreover, Fig. 10 shows an additional mode in the region of 10 − 15 MeV that exhausts roughly 20% of the total sum rule. This peak does not correspond to a compression mode, but as discussed in Ref. [19] rather to a kind of toroidal motion. The toroidal dipole mode is understood as a transverse zero-sound wave and its experimental observation would invalidate the hydrodynamical picture of the nuclear medium, since there is no restoring force for such modes in an ideal fluid.
In conclusion, continuum RPA calculations manage not only to predict the existence of the toroidal and the compression mode, but also to achieve a reasonable agreement of the corresponding centroid energies to other models focusing on the same problem, as well as to recent experimental data [48,47]. In Table 5, these results are presented for the case of the well studied nucleus 208 Pb.

Conclusions
Starting from a point coupling Lagrangian, we have used the non-spectral relativistic RPA approach to examine the corresponding excitation spectra and we have compared the results with spectral calculations based on the same Lagrangian. This non-spectral method has several advantages. The coupling to the continuum is treated consistently using the relativistic single particle Green's function at the appropriate energy. In this way, complicated sums over unoccupied states are avoided. This is particularly important for relativistic applications since the Dirac sea is now automatically treated properly and the unphysical transitions from holes in the Fermi sea to particles in the Dirac sea is avoided as long as we restrict our investigations to positive energies.
The ground state phenomena are calculated using the same Lagrangian by a self-consistent solution of the relativistic mean field equations in r-space. The residual particle-hole interaction used in the RPA calculations is derived in a fully self-consistent way from the second derivative of the corresponding energy density functional. In this way no additional parameters are required and one is able to reproduce the collective properties, namely the multipole giant resonances for various doubly closed shell spherical nuclei over the entire periodic table.
The calculations are carried out by using a new relativistic continuum RPA program for point-coupling models, that includes all the terms in the Lagrangian, in particular the two-body interactions with zero range, the density dependent parts with all the rearrangement terms, the derivative terms, the various current-current terms and the Coulomb interaction. As applications the nuclei 16 O, 40 Ca, 90 Zr, 132 Sn, and 208 Pb have been investigated proving that a hight level of accuracy is achieved, as compared to the discrete methods. Comparing calculations with spectral and non-spectral representations of the response function for the same Lagrangian we find, that in general the spectra are well reproduced within the spectral approximation, if an appropriate phenomenological smearing parameter is used and if a sufficiently large number of ph-configurations is taken into account in the latter case. We find, however, differences in neighborhood of the neutron threshold, where the coupling to the continuum is not properly reproduced in the spectral method.
As compared to the discrete case the non-spectral representation has the advantage of (i) a precise treatment of the coupling to the continuum and a fully consistent determination of the escape width without a phenomenological smearing parameter, (ii) a faster evaluation of the cross section, because one needs for fixed energy only two scattering solutions instead of the thousands of ph-configurations in the discrete case and (iii) a proper treatment of the Dirac sea without any further ah-configurations.
Relativistic CRPA describes very well the position of resonances in doubly magic spherical nuclei. Provided that proper pairing correlations are taken into account, a similar method can also be applied in open-shell nuclei. This requires the development of the relativistic continuum quasiparticle random phase approximation (CQRPA). This approach accounts on equal footing for the influence of the residual particle-hole (ph) as well as the particle-particle (pp) correlations. In analogy to non-relativistic calculations [58,59,60,61] this can be achieved on the basis of relativistic CRPA theory developed in this manuscript either by treating the pairing correlations in the BCS approach for nuclei far from the drip lines where no level in the continuum is occupied, or in the Hartree-Bogoliubov approximation valid for all nuclei up to the drip line. Investigations in this direction are in progress.
Of course, the present approach is based on the RPA and includes only 1p1h-configurations. Therefore only the escape width of the resonances can be reproduced properly. For heavy nuclei the decay width resulting from a coupling to more complex configurations is very important. In fact, such couplings have been introduced successfully in the relativistic scheme using the spectral representation in Refs. [24]. On the non-relativistic side, such techniques have also been used in the context of the non-spectral representation without [62,63] and with [64] pairing. So far, however, fully self-consistent relativistic applications including complex configurations with a proper treatment of the continuum are still missing.
Helpful discussions with G. Lalazissis, E. Litvinova, T. Nikšić, N. Paar, V. Tselyaev, and D. Vretenar are gratefully acknowledged. This research has been supported the Gesellschaft für Schwerionenforschung (GSI), Darmstadt, the Bundesministerium für Bildung und Forschung, Germany under project 06 MT 246 and by the DFG cluster of excellence "Origin and Structure of the Universe" (www.universe-cluster.de).
A The effective interaction in density dependent point-coupling models In Eq. (40) the effective interaction for RPA calculations is defined as the second derivative of the energy functional with respect to the density matrix: In coordinate representation the indices α,β, . . . are an abbreviation for the "coordinates" 1 = (r 1 , s 1 , d 1 , t 1 ), where s is the spin and t the isospin coordinate, and d = 1, 2 is the Dirac-index for large and small components. Starting from the energy density functional (7) and neglecting for the moment the Coulomb force, we find the density dependent zero range force where the vertices Γ c are 8×8 matrices acting on the indices s, d, t and reflect the different covariant structures of the fields including spin and isospin degrees of freedom. We express the 4×4 Dirac matrices as a direct product of spin matrices σ and 2×2 matrices γ D acting on large and small components and the spin matrices σ S=0 = 1 and σ S=1 = σ µ with the spherical coordinates of the Pauli spin matrices. In this way we obtain the vertices Γ c = γ D ×σ S ×τ T as direct products of 2-dimensional Dirac-, spin-and isospin matrices (see also the second column of Table 6). Finally, in Eq. (79) the quantities υ c (r) describe the strengths of all the various parts of the interaction derived in a consistent way from the Lagrangian. The ones derived from the four-fermion terms (3) are constants. Furthermore, due to a density dependence of the higher order terms (4) as well as the corresponding rearrangement terms, υ c (r) depends on the static density and therefore on the coordinate r. In addition, because of the derivative terms (5), they also contain Laplace operators. Summarizing, we have: c υ c (r) = scalar: α S + 2β S ρ S (r) + 3γ S ρ 2 S (r) + δ S ∆ time-like vector: In the isovector case the constants α S , α V , δ S and δ V are replaced by α T S , α T V , δ T S and δ T V . As we see in Table 1 the corresponding values β T S = γ T S = γ T V vanish. For spherical nuclei, the densities and currents in the Lagrangian depend only on the radial coordinate r. Therefore we expand the δ-function in Eq. (79) in terms of spherical harmonics Combining spin (S) and orbital (L) degrees of freedom we find by re-coupling to total angular momentum J Inserting this expression into Eq. (81) we obtain for the interaction a sum (or integral) of separable terms (channels) Each channel is characterized by a continuous parameter r and the discrete numbers c = (D, S, L, J, T ). The corresponding channel operators Q (1) c (r) are local single particle operators and the upper indices (1) and (2) in Eq. (84) indicate that these operators act on the "coordinates" 1 = (r 1 Ω 1 s 1 d 1 t 1 ) and 2 = (r 2 Ω 2 s 2 d 2 t 2 ). The total angular momentum is a good quantum number and for fixed J the sum over c in Eq. 84 runs only over specific numbers c = (D, S, L, T ) determined by the selection rules. We concentrate in this manuscript on states with natural parity, i.e. π = (−) L = (−) J . Considering that S = 0 for the scalar and the time-like vector and that S = 1 for the space-like vector we therefore have Finally we have eight discrete channels. Their quantum numbers are shown in Table 6. An essential feature of the effective interaction (81) is that it contains derivative terms in the form of Laplacians ∆ (retardation effects are neglected). In spherical coordinates, they contain radial derivatives as well as angular derivatives. The latter can be expressed by the angular momentum operators acting on spherical harmonics Y LM . Therefore we obtain: Here the radial derivatives ← − ∂ r and − → ∂ r act on the right and on the left side in Eq. (67), i.e. on R 0 c c (r r) and on R cc (r, r ). Since the integration is discretized r → r n = nh the operator − → ∂ r is represented by a matrix in r-space as for instance by the tree-point formula: This means that the term υ c (r) in Eq. (68) is no more diagonal in the coordinate r and it must be replaced by a matrix υ c (r, r ). The term which leads to off-diagonal terms in channel space is the Coulomb interaction. It brakes isospin symmetry and therefore it will be described by the general form υ cc (r, r ). In particular, we will have and the r dependance can be written as: and r < and r > are the smaller and the greater of r and r . This leads to a matrix υ cc (r, r ) in Eq. (68) as shown in Table 7.

B The continuum representation for the Green's function
In a non-spectral or continuum approach the relativistic single particle Green's function G κ (r, r ; E) obeys the equation: E −ĥ κ (r) G κ (r, r ; E) = δ(r − r ), whereĥ κ (r) is the radial Dirac-operator of Eq. (33) depending on the quantum number κ = (lj). This Green's function can be constructed at each energy E from two linearly independent solutions |u(r) = f u (r) g u (r) , |w(r) = f w (r) g w (r) (92) u * (r)| = (f u (r) g u (r)), w * (r)| = (f w (r) g w (r)) (93) of the Dirac equation with the same energy E E −ĥ κ (r) |u(r) = 0, E −ĥ κ (r) |w(r) = 0, but with different boundary conditions. The functions u(r) and w(r) are normalized in such a way that the Wronskian is equal to: Of course these scattering solutions depend on the energy E and on the quantum number κ, i.e. we have |u κ (r, E) and |w κ (r, E) . The Dirac-equation in r-space is a two-dimensional equation and therefore the corresponding single particle Green's function is a 2×2 matrix. Using the bracket notation of Dirac for the 2-dimensional spinors and following Ref. [33] we can express this Green's function as: G κ (r, r ; E) = |w κ (r; E) u * κ (r ; E)| for r > r |u κ (r; E) w * κ (r ; E)| for r < r with G κ (r , r; E) = G κ (r, r ; E) The solution u κ (r) is regular at the origin, i.e. following Ref. [34] we have for E > V + S in the limit r → 0: with k 2 = (E − V − S)(E − V + S + 2m) > 0 and j l (z) is a spherical Bessel function of the first kind. The wave function w κ (r) represents at large distances for E > 0 an outgoing wave, i.e. we have for r → ∞ w(r) → rh where h (1) l (z) is the spherical Hankel function of the first kind and for E < 0 an exponentially decaying state, i.e. we have for r → ∞ w(r) → r 2Kr π K l+ 1 2 (Kr) −Kr E+2m 2Kr π Kl + 1 2 (Kr) where K 2 = (V − S − E)(E − V + S + 2m) > 0 and j l (z) and K l+1/2 (z) are modified spherical Bessel functions [65]. For E < 0 the two scattering solutions are both real. This absence of any imaginary term will eventually give no contribution to the cross section of Eq. (34). We have to keep in mind, however, that at energies that correspond to eigen energies of a bound state, the solutions u κ (r, E) and w κ (r, E) coincide up to a factor, which means that the Wronskian vanishes at this energy. This corresponds to a pole in the response function on the real energy axis. By adding a small imaginary part to the energy E → E + i∆ we obtain a sharp peak in the strength distribution.

C The free response function in r-space
The reduced free response depends on the energy E and the channel indices c, c . The operators Q c given by Eq. (85) are characterized by the channel index c = (r, DSLT ). Each single particle matrix element of the form p|Q c |h in Eq. (101) separates into an angular, an isospin and a radial part.
Since we consider in this paper only ph-RPA in the same nucleus, the particle states have the same isospin as the hole states and thus the isospin matrix element p|τ T |h is simply a phase ±1.
Considering that this channel operator has a δ-function in the radial coordinate, the radial matrix elements p|γ c |h r = p(r)|γ D |h(r) then depend on r. They are found as sums over the large and small components in the radial spinors |h(r) and |p(r) for fixed values of r.
The angular matrix elements depend on the quantum numbers κ of particle and hole states, and, of course, on the channel quantum numbers S and L. In particular, we find for S = 0: Using for the angular and isospin part the abbreviation we obtain for the reduced response function of Eq. (57) in r-space: R 0 cc (r, r ; ω) = ph Q * c ph Q c ph h|γ + c |p r p|γ c |h r ω − ε p + ε h − Q * c hp Q c hp h|γ c |p r p|γ + c |h r ω + ε p − ε h (106) As in Eq. (58) we extend the sum over p over the full space and use completeness in the radial wave functions: R 0 cc (r, r ; ω) = hκ Q * c κh Q c κh h(r)|γ + c G κ (r, r ; ω + ε h )γ c |h(r ) + Q * c hκ Q c hκ h(r )|γ c G κ (r , r; −ω + ε h )γ + c |h(r) . (107) Since the angular matrix elements depend only on the quantum numbers κ the sum over p is here replaced by a sum over the quantum numbers κ, which is restricted by the selection rules of the reduced matrix elements (105). Having the exact form of the Green's function for the static radial Dirac equation (33), one can finally construct the non-spectral or continuum reduced response function (62): − Q * c hκ Q c hκ γ c hw (r; ω − ε h )γ c uh (r ; ω − ε h ) for r > r Using Eq. (97) we find R 0 c c (r , r; ω) = R 0 cc (r, r ; ω) It becomes clear now that the undeniable advantage of the non-spectral approach as compared to the spectral one, is the fact that the sum over the unoccupied states (particle states) is replaced by a sum over the quantum number κ, which is restricted by the selection rules for the reduced matrix elements Q c κh . For each κ, one has to determine only the pairs of the scattering wave functions |u and |w for the forward and backward term. In particular the sum over κ does not have to be extended over the states in the Dirac sea as in the spectral representation (for details see Ref. [8]). Therefore, not only the size of the configuration space is significantly reduced, but, more notably, the particle-hole as well as the antiparticle-hole basis is taken into account fully and without any approximation.