Introduction

The study of nuclei far from stability is an increasingly important part of nuclear physics [1,2]. As radioactive beams allow more experiments on these nuclei, theoretical modeling is changing in significant ways. New ideas and progress in computer technology have allowed nuclear theorists to understand bits and pieces of nuclear structure quantitatively [3]. Short-lived exotic nuclei offer unique tests of those aspects of our developing many-body theories that depend on neutron excess [4]. The major challenge is to predict or describe in detail exotic new properties of nuclei far from the stability valley, and to understand the origins of these properties.

For medium-mass and heavy nuclei, an important goal is obtaining a universal energy-density functional, which will be able to describe static and dynamic properties of finite nuclei and extended nucleonic matter with arbitrary neutron-to-proton ratio. Self-consistent methods based on density-functional theory are already sophisticated enough to allow precise analysis of ground-state properties (e.g. binding energies) in heavy nuclei [5,6,7]. They can also help describe nuclear decays and excited states. Their predictions for collective excitations as we approach the neutron drip line are especially interesting. What happens to low- and high-frequency multipole modes when the neutron excess is unusually large?

To address these questions we use the quasiparticle random-phase approximation (QRPA), a powerful tool for understanding both low-lying vibrational states and giant resonances [8]. The QRPA is a microscopic approach that is nevertheless simple enough to allow ``no-core'' calculations. The approximation, which should be good for collective vibrations as long as their amplitudes are small, is especially effective in conjunction with Skyrme energy functionals. Our work is a part of a broad program to test and improve these functionals, which thus far have been fitted mainly to ground-state observables, by applying them to collective excitations, particularly near the drip line. This paper lays out our approach and evaluates its accuracy. For these purposes we restrict ourselves to a single Skyrme functional, SkM$^*$. A forthcoming study will examine the performance of Skyrme functionals more generally.

The QRPA is a standard method for describing collective excitations in open-shell superconducting nuclei with stable mean-field solutions, either spherical or deformed. What is not standard, and at the same time is extremely important for weakly bound nuclei, is the treatment of the particle continuum. Continuum extensions of the random phase approximation (RPA) or QRPA are usually carried out in coordinate space, facilitating treatment of decay channels and guaranteeing correct asymptotics. Surprisingly, as we discuss below, the rich literature on the RPA and QRPA, which includes many coordinate-space calculations, contains few treatments of the continuum that exploit the entire Skyrme functional in a fully self-consistent way.

To avoid confusion, we state what we mean by a fully self-consistent RPA or QRPA calculation. First, the underlying mean-field calculation must be self-consistent in the usual sense. Next, the residual interaction used in the RPA or QRPA must be derived from the same force or energy functional that determines the mean field. An important consequence of this condition, and of other more detailed technical conditions discussed below, is that spurious excitations arising from symmetry breaking by the mean field have zero or nearly zero energy, leaving the physical intrinsic excitations completely uncontaminated by spurious motion. Finally, energy-weighted sum rules must be satisfied to high accuracy. We elaborate on these requirements below; Refs. [9,10,11] discuss ways in which RPA calculations commonly violate them.

The literature applying RPA or QRPA to nuclear structure is huge, and a complete review is beyond the scope of our paper. We do, however, present an overview of the studies that are related in one way or another to nuclear density functionals, self consistency, pairing, and the key issue of the particle continuum.

The standard version of QRPA, the so-called matrix formulation, is carried out in the configuration space [12,13] of single-quasiparticle states. A number of papers treat collective states in spherical nuclei in the Skyrme-RPA and QRPA matrix formulation (see Refs. [13,14] and references cited therein), in which the positive-energy continuum is discretized, e.g. by solving the Hartree-Fock-Bogoliubov (HFB) and QRPA equations in a harmonic-oscillator single-particle basis. Within this group, the first fully self-consistent calculations that properly account for continuum effects are those of Refs. [15,16], in which the localized canonical basis of coordinate-space HFB is used to calculate beta-decay rates of neutron-rich r-process nuclei and Gamow-Teller strength distributions. Recently, fully self-consistent HFB+QRPA calculations have also been carried out with the finite-range Gogny force [17]. Unlike many previous Gogny+HFB studies that employed a harmonic oscillator basis, Ref. [17] solves the HFB equations in the eigenbasis of a Woods-Saxon potential, the particle continuum of which is discretized by enclosing the system in a box.

