text:
# Quasiparticle random phase approximation based on the relativistic Hartree-Bogoliubov model

## I Introduction

The multipole response of unstable nuclei far from the line of $`\beta `$-stability presents a very active field of research, both experimental and theoretical. These nuclei are characterized by unique structure properties: the weak binding of the outermost nucleons and the effects of the coupling between bound states and the particle continuum. On the neutron rich side, in particular, the modification of the effective nuclear potential leads to the formation of nuclei with very diffuse neutron densities, to the occurrence of the neutron skin and halo structures. These phenomena will also affect the multipole response of unstable nuclei, in particular the electric dipole and quadrupole excitations, and new modes of excitations might arise in nuclei near the drip line.

A quantitative description of ground-states and properties of excited states in nuclei characterized by the closeness of the Fermi surface to the particle continuum, necessitates a unified description of mean-field and pairing correlations, as for example in the framework of the Hartree-Fock-Bogoliubov (HFB) theory. In order to describe transitions to low-lying excited states in weakly bound nuclei, in particular, the two-quasiparticle configuration space must include states with both nucleons in the discrete bound levels, states with one nucleon in a bound level and one nucleon in the continuum, and also states with both nucleons in the continuum. This cannot be accomplished in the framework of the BCS approximation, since the BCS scheme does not provide a correct description of the scattering of nucleonic pairs from bound states to the positive energy particle continuum. Collective low-lying excited states in weakly bound nuclei are best described by the quasiparticle random phase approximation (QRPA) based on the HFB framework. The HFB based QRPA has been investigated in a number of recent theoretical studies. In Ref. Eng.99  a fully self-consistent QRPA has been formulated in the HFB canonical single-particle basis. The Hartree-Fock-Bogoliubov formalism in coordinate state representation has also been used as a basis for the continuum linear response theory Mat.01 ; Mat.02 . In Ref. Kha.02  the HFB energy functional has been used to derive the continuum QRPA response function in coordinate space. HFB based continuum QRPA calculations have been performed for the low-lying excited states and giant resonances, as well as for the $`\beta `$ decay rates in neutron rich nuclei.

In this work we formulate the relativistic QRPA in the canonical single-nucleon basis of the relativistic Hartree-Bogoliubov (RHB) model. The RHB model is based on the relativistic mean-field theory and on the Hartree-Fock-Bogoliubov framework. It has been very successfully applied in the description of a variety of nuclear structure phenomena, not only in nuclei along the valley of $`\beta `$-stability, but also in exotic nuclei with extreme isospin values and close to the particle drip lines. Another relativistic model, the relativistic random phase approximation (RRPA), has been recently employed in quantitative analyses of collective excitations in nuclei. Two points are essential for the successful application of the RRPA in the description of dynamical properties of finite nuclei: (i) the use of effective Lagrangians with nonlinear self-interaction terms, and (ii) the fully consistent treatment of the Dirac sea of negative energy states.

The RRPA with nonlinear meson interaction terms, and with a configuration space that includes the Dirac sea of negative-energy state, has been very successfully employed in studies of nuclear compressional modes Vre.00 ; Pie.01 ; MGW.01 , of multipole giant resonances and of low-lying collective states in spherical nuclei Mawa.01 , of the evolution of the low-lying isovector dipole response in nuclei with a large neutron excess Vrepyg1.01 ; Vrepyg2.01 , and of toroidal dipole resonances Vre.02 .

In Sec. II we present the formalism and formulate the matrix equations of the relativistic QRPA in the canonical basis of the relativistic Hartree-Bogoliubov (RHB) framework for spherical even-even nuclei. In Sec. III the RHB+RQRPA approach is tested in the example of the isoscalar monopole, isovector dipole and isoscalar quadrupole response of 22O, and the results are compared with recent non-relativistic QRPA calculations of the multipole response of neutron rich oxygen isotopes. In Sec. IV the RQRPA framework is applied in the analysis of the evolution of the low-lying isovector dipole strength in Sn isotopes and N=82 isotones. The results are compared with recent experimental data. Section V contains the summary and the conclusions.

## II The relativistic quasiparticle random phase approximation

In this section the matrix equations of the relativistic quasiparticle random phase approximation (RQRPA) are formulated in the canonical basis of the relativistic Hartree-Bogoliubov (RHB) framework for spherical even-even nuclei.

### II.1 The relativistic mean-field Lagrangian and the equations of motion

The nuclear matter equation of state and detailed properties of finite nuclei have been very successfully described by relativistic mean-field (RMF) models Ser.86 ; Rei.89 ; Rin.96 . In this framework the nucleus is described as a system of Dirac nucleons that interact in a relativistic covariant manner by meson exchange. In particular, the isoscalar scalar $`\sigma `$-meson, the isoscalar vector $`\omega `$-meson, and the isovector vector $`\rho `$-meson build the minimal set of meson fields that is necessary for a quantitative description of bulk and single-particle nuclear properties. The model is defined by the Lagrangian density

$$=_N+_m+_{int}.$$
(1)
$`_N`$ denotes the Lagrangian of the free nucleon

$$_N=\overline{\psi }\left(i\gamma ^\mu _\mu m\right)\psi ,$$
(2)
where $`m`$ is the bare nucleon mass and $`\psi `$ denotes the Dirac spinor. $`_m`$ is the Lagrangian of the free meson fields and the electromagnetic field

$`_m`$ $`=`$ $`{\displaystyle \frac{1}{2}}_\mu \sigma ^\mu \sigma {\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2{\displaystyle \frac{1}{4}}\mathrm{\Omega }_{\mu \nu }\mathrm{\Omega }^{\mu \nu }+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _\mu \omega ^\mu `$ (3)
$`{\displaystyle \frac{1}{4}}\stackrel{}{R}_{\mu \nu }\stackrel{}{R}^{\mu \nu }+{\displaystyle \frac{1}{2}}m_\rho ^2\stackrel{}{\rho }_\mu \stackrel{}{\rho }^\mu {\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu },`$
with the corresponding masses $`m_\sigma `$, $`m_\omega `$, $`m_\rho `$, and $`\mathrm{\Omega }_{\mu \nu }`$, $`\stackrel{}{R}_{\mu \nu }`$, $`F_{\mu \nu }`$ are field tensors

$$\begin{array}{ccc}\hfill \mathrm{\Omega }_{\mu \nu }& =& _\mu \omega _\nu _\nu \omega _\mu \hfill \\ \hfill \stackrel{}{R}_{\mu \nu }& =& _\mu \stackrel{}{\rho }_\nu _\nu \stackrel{}{\rho }_\mu \hfill \\ \hfill F_{\mu \nu }& =& _\mu A_\nu _\nu A_\mu .\hfill \end{array}$$
(4)
The model Lagrangian density contains also the interaction terms

$$_{int}=\overline{\psi }\mathrm{\Gamma }_\sigma \sigma \psi \overline{\psi }\mathrm{\Gamma }_\omega ^\mu \omega _\mu \psi \overline{\psi }\stackrel{}{\mathrm{\Gamma }}_\rho ^\mu \stackrel{}{\rho }_\mu \psi \overline{\psi }\mathrm{\Gamma }_e^\mu A_\mu \psi .$$
(5)
The vertices read

$$\mathrm{\Gamma }_\sigma =g_\sigma ,\mathrm{\Gamma }_\omega ^\mu =g_\omega \gamma ^\mu ,\stackrel{}{\mathrm{\Gamma }}_\rho ^\mu =g_\rho \stackrel{}{\tau }\gamma ^\mu ,\mathrm{\Gamma }_e^m=e\frac{1\tau _3}{2}\gamma ^\mu ,$$
(6)
with the coupling constants $`g_\sigma `$, $`g_\omega `$, $`g_\rho `$ and $`e`$. This simple linear model, however, does not provide a quantitative description of complex nuclear systems. An effective density dependence has been introduced Bog.77  by replacing the quadratic $`\sigma `$-potential $`{\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2`$ with a quartic potential $`U(\sigma )={\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2+{\displaystyle \frac{g_2}{3}}\sigma ^3+{\displaystyle \frac{g_3}{4}}\sigma ^4`$. This potential includes the nonlinear $`\sigma `$ self-interactions with two additional parameters $`g_2`$ and $`g_3`$. It has been shown that one can describe the properties of nuclear matter and finite nuclei with high accuracy using density dependent coupling constants $`g_m(\rho )`$, instead of nonlinear $`\sigma `$ self-interaction TDP.02 .

From the model Lagrangian density the classical variation principle leads to the equations of motion. The time-dependent Dirac equation for the nucleon reads

$$\left[\gamma ^\mu (i_\mu +V_\mu )+m+S\right]\psi =0.$$
(7)
If one neglects retardation effects for the meson fields, a self-consistent solution is obtained when the time-dependent mean-field potentials

$`S(𝐫,t)`$ $`=`$ $`g_\sigma \sigma (𝐫,t),`$
$`V_\mu (𝐫,t)`$ $`=`$ $`g_\omega \omega _\mu (𝐫,t)+g_\rho \stackrel{}{\tau }\stackrel{}{\rho }_\mu (𝐫,t)+eA_\mu (𝐫,t){\displaystyle \frac{(1\tau _3)}{2}},`$ (8)
are calculated at each step in time from the solution of the stationary Klein-Gordon equations

$$\mathrm{\Delta }\varphi _m+U^{}(\varphi _m)=\pm \overline{\psi }\mathrm{\Gamma }_m\psi ,$$
(9)
where the upper sign holds for vector fields and the lower sign for the scalar field. The index $`m`$ denotes mesons and the photon, i.e. $`\varphi _m=\sigma ,\omega ^\mu ,\stackrel{}{\rho }^\mu ,A^\mu `$. This approximation is justified by the large meson masses. The corresponding meson exchange forces are of short range and therefore retardation effects can be neglected.

In practical applications to nuclear matter and finite nuclei, the relativistic models are used in the no-sea approximation: the Dirac sea of states with negative energies does not contribute to the densities and currents. For a nucleus with A nucleons

$$\overline{\psi }\mathrm{\Gamma }_m\psi =\underset{i=1}{\overset{A}{}}\overline{\psi }_i(𝐫,t)\mathrm{\Gamma }_m\psi _i(𝐫,t),$$
(10)
where the summation is performed only over the occupied orbits in the Fermi sea of positive energy states. The set of coupled equations (7) and (9) define the relativistic mean field (RMF) model. In the stationary case they reduce to a nonlinear eigenvalue problem, and in the time-dependent case they describe the nonlinear propagation of the Dirac spinors in time VBR.95 .

The mean-field approximation represents the lowest order of the quantum field theory: the meson field operators are replaced by their expectation values. The A nucleons, described by a Slater determinant $`|\mathrm{\Phi }`$ of single-particle spinors, move independently in the classical meson fields. The couplings of the meson fields to the nucleon are adjusted to reproduce the properties of nuclear matter and finite nuclei. The $`\sigma `$-meson approximates a large attractive scalar field, that is produced by very complicated microscopic processes, such as uncorrelated and correlated pion-exchange. The $`\omega `$-meson describes the short range repulsion between the nucleons, and the $`\rho `$-meson carries the isospin quantum number. The latter is required by the large empirical asymmetry potential in finite nuclear systems. The basic ingredient of the microscopic nuclear force is the pion. In relativistic mean-field models it does not contribute on the Hartree level because of parity conservation. The pion field has been included in the relativistic Hartree-Fock model. However, the resulting equations of motion are rather complicated and this model has been rarely used. Many effects that go beyond the mean-field level are apparently neglected in the RMF model. Among them are the Fock terms, the vacuum polarization effects and the short range Brueckner-type correlations. The experimental data to which the meson-nucleon couplings are adjusted, however, contain all these effects and much more. It follows that these effects are not completely neglected. On the contrary, they are taken into account in an effective way. The concept behind the RMF model is therefore equivalent to that of the density functional theory, which is widely used in solid state physics, molecular physics, chemistry and also in non-relativistic nuclear physics. The RMF model represent the covariant form of this method.

### II.2 Covariant density functional theory

The equations of motion of the relativistic mean field model can also be derived starting from a density functional. From the energy-momentum tensor one writes the total energy of the nuclear system

$`E_{RMF}[\psi ,\overline{\psi },\sigma ,\omega ^\mu ,\stackrel{}{\rho }^\mu ,A^\mu ]`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{A}{}}}{\displaystyle \psi _i^+(𝜶𝒑+\beta m)\psi _i}`$ (11)
$`+`$ $`{\displaystyle }[{\displaystyle \frac{1}{2}}(\sigma )^2+U(\sigma )]d^3r{\displaystyle \frac{1}{2}}{\displaystyle }[(\omega )^2+m_\omega ^2\omega ^2)]d^3r`$
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle }[(\rho )^2+m_\rho ^2\rho ^2)]d^3r{\displaystyle \frac{1}{2}}{\displaystyle }(A)^2d^3r`$
$`+`$ $`{\displaystyle \left[g_\sigma \rho _s\sigma +g_\omega j_\mu \omega ^\mu +g_\rho \stackrel{}{j}_\mu \stackrel{}{\rho }^\mu +ej_{c\mu }A^\mu \right]d^3r}.`$
By using the definition of the relativistic single nucleon density matrix

$$\widehat{\rho }(𝐫,𝐫^{},t)=\underset{i=1}{\overset{A}{}}|\psi _i(𝐫,t)\psi _i(𝐫^{},t)|,$$
(12)
the total energy can be written as a functional of the density matrix $`\widehat{\rho }`$ and the meson fields

$`E_{RMF}[\widehat{\rho },\varphi _m]`$ $`=`$ $`\text{Tr}\left[(𝜶𝒑+\beta m)\widehat{\rho }\right]`$ (13)
$`\pm {\displaystyle \left[\frac{1}{2}(\varphi _m)^2+U(\varphi _m)\right]d^3r}`$
$`+\text{Tr}\left[(\mathrm{\Gamma }_m\varphi _m)\widehat{\rho }\right].`$
The trace operation involves a sum over the Dirac indices and an integral in coordinate space. The index m is used as generic notation for all mesons and the photon. From the classical time-dependent variational principle

$$\delta _{t_1}^{t_2}𝑑t\left\{\mathrm{\Phi }|i_t|\mathrm{\Phi }E[\widehat{\rho },\varphi _m]\right\}=0$$
(14)
the equations of motion (7) and (9) are obtained. The equation of motion for the density matrix reads

$`i_t\widehat{\rho }`$ $`=`$ $`[\widehat{h}(\widehat{\rho },\varphi _m),\widehat{\rho }].`$ (15)
The single particle Hamiltonian $`\widehat{h}`$ is the functional derivative of the energy with respect to the single particle density matrix $`\widehat{\rho }`$

$$\widehat{h}=\frac{\delta E}{\delta \widehat{\rho }}.$$
(16)
### II.3 Pairing correlations and the Relativistic Hartree-Bogoliubov theory

The inclusion of pairing correlations is essential for a quantitative description open shell nuclei. In Ref. KR.91  a fully microscopic derivation of the relativistic Hartree-Bogoliubov theory has been developed. Using the Gorkov factorization technique, it has been shown that the pairing interaction results from the one-meson exchange ($`\sigma `$-, $`\omega `$\- and $`\rho `$-mesons). In practice, however, it turns out that the pairing correlations calculated in this way, with coupling constants taken from the standard parameter sets of the RMF model, are too strong. The repulsion produced by the exchange of vector mesons at short distances results in a pairing gap at the Fermi surface that is by a factor three too large. However, as has been argued in many applications of the Hartree-Fock-Bogoliubov theory, there is no real reason to use the same effective forces in both the particle-hole and particle-particle channels.

Pairing correlations can be easily included in the framework of density functional theory, by using a generalized Slater determinant $`|\mathrm{\Phi }`$ of the Hartree-Bogoliubov type. The ground state of a nucleus $`|\mathrm{\Phi }`$ is represented as the vacuum with respect to independent quasiparticle operators

$$\alpha _k^+=\underset{l}{}U_{lk}c_l^++V_{lk}c_l,$$
(17)
where $`U_{lk}`$, $`V_{lk}`$ are the Hartree-Bogoliubov coefficients. They determine the Hermitian single particle density matrix

$$\widehat{\rho }=V^{}V^T,$$
(18)
and the antisymmetric pairing tensor

$$\widehat{\kappa }=V^{}U^T.$$
(19)
The energy functional depends not only on the density matrix $`\widehat{\rho }`$ and the meson fields $`\varphi _m`$, but in addition also on the pairing tensor. It has the form

$$E[\widehat{\rho },\widehat{\kappa },\varphi _m]=E_{RMF}[\widehat{\rho },\varphi _m]+E_{pair}[\widehat{\kappa }],$$
(20)
where $`E_{RMF}[\widehat{\rho },\varphi ]`$ is the $`RMF`$-functional defined in Eq. (13). The pairing energy $`E_{pair}[\widehat{\kappa }]`$ is given by

$$E_{pair}[\widehat{\kappa }]=\frac{1}{4}\text{Tr}\left[\widehat{\kappa }^{}V^{pp}\widehat{\kappa }\right].$$
(21)
$`V^{pp}`$ is a general two-body pairing interaction. Finally, the total energy can be written as a functional of the generalized density matrix Va.61

$$=\left(\begin{array}{cc}\rho & \kappa \\ \kappa ^{}& 1\rho ^{}\end{array}\right),$$
(22)
which obeys the equation of motion

$$i_t=[(),].$$
(23)
The generalized Hamiltonian $``$ is a functional derivative of the energy with respect to the generalized density