Coordinate-space Green's functions as a method of implementing the RPA through linear response were first used in Ref. [18] and subsequently applied to the description of low- and high-energy nuclear modes (see, e.g., Refs. [19,20,21,22,23,24,25,26,27,9,28]). Many of those calculations are not realistic enough, however, because they ignore the spin-orbit and Coulomb residual interactions in the RPA [10,11]. Coordinate-space Green's-function QRPA was studied in Ref. [29], in the BCS approximation, with a phenomenological Woods-Saxon average potential. Coordinate-space HFB+QRPA for spherical nuclei was formulated in Refs. [30,31,32,33] and applied to excitations of neutron-rich nuclei. As in [29], the Hartree-Fock (HF) field in Refs. [30,31] was approximated by a Woods-Saxon potential. While the calculations of Refs. [32,33] are based on Skyrme-HFB fields, they violate full self consistency by replacing the residual velocity-dependent terms of the Skyrme force by the Landau-Migdal force in the QRPA, and neglecting spin-spin, spin-orbit, and Coulomb residual interactions entirely. Within this approach, extensive Skyrme-HF+BCS QRPA calculations of E1-strength in neutron-rich nuclei were carried out in Refs. [34,35].

An alternative coordinate-representation approach, also based on Green's functions, was formulated in Refs. [36,37] within Migdal's finite-Fermi-systems theory. Most of practical applications of this method, however, involve approximations that break self consistency in one way or another, including the use of highly truncated pairing spaces, different interactions in HFB and QRPA, and the so-called diagonal pairing approximation [36,38,39,40,41,42,43,44,45]. Properties of excited states and strength functions have also been investigated within the relativistic RPA [46,47,48,49,50,51,52] or QRPA [53,54]. The QRPA work employs the matrix formulation and is fully self-consistent, since it uses the same Lagrangian in the relativistic Hartree-Bogoliubov calculation of the ground state and in the QRPA matrix equations, which are solved in the canonical basis.

At the present, no fully self-consistent continuum HFB+QRPA calculations exist in deformed nuclei. Refs. [55,56] studied giant resonances in deformed nuclei within time-dependent HF theory, formulated in coordinate space with a complex absorbing boundary condition imposed. Symmetry-unrestricted RPA calculations, with no pairing, were carried out in Ref. [57] in a ``mixed representation'' [58] on a Cartesian mesh in a box, while Ref. [59] contains examples of BCS+QRPA calculations in the single-particle basis of a deformed Woods-Saxon potential.

The work described in this paper is fully self-consistent: among other things we use precisely the same interaction in the HFB and QRPA calculations so as to preserve the small-amplitude limit of time-dependent HFB. We formulate the QRPA in the canonical eigenbasis of the one-body particle-density matrix [60] which is calculated in the coordinate representation in a large spherical box. As mentioned above, the canonical basis has been used previously to study $\beta$ decay and Gamow-Teller strength [15,16]; its use in charge-conserving modes near the drip line is more challenging, however, because of the existence of spurious states in the monopole and dipole channels.¹ These zero-energy modes can mix with physical states unless the QRPA equations are solved with high accuracy. A less precise implementation of our approach was used to calculate neutrino-nucleus cross sections in $^{208}$Pb in Ref. [62].

This paper is organized as follows. Section 2 below presents our approach. In Sec. 3 we check the QRPA solutions carefully, focusing on spurious modes. Section 4 contains the main conclusions of our work. Mathematical details are in two appendices, the first of which is on the QRPA equations and the second on calculating the derivatives of the Skyrme functionals that enter the formalism.