$$=\frac{\delta E}{\delta }=\left(\begin{array}{cc}\widehat{h}_Dm\lambda & \widehat{\mathrm{\Delta }}\\ \widehat{\mathrm{\Delta }}^{}& \widehat{h}_D+m+\lambda \end{array}\right).$$
(24)
It contains two average potentials: the self-consistent mean field $`\widehat{h}_D`$, which encloses all the long range particle-hole (ph) correlations, and the pairing field $`\widehat{\mathrm{\Delta }}`$, which includes the particle-particle (pp) correlations. The single particle potential $`\widehat{h}_D`$ results from the variation of the energy functional with respect to the Hermitian density matrix $`\widehat{\rho }`$

$$\widehat{h}_D=\frac{\delta E}{\delta \widehat{\rho }},$$
(25)
and the pairing field is obtained from the variation of the energy functional with respect to the pairing tensor

$$\widehat{\mathrm{\Delta }}=\frac{\delta E}{\delta \widehat{\kappa }}.$$
(26)
The pairing field is an integral operator with the kernel

$$\mathrm{\Delta }_{ab}(𝐫,𝐫^{})=\frac{1}{2}\underset{c,d}{}V_{abcd}^{pp}(𝐫,𝐫^{})\kappa _{cd}(𝐫,𝐫^{}),$$
(27)
where $`a,b,c,d`$ denote quantum numbers that specify the Dirac indices of the spinors, and $`V_{abcd}^{pp}(𝐫,𝐫^{})`$ are the matrix elements of a general two-body pairing interaction.

The stationary limit of Eq. (23) describes the ground state of an open-shell nucleus Gon.96 ; Lal1.98 . It is determined by the solutions of the Hartree-Bogoliubov equations

$$\left(\begin{array}{cc}\widehat{h}_Dm\lambda & \widehat{\mathrm{\Delta }}\\ \widehat{\mathrm{\Delta }}^{}& \widehat{h}_D+m+\lambda \end{array}\right)\left(\begin{array}{c}U_k(𝐫)\\ V_k(𝐫)\end{array}\right)=E_k\left(\begin{array}{c}U_k(𝐫)\\ V_k(𝐫)\end{array}\right).$$
(28)
The chemical potential $`\lambda `$ is determined by the particle number subsidiary condition in order that the expectation value of the particle number operator in the ground state equals the number of nucleons. The column vectors denote the quasiparticle wave functions, and $`E_k`$ are the quasiparticle energies. The dimension of the RHB matrix equation is two times the dimension of the corresponding Dirac equation. For each eigenvector $`(U_k,V_k)`$ with positive quasiparticle energy $`E_k>0`$, there exists an eigenvector $`(V_k^{},U_k^{})`$ with quasiparticle energy $`E_k`$. Since the baryon quasiparticle operators satisfy fermion commutation relations, the levels $`E_k`$ and $`E_k`$ cannot be occupied simultaneously. For the solution that corresponds to a ground state of a nucleus with even particle number, one usually chooses the eigenvectors with positive eigenvalues $`E_k`$.

The RHB equations are solved self-consistently, with potentials determined in the mean-field approximation from solutions of static Klein-Gordon equations

$`\left[\mathrm{\Delta }+m_\sigma ^2\right]\sigma (𝐫)`$ $`=`$ $`g_\sigma \rho _s(𝐫)g_2\sigma ^2(𝐫)g_3\sigma ^3(𝐫)`$ (29)
$`\left[\mathrm{\Delta }+m_\omega ^2\right]\omega ^0(𝐫)`$ $`=`$ $`g_\omega \rho _v(𝐫)`$ (30)
$`\left[\mathrm{\Delta }+m_\rho ^2\right]\rho ^0(𝐫)`$ $`=`$ $`g_\rho \rho _3(𝐫)`$ (31)
$`\mathrm{\Delta }A^0(𝐫)`$ $`=`$ $`e\rho _p(𝐫)`$ (32)
for the $`\sigma `$-meson, the $`\omega `$-meson, the $`\stackrel{}{\rho }`$-meson and the photon field, respectively. Because of charge conservation, only the 3-rd component of the isovector $`\rho `$-meson contributes. In the ground-state solution for an even-even nucleus there are no currents (time reversal invariance) and the spatial components $`𝝎`$, $`𝝆_3`$, $`𝑨`$ of the vector fields vanish. In nuclei with an odd number of protons or neutrons time reversal symmetry is broken, and the resulting spatial components of the meson fields play an essential role in the description of magnetic moments and of moments of inertia in rotating nuclei. The equation for the isoscalar scalar $`\sigma `$-meson field contains nonlinear terms. The inclusion of nonlinear meson self-interaction terms in meson-exchange RMF models is absolutely necessary for a quantitative description of ground-state properties of spherical and deformed nuclei Rin.96 . The source terms in equations (29) to (32) are sums of bilinear products of baryon amplitudes

$`\rho _s(𝐫)`$ $`=`$ $`{\displaystyle \underset{k>0}{}}V_k^{}(𝐫)\gamma ^0V_k(𝐫)`$ (33)
$`\rho _v(𝐫)`$ $`=`$ $`{\displaystyle \underset{k>0}{}}V_k^{}(𝐫)V_k(𝐫)`$ (34)
$`\rho _3(𝐫)`$ $`=`$ $`{\displaystyle \underset{k>0}{}}V_k^{}(𝐫)\tau _3V_k(𝐫)`$ (35)
$`\rho _{\mathrm{em}}(𝐫)`$ $`=`$ $`{\displaystyle \underset{k>0}{}}V_k^{}(𝐫){\displaystyle \frac{1\tau _3}{2}}V_k(𝐫),`$ (36)
where $`_{k>0}`$ is a shorthand notation for the no-sea approximation. The self-consistent solution of the Dirac-Hartree-Bogoliubov integro-differential equations and Klein-Gordon equations for the meson fields determines the ground state of a nucleus. In the present implementation of the RHB model the coupled system of equations is solved by expanding the nucleon spinors $`U_k(𝐫)`$ and $`V_k(𝐫)`$, and the meson fields in the spherical harmonic oscillator basis GRT.90 .

### II.4 The Relativistic Quasi-Particle Random Phase Approximation

In this section we will derive the RQRPA from the time-dependent RHB model in the limit of small amplitude oscillations. The generalized density matrix $``$ and the fields $`\varphi _m=\sigma ,\omega ^\mu ,\stackrel{}{\rho }^\mu ,A^\mu `$ have been considered as independent variables related only by the equations of motion. One can use the Klein-Gordon equations to eliminate the meson degrees of freedom, but this is only possible in the small amplitude limit. The time-dependent meson field can be written as

$$\varphi _m=\varphi _m^{(0)}+\delta \varphi _m,$$
(37)
where $`\varphi _m^{(0)}`$ is the meson field that corresponds to the stationary ground state, and $`\delta \varphi _m`$ is a small variation of the meson field around the stationary state solution. In the linear approximation the corresponding Klein-Gordon equation reads

$$\left[\mathrm{\Delta }+U^{\prime \prime }(\varphi _m^{(0)})\right]\delta \varphi _m(𝐫)=\pm g_m\delta \rho _m(𝐫),$$
(38)
where $`\delta \rho _m(𝐫)`$ are the various densities and currents (see Eq. (10)). If there are no nonlinear meson self-interaction terms, $`U^{\prime \prime }(\varphi _m^{(0)})=m_m^2`$. The propagator $`G_m(𝐫,𝐫^{})`$ can be obtained analytically and it has the Yukawa form. In the case of nonlinear meson self-interaction terms $`U^{\prime \prime }(\varphi _m^{(0)})`$ depends on the field $`\varphi _m^{(0)}`$, and an analytical solution is no longer possible. The propagator $`G_m(𝐫,𝐫^{})`$ has to be calculated numerically (for details see Ref. MGT.97 ). In both cases we find a linear relation between $`\delta \varphi _m`$ and $`\delta \rho _m`$

$$\delta \varphi _m(𝐫)=\pm g_md^3r^{}G_m(𝐫,𝐫^{})\delta \rho _m(𝐫^{}).$$
(39)
The generalized Hamiltonian $``$ can now be expressed as a functional of the generalized density $``$ only. In the linear approximation the generalized density matrix is expanded

$$=_0+\delta (t),$$
(40)
where $`_0`$ is the stationary ground-state generalized density. Since $`(t)`$ is a projector at all times, in linear order

$$_0\delta +\delta _0=\delta .$$
(41)
In the quasiparticle basis the matrices $`_0`$ and $`_0=(_0)`$ are diagonal

$$_0=\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)\text{and }_0=\left(\begin{array}{cc}E_n& 0\\ 0& E_n\end{array}\right).$$
(42)
From Eq. (41) it follows that the matrix $`\delta `$ has the form

$$\delta =\left(\begin{array}{cc}0& \delta R\\ \delta R^{}& 0\end{array}\right).$$
(43)
The linearized equation of motion (23) reduces to

$$i_t=[_0,\delta ]+[\frac{\delta }{\delta }\delta ,_0].$$
(44)
Assuming an oscillatory solution

$$\delta (t)=\underset{\nu }{}\delta ^{(\nu )}e^{i\omega _\nu t}+h.c.,$$
(45)
the RQRPA equation is obtained

$$\left(\begin{array}{cc}A& B\\ B^{}& A^{}\end{array}\right)\left(\begin{array}{c}X^\nu \\ Y^\nu \end{array}\right)=\omega _\nu \left(\begin{array}{c}X^\nu \\ Y^\nu \end{array}\right).$$
(46)
For $`k<k^{}`$, $`l<l^{}`$ the RQRPA matrix elements read

$$A_{kk^{}ll^{}}=(E_k+E_k^{})\delta _{kl}\delta _{k^{}l^{}}+\frac{\delta ^2E}{\delta R_{kk^{}}^{}\delta R_{ll^{}}}\text{ and }B_{kk^{}ll^{}}=\frac{\delta ^2E}{\delta R_{kk^{}}^{}\delta R_{ll^{}}^{}}.$$
(47)
If the two-body Hamiltonian is density independent the matrices $`A`$ and $`B`$ have the simple forms Rin.80

$`A_{kk^{},ll^{}}`$ $`=`$ $`\mathrm{\Phi }\left|[\alpha _k^{}\alpha _k,[\widehat{H},\alpha _l^+\alpha _l^{}^+]]\right|\mathrm{\Phi }`$
$`B_{kk^{},ll^{}}`$ $`=`$ $`\mathrm{\Phi }\left|[\alpha _k^{}\alpha _k,[\widehat{H},\alpha _l^{}\alpha _l]]\right|\mathrm{\Phi }.`$ (48)
Using the representation of the Hamiltonian in the quasiparticle basis

$`\widehat{H}`$ $`=`$ $`E_0+{\displaystyle \underset{kk^{}}{}}H_{kk^{}}^{11}\alpha _k^+\alpha _k^{}+{\displaystyle \frac{1}{4}}{\displaystyle \underset{kk^{}ll^{}}{}}H_{kk^{}ll^{}}^{22}\alpha _k^+\alpha _k^{}^+\alpha _l^{}\alpha _l`$ (49)
$`+`$ $`{\displaystyle \underset{kk^{}ll^{}}{}}(H_{kk^{}ll^{}}^{40}\alpha _k^+\alpha _k^{}^+\alpha _l^{}^+\alpha _l^++h.c.)+{\displaystyle \underset{kk^{}ll^{}}{}}(H_{kk^{}ll^{}}^{31}\alpha _k^+\alpha _k^{}^+\alpha _l^{}^+\alpha _l+h.c.)`$
we find

$`A_{kk^{}ll^{}}`$ $`=`$ $`H_{kl}^{11}\delta _{k^{}l^{}}H_{k^{}l}^{11}\delta _{kl^{}}H_{kl^{}}^{11}\delta _{k^{}l}+H_{k^{}l^{}}^{11}\delta _{kl}+H_{kk^{}ll^{}}^{22}`$
$`B_{kk^{}ll^{}}`$ $`=`$ $`4H_{kk^{}ll^{}}^{40}.`$ (50)
In the quasiparticle representation the matrix $`H^{11}`$ is diagonal, i.e. $`H_{kl}^{11}=E_k\delta _{kl}`$. The matrices $`H^{22}`$ and $`H^{40}`$ are rather complicated expressions containing the two-body $`ph`$\- and $`pp`$-matrix elements and the coefficients $`U`$ and $`V`$ (for details see Ref. Rin.80 )
In the more general case of a density dependent Hamiltonian the same expressions can be used, but one has to take into account the rearrangement terms originating from the variation of the interaction with respect to the density $`\widehat{\rho }`$.

### II.5 The Relativistic QRPA in the canonical basis

The full RQRPA equations are rather complicated, because they require the evaluation of the matrix elements $`H_{kk^{}ll^{}}^{22}`$ and $`H_{kk^{}ll^{}}^{40}`$ in the basis of the Hartree-Bogoliubov spinors $`U_k(𝐫)`$ and $`V_k(𝐫)`$. It is considerably simpler to solve these equations in the canonical basis, in which the relativistic Hartree-Bogoliubov wave functions can be expressed in the form of BCS-like wave functions. In this case one needs only the matrix elements $`V_{\kappa \lambda ^{}\kappa ^{}\lambda }^{ph}`$ of the residual $`ph`$-interaction, and the matrix elements $`V_{\kappa \kappa ^{}\lambda \lambda ^{}}^{pp}`$ of the pairing $`pp`$-interaction, as well as certain combinations of the occupation factors $`u_\kappa `$, $`v_\kappa `$. The numerical details are described in the Appendix. In the following we use the indices $`\kappa `$, $`\lambda `$, $`\kappa ^{}`$ and $`\lambda ^{}`$ to denote states in the canonical basis. We emphasize that the solution of the relativistic quasi-particle RPA equations in the canonical basis does not represent an approximation. We obtain a full solution and the results do not depend on this special choice of the basis.

Taking into account the rotational invariance of the nuclear system, the quasiparticle pairs can be coupled to good angular momentum and the matrix equations of the relativistic quasiparticle random phase approximation (RQRPA) read

$$\left(\begin{array}{cc}A^J& B^J\\ B^{}_{}{}^{}J& A^{}_{}{}^{}J\end{array}\right)\left(\begin{array}{c}X^{\nu ,JM}\\ Y^{\nu ,JM}\end{array}\right)=\omega _\nu \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}X^{\nu ,JM}\\ Y^{\nu ,JM}\end{array}\right).$$
(51)
For each RQRPA energy $`\omega _\nu `$, $`X^\nu `$ and $`Y^\nu `$ denote the corresponding forward- and backward-going two-quasiparticle amplitudes, respectively. The coupled RQRPA matrices in the canonical basis read

$`A_{\kappa \kappa ^{}\lambda \lambda ^{}}^J`$ $`=`$ $`H_{\kappa \lambda }^{11(J)}\delta _{\kappa ^{}\lambda ^{}}H_{\kappa ^{}\lambda }^{11(J)}\delta _{\kappa \lambda ^{}}H_{\kappa \lambda ^{}}^{11(J)}\delta _{\kappa ^{}\lambda }+H_{\kappa ^{}\lambda ^{}}^{11(J)}\delta _{\kappa \lambda }`$ (52)
$`+{\displaystyle \frac{1}{2}}(\xi _{\kappa \kappa ^{}}^+\xi _{\lambda \lambda ^{}}^++\xi _{\kappa \kappa ^{}}^{}\xi _{\lambda \lambda ^{}}^{})V_{\kappa \kappa ^{}\lambda \lambda ^{}}^{ppJ}`$
$`+\zeta _{\kappa \kappa ^{}\lambda \lambda ^{}}V_{\kappa \lambda ^{}\kappa ^{}\lambda }^{phJ}`$
$`B_{\kappa \kappa ^{}\lambda \lambda ^{}}^J`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\xi _{\kappa \kappa ^{}}^+\xi _{\lambda \lambda ^{}}^+\xi _{\kappa \kappa ^{}}^{}\xi _{\lambda \lambda ^{}}^{})V_{\kappa \kappa ^{}\lambda \lambda ^{}}^{ppJ}`$ (53)
$`+\zeta _{\kappa \kappa ^{}\lambda \lambda ^{}}(1)^{j_\lambda j_\lambda ^{}+J}V_{\kappa \lambda \kappa ^{}\lambda ^{}}^{phJ}.`$
$`H^{11}`$ denotes the one-quasiparticle terms

$$H_{\kappa \lambda }^{11}=(u_\kappa u_\lambda v_\kappa v_\lambda )h_{\kappa \lambda }(u_\kappa v_\lambda +v_\kappa u_\lambda )\mathrm{\Delta }_{\kappa \lambda },$$
(54)
i.e. the canonical RHB basis does not diagonalize the Dirac single-nucleon mean-field Hamiltonian $`\widehat{h}_D`$ and the pairing field $`\widehat{\mathrm{\Delta }}`$. The occupation amplitudes $`v_k`$ of the canonical states are eigenvalues of the density matrix. $`V^{ph}`$ and $`V^{pp}`$ are the particle-hole and particle-particle residual interactions, respectively. Their matrix elements are multiplied by the pairing factors $`\xi ^\pm `$ and $`\zeta `$, defined below by the occupation amplitudes of the canonical states. The relativistic particle-hole interaction $`V^{ph}`$ is defined by the same effective Lagrangian density as the mean-field Dirac single-nucleon Hamiltonian $`\widehat{h}_D`$. $`V^{ph}`$ includes the exchange of the isoscalar scalar $`\sigma `$-meson, the isoscalar vector $`\omega `$-meson, the isovector vector $`\rho `$-meson, and the electromagnetic interaction. The two-body matrix elements include contributions from the spatial components of the vector fields.

$$\zeta _{\kappa \kappa ^{}\lambda \lambda ^{}}=\{\begin{array}{cc}\eta _{\kappa \kappa ^{}}^+\eta _{\lambda \lambda ^{}}^+\hfill & \text{for }\sigma \text{, and the time-components }\omega ^0\text{}\rho ^0\text{}A^0\text{ if J is even}\hfill \\ & \text{for the space-components }𝝎\text{}𝝆\text{}𝑨\text{ if J is odd}\hfill \\ \eta _{\kappa \kappa ^{}}^{}\eta _{\lambda \lambda ^{}}^{}\hfill & \text{for }\sigma \text{, and the time-components }\omega ^0\text{}\rho ^0\text{}A^0\text{ if J is odd}\hfill \\ & \text{for the space components }𝝎\text{}𝝆\text{}𝑨\text{ if J is even}\hfill \end{array}$$

with the $`\eta `$-coefficients defined by

$$\eta _{\kappa \kappa ^{}}^\pm =u_\kappa v_\kappa ^{}\pm v_\kappa u_\kappa ^{},$$

and

$$\xi _{\kappa \kappa ^{}}^\pm =u_\kappa u_\kappa ^{}v_\kappa v_\kappa ^{}.$$

The RQRPA configuration space includes the Dirac sea of negative energy states. In addition to the configurations built from two-quasiparticle states of positive energy, the RQRPA configuration space must also contain pair-configurations formed from the fully or partially occupied states of positive energy and the empty negative-energy states from the Dirac sea. The inclusion of configurations built from occupied positive-energy states and empty negative-energy states is essential for current conservation and the decoupling of spurious states Daw.90 . In recent applications of the relativistic RPA it has been shown that the fully consistent inclusion of the Dirac sea of negative energy states in the RRPA configuration space is essential for a quantitative comparison with the experimental excitation energies of giant resonances Vre.00 ; Ring.01 .

It should be emphasized that the present RQRPA model is fully consistent: the same interactions, both in the particle-hole and particle-particle channels, are used in the RHB equation (28) that determines the canonical quasiparticle basis, and in the RQRPA equation (51). In both channels the same strength parameters of the interactions are used in the RHB and RQRPA calculations. No additional adjustment of the parameters is needed in RQRPA calculations. This is an essential feature of our calculations and it ensures that RQRPA amplitudes do not contain spurious components associated with the mixing of the nucleon number in the RHB ground state (for $`0^+`$ excitations), or with the center-of-mass translational motion (for $`1^{}`$ excitations).

In the next section we present results of illustrative RQRPA calculations of the multipole response in spherical nuclei. For the multipole operator $`\widehat{Q}_{\lambda \mu }`$ the response function $`R(E)`$ is defined

$$R(E,J)=\underset{\nu }{}B(J,\omega _\nu )\frac{1}{\pi }\frac{\mathrm{\Gamma }/2}{(E\omega _\nu )^2+(\mathrm{\Gamma }/2)^2},$$
(55)
where $`\mathrm{\Gamma }`$ is the width of the Lorentzian distribution, and

$`B(J,\omega _\nu )`$ $`=`$ $`|{\displaystyle \underset{\kappa \kappa ^{}}{}}\{X_{\kappa \kappa ^{}}^{\nu ,J0}\kappa ||\widehat{Q}_J||\kappa ^{}`$ (56)
$`+`$ $`(1)^{j_\kappa j_\kappa ^{}+J}Y_{\kappa \kappa ^{}}^{\nu ,J0}\kappa ^{}||\widehat{Q}_J||\kappa \}(u_\kappa v_\kappa ^{}+(1)^Jv_\kappa u_\kappa ^{})|^2.`$
In all the examples considered in Sec. III, the discrete strength distribution are folded by a Lorentzian of width $`\mathrm{\Gamma }`$=1MeV . For the state $`|J,\nu `$, the RQRPA transition density reads

$`\delta \rho _J^\nu (r)`$ $`=`$ $`{\displaystyle \underset{\kappa \kappa ^{}}{}}\left\{\kappa Y_J\kappa ^{}f_\kappa (r)f_\kappa ^{}(r)+\widehat{\kappa }Y_J\widehat{\kappa }^{}g_\kappa (r)g_\kappa ^{}(r)\right\}`$ (57)
$`\left(X_{\kappa \kappa ^{}}^{\nu ,J0}+(1)^JY_{\kappa \kappa ^{}}^{\nu ,J0}\right)(u_\kappa v_\kappa ^{}+(1)^Jv_\kappa u_\kappa ^{}),`$
where $`\kappa `$ and $`\widehat{\kappa }`$ denote the quantum numbers of the large and small components of the Dirac spinors, respectively. $`f_\kappa (r)`$ and $`g_\kappa (r)`$ are the corresponding large and small radial components.

## III Illustrative calculations and tests of the RQRPA

Nuclear properties calculated with the RHB+RQRPA model will, of course, crucially depend on the choice of the effective RMF Lagrangian in the $`ph`$-channel, as well as on the treatment of pairing correlations. The most successful RMF effective interactions are purely phenomenological, with parameters adjusted to reproduce the nuclear matter equation of state and a set of global properties of spherical closed-shell nuclei. In most applications of the RHB model, in particular, we have used the NL3 effective interaction Lal.97  for the RMF effective Lagrangian. Properties calculated with NL3 indicate that this is probably the best nonlinear effective interaction so far, both for nuclei at and away from the line of $`\beta `$-stability. In the $`pp`$-channel of the RHB model we have used a phenomenological pairing interaction, the pairing part of the Gogny force,

$$V^{pp}(1,2)=\underset{i=1,2}{}e^{((𝐫_1𝐫_2)/\mu _i)^2}(W_i+B_iP^\sigma H_iP^\tau M_iP^\sigma P^\tau ),$$
(58)
with the set D1S Ber.84  for the parameters $`\mu _i`$, $`W_i`$, $`B_i`$, $`H_i`$ and $`M_i`$ $`(i=1,2)`$. This force has been very carefully adjusted to the pairing properties of finite nuclei all over the periodic table. In particular, the basic advantage of the Gogny force is the finite range, which automatically guarantees a proper cut-off in momentum space. All RHB+RQRPA calculations presented in this work have been performed with the NL3+D1S combination of effective interactions.

In order to illustrate the RHB+RQRPA approach and to test the numerical implementation of the RQRPA equations, in this section we calculate the isoscalar monopole, isovector dipole and isoscalar quadrupole response of 22O. Similar calculations for the neutron-rich oxygen isotopes were recently performed by Matsuo Mat.01 ; Mat.02  in the framework of the non-relativistic continuum linear response theory based on the Hartree-Fock-Bogoliubov formalism in coordinate state representation. The two theoretical frameworks differ, of course, both in the physical contents, as well as in the numerical implementation. The results can, nevertheless, be compared at least at the qualitative level. In the HFB+QRPA model of Refs. Mat.01 ; Mat.02  a Woods-Saxon parameterization is adopted for the single-particle potential, and a Skyrme-type density dependent delta force is used for the residual interaction in the $`ph`$-channel of the QRPA. Since the calculation of the single-particle potential and $`ph`$-interaction is not self-consistent, the interaction strength of the residual interaction is renormalized for each nucleus in such a way that the dipole response has a zero-energy mode corresponding to the spurious center of mass motion. For the pairing interaction, a density-dependent delta force is used both in the calculation of the HFB pairing field for the ground state, and in the linear response equation for the excitations. The calculation is consistent in the $`pp`$-channel. The present RHB+RQRPA calculations are fully self-consistent: the same combination of effective interactions, NL3 in the $`ph`$-channel and Gogny D1S in the $`pp`$-channel, are used both in the RHB calculation of the ground state and as RQRPA residual interactions. The parameters of the RQRPA residual interactions have exactly the same values as those used in the RHB calculation.

In the analysis of Refs. Mat.01 ; Mat.02 , Matsuo has illustrated the importance of a consistent treatment of pairing correlations in the HFB+QRPA framework. The residual pairing interaction in the QRPA generates pronounced dynamical correlation effects on the responses through pair density fluctuations. Moreover, the energy weighted sum rules are only satisfied if the pairing interaction is consistently included both in the static HFB and in the dynamical linear response. We have verified that the results obtained in the HFB+QRPA framework are also reproduced in the RHB+RQRPA calculations.

In the left panel of Fig. 1 we display the monopole strength function of the neutron number operator in 22O. There should be no response to the number operator since it is a conserved quantity, i.e. the Nambu-Goldstone mode associated with the nucleon number conservation should have zero excitation energy. The dashed curve (no dynamical pairing) represents the strength function obtained when the pairing interaction is included only in the RHB calculation of the ground state, but not in the residual interaction of the RQRPA. The solid line (zero response) corresponds to the full RHB+RQRPA calculation, with the pairing interaction included both in the RHB ground state, and in the RQRPA residual interaction. The same result was also obtained in the HFB+QRPA calculation for 24O in Ref. Mat.02 : the spurious strength of the number operator appears when the pairing interaction is included only in the stationary solution for the ground state, i.e. when the dynamical QRPA pairing correlations are neglected.

The isoscalar strength functions of the monopole operator $`_{i=1}^Ar_i^2`$ in 22O, shown in the right panel of Fig. 1, correspond to three different calculations: a) the RMF+RRPA calculation without pairing, b) pairing correlations are included in the RHB calculation of the ground state, but not in the RQRPA residual interaction (no dynamical pairing), and c) the fully self-consistent RHB+RQRPA calculation. Just as in the case of the number operator, by including pairing correlations only in the RHB ground state a strong spurious response is generated below 10 MeV. The Nambu-Goldstone mode is found at zero excitation energy (in this particular calculation it was located below 0.2 MeV) only when pairing correlations are consistently included also in the residual RQRPA interaction. When the result of the full RHB+RQRPA is compared with the response calculated without pairing, one notices that, as expected, pairing correlations have relatively little influence on the response in the region of giant resonances above 20 MeV. A more pronounced effect is found at lower energies. The fragmentation of the single peak at $``$ 12.5 MeV reflects the broadening of the Fermi surface by the pairing correlations.

The isovector strength function ($`J^\pi =1^{}`$) of the dipole operator

$$\widehat{Q}_{1m}^{T=1}=\frac{N}{N+Z}\underset{p=1}{\overset{Z}{}}r_pY_{1m}\frac{Z}{N+Z}\underset{n=1}{\overset{N}{}}r_nY_{1m}$$
(59)
for 22O is displayed in the left panel of Fig. 2. In this example we also compare the results of the RMF+RRPA calculations without pairing, with pairing correlations included only in the RHB ground state (no dynamical pairing), and with the fully self-consistent RHB+RQRPA response. A large configuration space enables the separation of the zero-energy mode that corresponds to the spurious center of mass motion. In the present calculation for 22O this mode is found at $`E=0.04`$ MeV.

The isovector dipole response in neutron-rich oxygen isotopes has recently attracted considerable interest because these nuclei might be good candidates for a possible identification of the low-lying collective soft mode (pygmy state), that corresponds to the oscillations of excess neutrons out of phase with the core composed of an equal number of protons and neutrons Try.01 ; Lei.01 . The strength functions shown in Fig. 2 illustrate the importance of including pairing correlations in the calculation of the isovector dipole response. Pairing is, of course, particularly important for the low-lying strength below 10 MeV. The inclusion of pairing correlations in the full RHB+RQRPA calculation enhances the low-energy dipole strength near the threshold. For the main peak in the low-energy region ($``$ 8.65 MeV), in the right panel of Fig. 2 we display the proton and neutron transition densities. In contrast to the well known radial dependence of the IVGDR transition densities (proton and neutron densities oscillate with opposite phases, the amplitude of the isovector transition density is much larger than that of the isoscalar component), the proton and neutron transition densities for the main low-energy peak are in phase in the nuclear interior, there is no contribution from the protons in the surface region, the isoscalar transition density dominates over the isovector one in the interior, and the strong neutron transition density displays a long tail in the radial coordinate. A similar behavior has been predicted for the light neutron halo nuclei 6He, 11Li and 12Be in Ref. Sag.96 , where it has been shown that the long tails of the wave functions of the loosely-bound neutrons are responsible for the different radial dependence of the transition densities that correspond to the soft low-energy states as compared to those of the giant resonances.

The effect of pairing correlations on the isovector dipole response in 22O is very similar to the one obtained in the HFB+QRPA framework (Fig. 8 of Ref. Mat.02 ). In the low-energy region below 10 MeV, however, the pairing interaction used in the QRPA calculation produces a much stronger enhancement of the dipole strength, as compared to the results shown in Fig. 2. The reason probably lies in the choice of the pairing interaction. While we use the volume Gogny pairing, in Ref. Mat.02  a density-dependent delta force was used in the $`pp`$ channel. This interaction is surface peaked and therefore produces a stronger effect on the low-energy dipole strength near the threshold. Nevertheless, we emphasize that the RHB+RQRPA results for the low-lying dipole strength distribution in 22O are in very good agreement with recent experimental data Lei.01 .

In the left panel of Fig. 3 we display the RHB+RQRPA isoscalar and isovector quadrupole ($`J^\pi =2^+`$) strength distributions in 22O. The low-lying $`J^\pi =2^+`$ state is calculated at $`E=2.95`$ MeV, and this value should be compared with the experimental excitation energy of the first $`2^+`$ state: 3.2 MeV Bell.01 . The strong peak at $`E=22.3`$ MeV in the isoscalar strength function corresponds to the isoscalar giant quadrupole resonance (IS GQR). The isovector response, on the other hand, is strongly fragmented over the large region of excitation energies $`E1838`$ MeV. The effect of pairing correlations on the isoscalar response is illustrated in the right panel of Fig. 3, where again the full RHB+RQRPA strength function is compared to the RMF+RRPA calculation without pairing, and with the response obtained when the pairing interaction is included only in the RHB ground state (no dynamical pairing). As one would expect, the effect of pairing correlations is not particularly pronounced in the giant resonance region. The inclusion of pairing correlations, however, has a relatively strong effect on the low-lying $`2^+`$ state. This is seen more clearly in the left panel of Fig. 4, where only the low-energy portion of the isoscalar strength distributions in 22O is shown. With respect to the RRPA calculation, the inclusion of the pairing interaction in the static solution for the ground state increases the excitation energy of the lowest $`2^+`$ state by $``$ 3 MeV. The fully self-consistent RHB+RQRPA calculation lowers the excitation energy from $``$ 4.5 MeV to $`E=2.95`$ MeV. The inclusion of pairing correlations increases the collectivity of the low-lying $`2^+`$ state. A very similar result for the low-lying quadrupole state in 24O has been obtained by Matsuo in the HFB+QRPA framework Mat.02 . The proton and neutron transition densities for the $`2^+`$ state at $`E=2.95`$ MeV are shown in the right panel of Fig. 4. They display a characteristic radial dependence. Both transition densities are, of course, peaked in the surface region, but the proton contribution is much smaller. The RHB+RQRPA results for the $`2^+`$ excitations are in agreement with non-relativistic QRPA calculations of the quadrupole response in neutron rich oxygen isotopes Kha.00 ; Kha.00e ; Mat.02 ; Kha.02 .

## IV Evolution of the low-lying isovector dipole strength in Sn isotopes and N=82 isotones

The dipole response of very neutron-rich isotopes is characterized by the fragmentation of the strength distribution and its spreading into the low-energy region, and by the mixing of isoscalar and isovector modes. It appears that in most relatively light nuclei the onset of dipole strength in the low-energy region is due to non-resonant independent single particle excitations of the loosely bound neutrons. The structure of the low-lying dipole strength changes with mass. As we have shown in the RRPA analysis of Ref. Vrepyg2.01 , in heavier nuclei low-lying dipole states appear that are characterized by a more distributed structure of the RRPA amplitude. Among several peaks characterized by single particle transitions, a single collective dipole state is identified below 10 MeV, and its amplitude represents a coherent superposition of many neutron particle-hole configurations.

Very recently experimental data have been reported on the concentration of electric dipole strength below the neutron separation energy in $`N=82`$ semi-magic nuclei. The distribution of the electric dipole strength in 138Ba, 140Ce, and 144Sm displays a resonant structure between 5.5 MeV and 8 MeV, exhausting $``$ 1% of the isovector E1 EWSR Zil.02 . In 138Ba negative parity quantum numbers have been assigned to 18 dipole excitations between 5.5 MeV and 6.5 MeV Pie.02 .

In Figs. 5 and 6 we display the isovector dipole strength distributions in eight $`N=82`$ isotones, calculated in the RHB+RQRPA framework with the NL3+D1S combination of effective interactions. The calculation is fully self-consistent, with the Gogny finite-range pairing included both in the RHB ground state, and in the RQRPA residual interaction. The isovector dipole response is shown for even-Z nuclei from 146Gd to the doubly magic 132Sn. In addition to the characteristic peak of the isovector giant dipole resonance (IVGDR) at $``$ 15 MeV, the evolution of the low-lying dipole strength with decreasing proton number is clearly observed below 10 MeV. The strength of the low-lying dipole response increases with the relative increase of the neutron contribution, i.e. with reducing the number of protons. For the main peaks in the low-energy region below 10 MeV, in the panels on the right side of Figs. 5 and 6 we display the corresponding neutron and proton transition densities. The radial dependence is very different from that of the transition densities of the IVGDR peak. For all eight nuclei the main peak below 10 MeV does not correspond to an isovector excitation, i.e. the proton and neutron transition densities have the same sign. The relative contribution of the protons in the surface region decreases with reducing the proton number. In particular, for the nuclei shown in Fig. 6: 138Ba, 136Xe, 134Te and 132Sn, there is practically no proton contribution to the transition density beyond 6 fm. The dynamics is that of a pygmy resonance: the neutron skin oscillates against the core. In Ref. Zil.02  it was emphasized that the observed low-lying dipole states in the N=82 isotones are not just statistical E1 excitations sitting on the tail of the GDR, but represent a fundamental structure effect. In Fig. 7 we show that this is also the case for the RHB+RQRPA results. For the dipole strength distribution of 140Ce, shown in the left panel, in the right column we compare the neutron and proton transition densities for the IVGDR peak at 14.31 MeV, for the peak at 12.51 MeV, and for the main peak in the low-energy region at 8.22 MeV. The peak at 12.51 MeV, as well as other peaks in the interval 10-14 MeV, displays transition densities very similar to those of the GDR peak, i.e. these states belong to the tail of the GDR. The dynamics of the low-energy mode at 8.22 MeV, on the other hand, is very different: the proton and neutron transition densities are in phase in the nuclear interior, there is almost no contribution from the protons in the surface region, the isoscalar transition density dominates over the isovector one, and the peak of the strong neutron transition density in the surface region is shifted toward larger radii.

On a quantitative level, the present RHB+RQRPA calculation does not compare too well with the experimental data on the low-lying dipole strength in the N=82 isotones. First, while the observed low-energy dipole states in 138Ba, 140Ce, and 144Sm are concentrated between 5.5 MeV and 8 MeV, the calculated pygmy states in these nuclei are above 8 MeV. This can be partly explained by the low effective nucleon mass of the NL3 mean-field interaction VNR.02 . On the other hand, the excitation energies of the IV GDR are, as will be shown below in the example of Sn isotopes, rather well reproduced by the NL3 interaction. The fact that NL3 reproduces the IV GDR, but not the centroid of the low-energy dipole strength, might indicate that the isovector channel of this force needs a better parameterization. Second and more important, the number of RQRPA peaks below 10 MeV, for the operator (59), is much smaller than the number of observed dipole states in the low-energy region Pie.02 ; Zil.02 . The observed low-lying E1 strength consists of many states of different origin. This has been discussed in Ref. Zil.02 . In addition to the two-phonon and three-phonon states, and the soft pygmy state, in this energy region one could also expect some compressional low-lying isoscalar dipole strength CLY.01 , maybe mixed with toroidal states Vre.02 ; Rye.02 , as well as the E1 strength generated by the breaking of the isospin symmetry due to a clustering mechanism Iac.85 . A detailed investigation of the nature of all observed low-lying dipole states in N=82 nuclei is, of course, beyond the scope of the present analysis, since our model space does not include multi-phonon configurations.

The Sn isotopes present another very interesting example of the evolution of the low-lying dipole strength with neutron number Gov.98 . In Ref. Vrepyg2.01  we have performed an analysis of the isovector dipole response of neutron-rich Sn isotopes in the relativistic RPA framework. The RMF+RRPA calculation has shown that, among several dipole states in the low-energy region between 7 MeV and 9 MeV, and characterized by single particle transitions, a single state is found with a more distributed structure of the RRPA amplitude, exhausting approximately 2% of the EWSR. The results of the fully self-consistent RHB+RQRPA calculation, with the NL3+D1S combination of effective interactions, are shown in Figs. 8 and 9: the isovector dipole strength functions of the Sn isotopes (left panels), and the corresponding proton and neutron transition densities for the main peaks in the low-energy region (right panels). With the increase of the number of neutrons a relatively strong peak appears below 10 MeV, characterized by the dynamics of the pygmy resonance (see the transition densities). The low-energy pygmy peak is most pronounced in 124Sn. It does not become stronger by further increasing the neutron number, and additional fragmentation of the low-lying strength is observed in 132Sn. For the Sn isotopes we can compare the RHB+RQRPA results with experimental data on IV GDR. In the upper panel of Fig. 10 the experimental IVGDR excitation energies Ber.75  are shown in comparison with the calculated $`E_{\mathrm{GDR}}`$. The energy of the resonance is defined as the centroid energy

$$\overline{E}=\frac{m_1}{m_0},$$
(60)
with the energy weighted moments for discrete spectra

$$m_k=\underset{\nu }{}B(J,\omega _\nu )E_\nu ^k.$$
(61)
For $`k=1`$ this equation defines the energy weighted sum rule (EWSR). The calculated energies of the IV GDR are in excellent agreement with experimental data, and the mass dependence of the excitation energies is reproduced in detail. In the middle panel of Fig. 10 we plot the calculated energies of the pygmy states. In comparison with the IV GDR, the excitation energies of the pygmy states decrease more steeply with the mass number. The ratio of the energy weighted m1 moments calculated in the low (E$``$10 MeV) and high (E$`>`$10 MeV) energy regions, as function of the mass number, is plotted in the lower panel of Fig. 10. The relative contribution of the low-energy region increases with the neutron excess. The ratio m1,LOW/ m1,HIGH reaches a maximum $`0.06`$ for 124Sn, and it slowly decreases to $`0.05`$ for 132Sn.

## V Summary

In this work we have formulated the relativistic QRPA in the canonical single-nucleon basis of the relativistic Hartree-Bogoliubov (RHB) model. The RHB model presents the relativistic extension of the Hartree-Fock-Bogoliubov framework, and it provides a unified description of mean-field and pairing correlations. A consistent and unified treatment of the $`ph`$ and $`pp`$ channels is very important for weakly bound nuclei far from stability. In the RHB framework the ground state of a nucleus can be written either in the quasiparticle basis as a product of independent quasi-particle states, or in the canonical basis as a highly correlated BCS-state. By definition, the canonical basis diagonalizes the density matrix and it is always localized. It describes both the bound states and the positive-energy single-particle continuum. The QRPA model employed in this work is fully self-consistent. For the interaction in the particle-hole channel effective Lagrangians with nonlinear meson self-interactions are used, and pairing correlations are described by the pairing part of the finite range Gogny interaction. Both in the $`ph`$ and $`pp`$ channels, the same interactions are used in the RHB equations that determine the canonical quasiparticle basis, and in the matrix equations of the RQRPA. This is very important, because the energy weighted sum rules are only satisfied if the pairing interaction is consistently included both in the static RHB and in the dynamical RQRPA calculations. The two-quasiparticle configuration space includes states with both nucleons in the discrete bound levels, states with one nucleon in the bound levels and one nucleon in the continuum, and also states with both nucleons in the continuum. The RQRPA configuration space includes the Dirac sea of negative energy states. In addition to the configurations built from two-quasiparticle states of positive energy, the RQRPA configuration space contains pair-configurations formed from the fully or partially occupied states of positive energy and the empty negative-energy states from the Dirac sea. The inclusion of configurations built from occupied positive-energy states and empty negative-energy states is essential for the decoupling of spurious states.

The RHB+RQRPA approach has been tested in the example of the isoscalar monopole, isovector dipole and isoscalar quadrupole excitations of 22O. The NL3 parameterization has been used for the RMF effective Lagrangian, and the Gogny D1S finite range interaction has been employed in the $`pp`$ channel. In the present numerical implementation the RHB eigenvalue equations, the Klein-Gordon equations for the meson fields, and the RQRPA matrix equations are solved by expanding the nucleon spinors and the meson fields in a basis of eigenfunctions of a spherical harmonic oscillator. The calculations have illustrated the importance of a consistent treatment of pairing correlations in the RHB+RQRPA framework. The results have been compared with calculations performed in the non-relativistic continuum QRPA based on the coordinate state representation of the HFB framework. It has been shown that the RHB+RQRPA results are in agreement with recent experimental data and with non-relativistic QRPA calculations of the multipole response of neutron rich oxygen isotopes.

The RHB+RQRPA has been employed in the analysis of the evolution of the low-lying isovector dipole strength in Sn isotopes and N=82 isotones. The analysis is motivated by very recent data on the concentration of electric dipole strength below the neutron separation energy in $`N=82`$ semi-magic nuclei. It has been shown that in neutron rich nuclei a relatively strong peak appears in the dipole response below 10 MeV, with a QRPA amplitude characterized by a coherent superposition of many neutron quasiparticle configurations. The dynamics of this state corresponds to that of a pygmy dipole resonance: the oscillation of the skin of excess neutrons against the core composed of an equal number of protons and neutrons. It should be emphasized that, even though the IV GDR excitation energies calculated with the NL3 effective interaction are in excellent agreement with experimental data on Sn isotopes, the pygmy peaks in the low-energy region do not compare too well with the data on low-lying dipole strength in $`N=82`$ isotones. The calculated peaks are $``$ 2 MeV higher than the experimental weighted mean energies. This might indicate that there are problems with the isovector channel of the effective interaction and with the effective mass. Namely, if the pygmy resonance is directly related to the thickness of the neutron skin, the splitting between the excitation energies of the pygmy state and the IV GDR should be determined by the isovector channel of the effective force. A detailed quantitative analysis of the empirical low-lying isovector dipole response of neutron rich $`N=82`$ nuclei in the RHB+RQRPA framework will be included in a forthcoming publication.

Summarizing, the relativistic QRPA formulated in the canonical basis of the RHB model represents a significant contribution to the theoretical tools that can be employed in the description of the multipole response of unstable weakly bound nuclei far from stability.

## VI Appendix: Numerical details of the solution of the RQRPA equations in the canonical basis

The relativistic quasi-particle RPA equations can be simplified considerably by employing the canonical basis. According to the theorem of Bloch and Messiah Blo.62 , any RHB wave function can be expressed either in the quasiparticle basis as a product of independent quasiparticle states, or in the canonical basis as a highly correlated BCS-state. For systems with an even number of particles we have

$$|\mathrm{\Phi }=\underset{\kappa >0}{}(u_\kappa +v_\kappa a_\kappa ^{}a_{\overline{\kappa }}^{})|.$$
(62)
$`|`$ denotes the nucleon vacuum, the operators $`a_\kappa ^{}`$ and $`a_{\overline{\kappa }}^{}`$ create nucleons in the canonical basis. The occupation probabilities are given by

$$v_\kappa ^2=\frac{1}{2}\left(1\frac{\epsilon _\kappa m\lambda }{\sqrt{(\epsilon _\kappa m\lambda )^2+\mathrm{\Delta }_\kappa ^2}}\right).$$
(63)
$`\epsilon _\kappa =\kappa |\widehat{h}_D|\kappa `$ and $`\mathrm{\Delta }_\kappa =`$ $`\kappa |\widehat{\mathrm{\Delta }}|\overline{\kappa }`$ are the diagonal elements of the Dirac single-particle Hamiltonian and the pairing field in the canonical basis, respectively. In contrast to the BCS framework, however, neither of these fields is diagonal in the canonical basis. The basis itself is specified by the requirement that it diagonalizes the single-nucleon density matrix $`\widehat{\rho }(𝐫,𝐫^{})=_{k>0}V_k(𝐫)V_k^{}(𝐫^{})`$. The transformation to the canonical basis determines the energies and occupation probabilities of single-nucleon states, that correspond to the self-consistent solution for the ground state of a nucleus. Since it diagonalizes the density matrix, the canonical basis is localized. It describes both the bound states and the positive-energy single-particle continuum Dob.96 .

Many of the eigenvalues (63) of the density matrix are identically zero. In particular, those at very high energies in the continuum, but also those that correspond to the levels in the Dirac sea (no-sea approximation). Because of this degeneracy the levels in the canonical basis are not uniquely determined by the numerical diagonalization of the density matrix $`\widehat{\rho }(𝐫,𝐫^{})`$. In addition to the well defined eigenstates $`|\kappa `$ with non-degenerate eigenvalues $`0<v_\kappa ^2<1`$, there is one set of eigenstates with eigenvalues equal 0 and another set of eigenstates with eigenvalues equal 1. Any linear combination of eigenstates with eigenvalue 0 (1) is again an eigenstate with eigenvalue 0 (1). The diagonal pairing matrix elements $`\mathrm{\Delta }_\mu `$ vanish in these degenerate subspaces. The corresponding single particle energies $`\epsilon _\mu `$, however, are arbitrary and unphysical. Within these two subspaces the canonical basis is not uniquely defined.

We therefore introduce an additional requirement, that the canonical basis in each of these subspaces diagonalizes the single particle Hamiltonian $`\widehat{h}_D`$. In practical applications one thus first diagonalizes the matrix $`\widehat{\rho }`$. This gives all the canonical basis states with $`0<v_\kappa ^2<1`$, and in addition two sets of degenerate eigenstates with eigenvalues 0 and 1. Two eigenstates $`|\kappa `$ and $`|\lambda `$ are considered degenerate if the corresponding eigenvalues differ less than a given parameter $`ϵ_d`$:

$$|v_\kappa ^2v_\lambda ^2|<ϵ_d.$$
(64)
In the second step the single particle Hamiltonian $`\widehat{h}_D`$ is diagonalized in the subspace of degenerate eigenvectors of the density matrix with eigenvalues 0 (1). These new vectors are also eigenvectors of $`\widehat{\rho }`$ with eigenvalues 0 (1). This procedure uniquely determines the energies $`\epsilon _\kappa `$ and occupation probabilities $`v_\kappa ^2`$ of single-particle states, that correspond to the self-consistent solution for the ground state of a nucleus. An appropriate choice, of course, has to be made for the parameter $`ϵ_d`$. If it is too large, a linear combination of the eigenstates $`|\kappa `$ and $`|\lambda `$ that diagonalizes $`\widehat{h}_D`$, will no longer be an eigenvector of the density matrix $`\widehat{\rho }`$.

It is important to illustrate how the RQRPA results depend on the choice of the parameter $`ϵ_d`$ in Eq. (64). For the nucleus 22O, in Fig. 11(a) we display the isovector dipole strength distributions, calculated with $`ϵ_d=10^410^7`$. For any two values of $`ϵ_d>10^6`$ the corresponding strength distributions show pronounced differences. When $`ϵ_d10^6`$, the dipole response does not depend any longer on its precise numerical value, and the spurious Nambu-Goldstone $`1^{}`$ mode is found at an excitation energy $`0.1`$ MeV.

The RQRPA matrix is diagonalized in the finite dimensional two-quasiparticle ($`2qp`$) vector space. There are two types of $`2qp`$ states: 1) those built from $`qp`$ states of positive energy, and 2) those formed by one fully or partially occupied state of positive energy and one empty negative energy state from the Dirac sea. The dimension of the RQRPA configuration space is thus determined by two cut-off parameters: $`E_{Cp}`$ is the maximum value of the sum of the diagonal matrix elements of $`H^{11}`$ (54) for the first type of $`2qp`$ states, and $`E_{Ca}`$ is the maximum absolute value of the sum of the diagonal matrix elements of $`H^{11}`$ (54) for $`2qp`$ states with one quasiparticle in the Dirac sea. The choice of the two cut-off parameters $`E_{Cp}`$ and $`E_{Ca}`$ is restricted by the following conditions: (a) there should be no response to the number operator, i.e. the Nambu-Goldstone $`0^+`$ mode associated with the nucleon number conservation should have zero excitation energy, (b) the spurious excitation corresponding to the translation of the nucleus decouples as a zero-energy excitation mode, and (c) the response function does not depend on the precise numerical values of $`E_{Cp}`$ and $`E_{Ca}`$.

In Fig. 11(b) we show how the response to the neutron number operator for 22O varies with the cut-off parameter $`E_{Cp}`$ in the range 30 – 270 MeV. The choice $`E_{Ca}`$=1700 MeV includes the entire negative-energy Dirac spectrum. The response is obviously reduced as the number of $`2qp`$ configurations increases. Already for $`E_{Cp}=`$ 90 MeV the Nambu-Goldstone $`0^+`$ mode converges to $`0.1`$ MeV.

A large configuration space is also necessary in order to bring the spurious $`1^{}`$ state at zero excitation energy. In Fig. 11(c) we illustrate the convergence of the energy of the $`1^{}`$ spurious state in 22O and 120Sn. The excitation energies are plotted as functions of the energy cut-off parameter $`E_{Cp}`$. $`E_{Ca}`$ is kept at 1700 MeV.

The choice of the cut-off parameter $`E_{Ca}`$ has a pronounced influence on the calculated isoscalar monopole response. This is illustrated in Fig. 11(d), where we show how the energies of the giant monopole resonance (GMR) in 22O and 120Sn depend on the value of $`E_{Ca}`$. For $`E_{Ca}`$ 1150 MeV, only positive-energy $`2qp`$ states are included in the RQRPA basis and the excitation energies of the GMR peaks are simply too low. As $`E_{Ca}`$ is increased to include the negative energy states, the GMR excitation energies also increase and saturate for $`E_{Ca}`$1500 MeV.

ACKNOWLEDGMENTS

This work has been supported in part by the Bundesministerium für Bildung und Forschung under project 06 TM 979, and by the Gesellschaft für Schwerionenforschung (GSI) Darmstadt. T. N. acknowledges the support from the Alexander von Humboldt - Stiftung.

==========================================================================