CDCC calculations with the Lagrange-mesh technique

Nuclear Physics A 845 (2010) 88–105

www.elsevier.com/locate/nuclphysa

CDCC calculations with the Lagrange-mesh technique

T. Druet a,b, D. Baye a,b, P. Descouvemont b,∗,1, J.-M. Sparenberg a,b

a Physique Quantique, C.P. 165/82, Université Libre de Bruxelles (ULB), B 1050 Brussels, Belgiumb Physique Nucléaire Théorique et Physique Mathématique, C.P. 229, Université Libre de Bruxelles (ULB),

B 1050 Brussels, Belgium

Received 7 January 2010; received in revised form 19 May 2010; accepted 28 May 2010

Available online 3 June 2010

Abstract

We apply the Lagrange-mesh technique to the Continuum Discretized Coupled Channel (CDCC) theory.The CDCC equations are solved with the R-matrix method, using Lagrange functions as variational basis.The choice of Lagrange functions is shown to be efficient and accurate for elastic scattering as well as forbreakup reactions. We describe the general formalism for two-body projectiles, and apply it to the d + 58Nicollision at Ed = 80 MeV. Various numerical and physical aspects are discussed. Benchmark calculationson elastic scattering and breakup are presented.© 2010 Elsevier B.V. All rights reserved.

Keywords: CDCC method; R-matrix theory; Lagrange meshes; d + 58Ni scattering and breakup

1. Introduction

The Continuum Discretized Coupled Channel (CDCC) method has been introduced more thanthirty years ago [1] to describe deuteron induced reactions. Owing to the low binding energy ofthe deuteron, it was shown that including continuum channels significantly improves the descrip-tion of d + nucleus elastic cross sections [1,2]. The main goal of the CDCC method is to solve

* Corresponding author.E-mail addresses: [email protected] (T. Druet), [email protected] (D. Baye), [email protected] (P. Descouvemont),

[email protected] (J.-M. Sparenberg).1 Directeur de Recherches FNRS.

0375-9474/$ – see front matter © 2010 Elsevier B.V. All rights reserved.doi:10.1016/j.nuclphysa.2010.05.060

T. Druet et al. / Nuclear Physics A 845 (2010) 88–105 89

the Schrödinger equation for reactions where the projectile presents a cluster structure with a lowdissociation energy.

After the pioneering work of Rawitscher on d + 40Ca, the CDCC method was extensively de-veloped and used by several groups [3–7]. Many improvements to the original formulation havebeen proposed: three-body projectiles [8–10], inclusion of core excited states [11], applicationto breakup [2,12], reaction [13] and fusion [14]. The recent development of intense radioac-tive beams strongly revived the need for accurate reaction theories. The CDCC method providesa solution which, in principle, can be very close to the exact wave function, and is thereforewell adapted to low-energy reactions. Continuum states of the projectile are either described bysquare-integrable wave functions (“pseudostate” method, see Ref. [3]) or by averages of exactscattering wave functions (“bin” method, see Refs. [15,12]). In both cases the Schrödinger equa-tion is replaced by a set of coupled differential equations, which provides the physical quantitiessuch as the collision matrix.

In the present work, we propose a new approach to solve the CDCC equations. This methodis based on the R-matrix theory [16,17] associated with a Lagrange-mesh basis [18]. The R-matrix method is widely used in atomic and nuclear physics. It is based on a separation of theconfiguration space into two regions: the internal region, where the wave function is expandedover a basis, and the external region, where the interaction is limited to the Coulomb force, andwhere the wave function has reached its asymptotic behaviour. The matching between internaland external solutions provides the collision matrix, and therefore various cross sections.

In principle, the basis aimed at describing the internal wave function can be chosen freely.Here we again adopt a Lagrange basis, which was shown to be very efficient in bound-state [18,19] and scattering [20,21] calculations. We will show that the combination of both approachesprovides a fast and accurate technique to solve the CDCC equations, even for large systems,where traditional methods meet convergence problems [22,4].

In Section 2, we briefly present the CDCC formalism. Section 3 is devoted to the R-matrixtheory applied to the CDCC equations, with the Lagrange basis. The method is tested in Section 4with the d + 58Ni elastic scattering and breakup at Ed = 80 MeV. Concluding remarks andoutlook are presented in Section 5.

2. The CDCC formalism

The CDCC theory has been known for many years. We briefly present here the general for-malism, and derive the differential system to be solved. We refer the reader to Refs. [3,15] formore detail. We assume here a two-body projectile made of a core c and a fragment f , with r asrelative coordinate. The channel spin s is obtained from the coupling of the core and fragmentspins. As usual in CDCC calculations, the internal structure of the particles is neglected. Thetwo-body Hamiltonian associated with the projectile is therefore written as

H0 = Tr + Vcf (r), (1)

where potential Vcf (r) is in general real, but may contain a spin–orbit term. We limit the pre-sentation to local potentials for the sake of simplicity, but the model can be applied to non-localpotentials as well. In Eq. (1), Tr is the kinetic energy between the core and the fragment.

Let us consider the core-fragment Schrödinger equation

H0Φjm

(r) = EjΦ

jm(r), (2)
l l l

90 T. Druet et al. / Nuclear Physics A 845 (2010) 88–105

which presents a finite set of bound-state solutions (Ejl < 0) and a continuum spectrum (Ej

l > 0).In Eq. (2), the wave function with total angular momentum j is written as

Φjml (r) = 1

rφ

jl (r)il

[Yl(Ωr) ⊗ χs

]jm, (3)

where χs is a spinor associated with the projectile spin.Let us now consider the three-body system involving the spinless target with mass At and

charge Zt , and the projectile with mass Ap and charge Zp . The three-body Hamiltonian reads

H = H0 + TR + Vtc

(R + Af

Ap

r

)+ Vtf

(R − Ac

Ap

r

), (4)

where Ac and Af are the core and fragment masses (Ap = Ac +Af ), R = (R,ΩR) is the relativecoordinate between the target and projectile centers of mass, and (Vtc,Vtf ) are target-core andtarget-fragment optical potentials. In the CDCC formalism, the total wave function with totalangular momentum J and parity π is expanded over two-body solutions (3) as

Ψ JMπ(R, r) = 1

rR

∑ljLi

uJπljLi(R)φ

jli(r)Y

JMljL (ΩR,Ωr), (5)

where

YJMljL (ΩR,Ωr) = il+L

[[Yl(Ωr) ⊗ χs

]j ⊗ YL(ΩR)]JM

. (6)

In Eq. (5), index i stands for two-body bound states (with energy Ejli < 0) and for continuum

states. The treatment of the breakup channels can be done in two ways:

• The pseudostate method (PS). In this variant, functions φjli(r) are square-integrable positive-

energy (Ejli > 0) approximate solutions of Eq. (2). The two-body Schrödinger equation (2)

is solved over a finite basis, which provides a discrete spectrum.• The bin method (bin). The alternative is to construct φ

jli(r) by averaging continuum wave

functions φjl (kr , r) with wave number kr [23,24]. The two-body continuum is truncated at

some wave number kmax, and the [0, kmax] interval is divided in Nb subintervals [ki−1, ki]with width �i . The basis functions are defined as

φjli(r) =

√1

�i

ki∫ki−1

φjl (kr , r) dkr , (7)

where φjl (kr , r) presents the asymptotic behaviour

φjl (kr , r) →

√2

π

[cos δ

jl (kr )Fl(krr) + sin δ

jl (kr )Gl(krr)

], (8)

in order to be normalized such as⟨φ

jl (kr , r)

∣∣φjl

(k′r , r

)⟩ = δ(kr − k′

r

). (9)

In Eq. (8), Fl(x) and Gl(x) are the regular and irregular Coulomb functions, respectively,and δ

jl (kr ) is the phase shift. In definition (7) various weight functions (taken here as unity)

can be chosen according to the nature of the problem [12].


The PS and bin functions (7) are both square-integrable, and are normalized as⟨φ

jli

∣∣φj

li′⟩ = δii′ . (10)

The three-body wave functions uJπγ i (R) are then determined from the Hamiltonian (4) (in-

dex γ stands for γ = ljL). The central potentials Vtc and Vtf are first expanded numerically inmultipoles using

Vtc

(R + Af

Ap

r

)+ Vtf

(R − Ac

Ap

r

)=

∑λ

Vλ(r,R)Pλ(cos θRr), (11)

where θRr is the angle between R and r , and Pλ(x) a Legendre polynomial. This provides thecoupling potentials as

V Jπγ i,γ ′i′(R) = ⟨

φjli(r)Y

JMγ (ΩR,Ωr)

∣∣Vtc + Vtf

∣∣φj ′l′i′(r)Y

JMγ ′ (ΩR,Ωr)

⟩, (12)

which represent five-dimensional integrals. The four angular integrals are performed analytically,and the integral over the core-fragment coordinate r requires a numerical approach. As we willsee below [Eq. (33)], the choice of Lagrange functions as basis makes this integral very simple,since it does not require any numerical quadrature.

Using the multipole expansion (11), the potentials (12) can be written as

V Jπγ i,γ ′i′(R) = (−1)J−s il

′+L′−l−Lj j ′ l l′LL′ ∑λ

(−1)λ(

l l′ λ

0 0 0

)(L L′ λ

0 0 0

)

×{

j L J

L′ j ′ λ

}{j l s

l′ j ′ λ

} ∞∫0

φjli(r)Vλ(r,R)φ

j ′l′i′(r) dr, (13)

where x = √2x + 1. When the projectile has an internal spin s = 0, this equation is simplified.

Finally, the relative functions uJπγ i (R) are determined from a set of Nc coupled equations[

− h2

2μ

(d2

dR2− L(L + 1)

R2

)+ E

jli − E

]uJπ

γ i (R) +∑γ ′i′

V Jπγ i,γ ′i′(R)uJπ

γ ′i′(R) = 0, (14)

where μ is the projectile-target reduced mass. The number of equations Nc is obtained from theγ values and from the number of pseudostates (or bins). Of course the sum over the projectileangular momentum l must be truncated at some value lmax. However all λ values, consistent withangular-momentum couplings, are included in the summation.

3. The R-matrix method on a Lagrange mesh

3.1. Introduction

The system (14) is common to all CDCC calculations, but various methods have been pro-posed in the literature to solve it. For example, Ref. [24] uses the method proposed by Ichimuraet al. [25], based on an iterative procedure, and well adapted to long-range potentials. Nunesand Thompson [4] determine a finite-difference solution [22,26] of the system (14). In Ref. [13],Huu-Tai uses the B-spline method [27] to solve the CDCC system.

In the present work, we adapt the R-matrix theory [16,28,17] to the CDCC formalism. TheR-matrix method is well known in atomic and nuclear physics. It represents an efficient tool to


extend variational calculations, performed on a limited basis, to scattering problems. The methodis numerically stable and can be used for large systems.

3.2. The R-matrix method

The R-matrix theory has been developed for more than sixty years, first in nuclear physics[29] and then essentially in atomic physics [30]. The idea is to split the configuration space intwo regions: the internal region (with radius a), where the wave function is expanded over a finiteset of basis functions, and the external region, where the wave function has reached its Coulombasymptotic behaviour. Since the kinetic-energy operator is not Hermitian over a finite interval,the R-matrix theory makes use of the Bloch operator [31]

L =∑

c

|c〉Lc〈c|, (15)

where the channel index c stands for c = γ i, and where channel wave functions |c〉 are definedas

|c〉 = 1

rφ

jli(r)Y

JMγ (ΩR,Ωr). (16)

In Eq. (15), the partial Bloch operator Lc is given by

Lc = h2

2μδ(R − a)

(d

dR− Bc

R

), (17)

where Bc is a (real) boundary parameter, taken here as Bc = 0. The Bloch operator acts at thesurface (R = a) only, and makes the Hamiltonian Hermitian over the internal region. It ensuresthe continuity of the logarithmic derivative at the surface.

As mentioned before, the internal wave functions uJπc (R) are expanded over N basis functions

ϕn(R) as

uJπc,int(R) =

N∑n=1

AJπcn ϕn(R), for R � a (18)

where coefficients AJπcn are to be determined. At large distances R, the potentials in (14) tend to

the Coulomb interaction as

V Jπcc′ (R) → ZpZte

2

Rδcc′ . (19)

Accordingly, the external wave function can be expanded as a combination of Coulomb func-tions. For an entrance channel ω, the external part is given by

uJπc,ext(R) = 1√

vc

(Ic(kcR)δcω − Oc(kcR)UJπ

cω

), (20)

where vc and kc are the relative velocity and wave number in channel c, Ic(x) = Gc(x)− iFc(x)

and Oc(x) = I ∗c (x) are the incoming and outgoing Coulomb functions, and UJπ is the collision

matrix.Adding the Bloch operator on both sides of Eq. (14) provides, in a schematic notation, the

Bloch–Schrödinger system of equations

(Tc + Lc + Ec − E)uJπc,int +

∑′

V Jπcc′ uJπ

c′,int = LcuJπc,ext, (21)

c


which is solved in the internal region. As the Bloch operator acts at the surface, and since thewave functions uJπ

c (R) are continuous, the asymptotic expression (20) can be used in the r.h.s. ofEq. (21). In the l.h.s., expansion (18) is used for the internal wave function. Let us define matrixCJπ with dimension NcN as

CJπcn,c′n′ = 〈ϕn|(Tc + Lc + Ec − E)δcc′ + V Jπ

cc′ |ϕn′ 〉int, (22)

where subscript “int” means that the matrix element is calculated over the internal region only.Then, the R-matrix is obtained from

RJπcc′ = h2

2μa

∑nn′

ϕn(a)[(

CJπ)−1]

cn,c′n′ϕn′(a), (23)

which involves the basis functions at the channel radius. In general, potentials V Jπcc′ (R) are com-

plex, and the R matrix is complex as well. The collision matrix UJπ is then deduced from thecontinuity of the wave function at the channel radius as [17]

UJπ = (ZJπ

O

)−1ZJπ

I , (24)

where matrix ZJπO reads(

ZJπO

)cc′ = (kc′a)−1/2(Oc(kca)δcc′ − kc′aRJπ

cc′ O ′c′(kc′a)

), (25)

and a similar definition holds for ZJπI . Notice that for complex potentials, ZO = Z∗

I , and thecollision matrix UJπ is not unitary. The R matrix (23) depends on the channel radius a. However,the collision matrix, which determines the cross sections, should not depend on a, providedthat a is large enough to ensure condition (19), and that the basis functions ϕn(R) are able toaccurately describe the internal wave functions (18). Various choices exist for the basis functions(see Ref. [17] for a detailed discussion). Here we use Lagrange basis functions.

When the collision matrix is determined, coefficients AJπcn of the internal wave function (18)

can be obtained from the linear system∑c′n′

CJπcn,c′n′AJπ

c′n′ = ⟨ϕn

∣∣Lc

∣∣uJπc,ext

⟩. (26)

3.3. Lagrange basis

In the Lagrange-mesh method, the basis functions ϕn(R) are associated with a Gauss quadra-ture. For a finite interval [0, a], they are defined as [20]

ϕn(R) = (−1)N+n

√xn(1 − xn)

axn

RPN(2R/a − 1)

R − axn

, (27)

where PN(x) is the Legendre polynomial of degree N , and xn are the zeros of

PN(2xn − 1) = 0. (28)

These basis functions satisfy the Lagrange condition

ϕn(axm) = (aλn)−1/2δnm, (29)

where λn is the weight of the Gauss–Legendre quadrature corresponding to the [0,1] interval.If the matrix elements between basis functions (27) are computed at the Gauss approximation


of order N , consistent with the N mesh points, the calculation is strongly simplified. At thisapproximation, the overlap reads

〈ϕn|ϕn′ 〉 =a∫

0

ϕn(R)ϕn′(R)dR ≈ δnn′ , (30)

and matrix elements of a local potential V (R) are given by

〈ϕn|V |ϕn′ 〉 =a∫

0

ϕn(R)V (R)ϕn′(R)dR ≈ V (axn)δnn′ . (31)

They are simply given by the values of the potential at the mesh points. The matrix elementsof the kinetic energy are available in an analytical form [32]. Notice that the method can bedirectly applied to non-local potentials [33]. Eq. (31) is of course valid for any operator, suchas the electromagnetic operators, or the r.m.s. radius. In all cases, the calculation of the matrixelements is very fast: it does not need any numerical quadrature, and off-diagonal terms vanish.This property is quite appreciable in CDCC calculations, where computer times represent animportant issue. In spite of this simplicity, the Lagrange-mesh method is as accurate as standardvariational methods [18].

3.4. Description of two-body breakup channels

In the present work, Lagrange functions are used to describe the projectile-target motion, butalso to define the internal states of the projectile φ

jli(r) [see Eqs. (2), (3)]. This provides a further

simplification of the coupling potentials (13), since the integral over the core-fragment coordi-nate r is simply replaced by a calculation of the potential at the mesh points. For pseudostates,two different types of Lagrange functions can be used: Legendre functions, similar to (27) butinvolving the channel radius ar , or Laguerre functions [18], defined in the [0,∞] interval. Thelatter basis involves a scaling parameter h, adapted to the typical size of the system. In practice,the projectile wave functions must be independent of the precise value of this scaling parameter.

The projectile wave functions φjli(r) are therefore expanded as

φjli(r) =

Nr∑n=1

cjli,nϕn(r), (32)

where Nr is the number of mesh points associated with the coordinate r . Notice that, for theLegendre functions, Nr and the two-body channel radius ar do not need to be identical to thecorresponding values associated with coordinate R.

Coefficients cjli,n are determined from the two-body Schrödinger equation (2). For the Legen-

dre mesh, the integrals in (13) are obtained as

ar∫0

φjli(r)Vλ(r,R)φ

j ′l′i′(r) dr ≈

Nr∑n=1

cjli,nc

j ′l′i′,nVλ(arxn,R). (33)

A similar formula holds for Laguerre functions, where the integral is performed from zero toinfinity.


The treatment of breakup reactions, as well as the use of bin functions (7) require two-bodycontinuum wave functions. They can be either determined with finite-difference methods, or withthe R-matrix theory. In the R-matrix formalism, the internal wave function of the projectile atmomentum kr is expanded over Legendre functions and reads

φjl (kr , r) =

Nr∑n=1

djl,n(kr )ϕn(r), for r � ar , (34)

and the external part is given by Eq. (8). The calculation of the two-body phase shift δjl (kr )

and of coefficients djl,n(kr ) is analogous to the development of Section 3.2 [Eqs. (24) and (26),

respectively, in a single-channel problem]. With coefficients djl,n(kr ), the bin functions (7) are

defined, in the internal region, as in Eq. (32) with coefficients cjli,n given by

cjli,n =

√1

�i

ki∫ki−1

djl,n(kr ) dkr . (35)

This means that bin functions are defined on the same basis as PS functions, with different coef-ficients. The integral over the wave number kr is determined by a Simpson numerical quadrature.

This property is also useful to determine breakup cross sections. A CDCC breakup matrixelement from an entrance channel c = γ i to a breakup configuration with quantum numbers γ

and wave number kr can be approximated as [23,34]

UJπc,γ (kr ) ≈

∑i

UJπc,γ i

⟨φ

jli

∣∣φjl (kr )

⟩, (36)

where we have to determine the overlap between a pseudostate or a bin φjli and a scattering wave

function φjl (kr ). Using Eqs. (32) and (29), this overlap is given at the Gauss approximation by⟨

φjli

∣∣φjl (kr )

⟩ ≈ ∑n

(arλn)1/2c

jli,nφ

jl (kr , arxn) (37)

which is valid for any description of the scattering wave function. If both wave functions areconsistently expanded over the same Lagrange mesh, this overlap takes the simple form⟨

φjli

∣∣φjl (kr )

⟩ ≈ ∑n

cjli,nd

jl,n(kr ), (38)

with definition (34). A further simplification occurs for bin functions. If functions φjli are defined

as in Eq. (7) the matrix element (38) is given by⟨φ

jli

∣∣φjl (kr )

⟩ = {√1/�i if kr ∈ [ki−1, ki],

0 otherwise,(39)

and is therefore constant inside a bin interval. Eq. (39) is simple, but requires an interpolationprocedure. On the contrary, definition (38) provides a smooth function of the wave number kr

and does not need any further interpolation.

3.5. Extension to large systems of coupled equations

In reactions involving halo nuclei, long-range (dipole or quadrupole) Coulomb potentials arepresent, and the R-matrix radius should be extended to account for this property. At first sight


this requires the use of very large bases. This problem can be addressed within the R-matrixformalism by using propagation methods [35,17]. These developments of the R-matrix theoryare well known and widely used in atomic physics to deal with the long range of the Coulombinteraction.

For large numbers of coupled equations, the inversion of matrix CJπ (see Eq. (23)) representsthe main part of the computer time. The principle of propagation techniques is to split the [0, a]interval in several (NS ) subintervals. The number of basis functions in each subinterval is there-fore reduced by a factor NS (typically NS ≈ 5–10, but there is no limitation provided that thenumber of basis functions per interval N/NS is large enough to reproduce the wave function).A single inversion of a matrix of dimension NcN can be therefore replaced by NS inversionsof matrices of dimension (NcN/NS). In Ref. [17] we show that Lagrange functions are welladapted to propagation techniques. We give the matrix elements of the kinetic energy betweenLagrange functions in an interval [ai, ai+1] involving lower and upper boundaries ai and ai+1,respectively. The simplicity of propagation methods with Lagrange functions arises from the factthat the overlap and potential matrix elements remain as simple as in Eqs. (30) and (31).

4. Application: d + 58Ni scattering at Ed = 80 MeV

4.1. Conditions of the calculations

We apply the method to the d + 58Ni elastic scattering and breakup at Ed = 80 MeV. This sys-tem is considered in the literature as a test case. It has been investigated by several authors (seefor example Refs. [3,2,24,6,7,10]) who showed the importance of the breakup channels in theelastic cross sections. The nucleon + 58Ni central interactions are chosen as in these references(the spins of the nucleons are assumed to be zero). The p + n potential is taken as the Min-nesota force [36], fitted on the deuteron binding energy and on the low-energy nucleon–nucleonscattering properties. Our main goal here is to illustrate the R-matrix formalism in the CDCCframework. We essentially use the PS method to discretize the p + n continuum. As far as theR-matrix is concerned, the PS and bin methods are treated in an equivalent way. Throughout thepaper we use h2/2mN = 20.736 MeV fm2, e2 = 1.44 MeV fm, and use integer masses for thenucleons and 58Ni.

4.2. Influence of the p + n basis

The p+n ground state and pseudostates are described with a Laguerre mesh. Unless specifiedotherwise we use Nr = 30 and h = 0.4 fm. The scaling parameter h is optimized for the p + n

potential. However, the physical results, and in particular the ground-state energy, are insensitiveto that choice. Fig. 1 shows the PS energies up to lmax = 6 and Emax = 60 MeV. This truncationenergy approximately corresponds to kmax = 1.2 fm−1, a typical value in the literature and usedthroughout this work. Parameters lmax and kmax determine the size of the CDCC system (14).

Before applying the R-matrix to the d + 58Ni scattering, it is crucial to choose the channelradius a. This can be done by analysing the d + 58Ni potentials, and in particular their asymptoticbehaviours. Typical nuclear potentials are shown in Fig. 2 for J = 0+ (the centrifugal term is notincluded). The diagonal potential for the deuteron ground state is negligible beyond R = 10 fm.However, the potentials associated with some PS extend to larger distances (the 2+

1 potentialwhich presents the longest range is displayed). This is also true for the coupling potentials. The


Fig. 1. p + n pseudostates up to Emax = 60 MeV as a function of the deuteron angular momentum, with a Laguerremesh (Nr = 30, h = 0.4 fm).

Fig. 2. Real part (in MeV) of d + 58Ni diagonal (upper panel) and non-diagonal (lower panel) J = 0+ potentials forsome typical channels. The labels represent the quantum numbers (l,L, i; l′,L′, i′)

channel radius must be chosen such that, at relative distances R � a, all nuclear terms are neg-ligible. From Fig. 2, we conclude that the channel radius must be at least a ≈ 15 fm to deriveprecise phase shifts. The sensitivity with respect to a will be analysed later.

As we want to compare our results with the literature we use, for elastic scattering (Sec-tions 4.2 and 4.3), the same conditions, i.e. we neglect odd partial waves in the p + n relativemotion. This is justified for elastic cross sections which are hardly affected by odd partial waves.Table 1 illustrates the sensitivity of the elastic collision matrix for different p +n bases. We con-sider variations of lmax and of Nr for typical J values. The convergence with respect to lmax isvery fast; in order to save computer time lmax = 4 will be used in the calculations. For complete-ness we also display the collision matrix obtained with the bin method, for different discretization


Table 1Convergence of the elastic element UJπ

11 (multiplied by 102) as a function of the deuteron maximum angular momentumlmax (top, Nr = 30), of the basis size Nr with the PS method (middle, lmax = 4), and of the number of bins Nb (bottom,lmax = 4). The R-matrix parameters are a = 15 fm and N = 50.

102 × UJπ11 lmax = 0 lmax = 2 lmax = 4 lmax = 6

J = 0+ 4.42 5.23 5.07 5.03+9.85i +9.88i +9.87i +9.84i

J = 5− −4.40 −5.01 −4.83 −4.79−9.41i −10.34i −10.35i −10.36i

J = 20+ 86.14 86.39 86.50 86.51+22.64i +18.84i +18.56i +18.55i

PS Nr = 15 Nr = 20 Nr = 25 Nr = 30

J = 0+ 5.07 5.02 5.06 5.07+10.00i +9.83i +9.94i +9.87i

J = 5− −4.82 −4.81 −4.82 −4.83−10.34i −10.33i −10.34i −10.35i

J = 20+ 86.51 86.50 86.50 86.50+18.58i +18.54i +18.56i +18.56i

bin Nb = 15 Nb = 20 Nb = 25 Nb = 30

J = 0+ 5.19 5.12 5.09 5.08+9.94i +9.93i +9.93i +9.93i

J = 5− −4.92 −4.87 −4.85 −4.84−10.35i −10.35i −10.34i −10.34i

J = 20+ 86.46 86.48 86.49 86.49+18.59i +18.57i +18.57i +18.56i

numbers Nb. For converged values, the differences between the PS and bin methods are small(less than 1%), and have a negligible impact on the elastic-scattering cross section. Higher preci-sion could still be obtained by increasing the truncation wave number kmax (or, equivalently, theenergy Emax).

Table 1 shows that the J = 0+ partial wave is the most sensitive to the conditions of thecalculations. This is expected, since it presents the lowest Coulomb barrier. The internal partof the wave function is therefore more important than for other partial waves. The number ofchannels or, in other words, the size of the system (14) strongly depends on lmax. It varies from 15for lmax = 0 to 198 for lmax = 6. These numbers must be multiplied by the number of basisfunctions N to obtain the dimension of matrix CJπ (22).

The elastic cross sections are displayed in Fig. 3 and compared with available experimen-tal data [37]. The calculations are performed with a channel radius a = 15 fm and N = 50. Asobserved in previous calculations, the CDCC approximation provides a fair agreement with ex-periment. In particular the locations of the minima are well reproduced. Including odd partialwaves in the p +n wave functions does not affect the results at the scale of the figure. For anglesθ > 20◦, the influence of breakup channels is not negligible.

4.3. Influence of the d + 58Ni R-matrix parameters

In this subsection, we analyze the sensitivity of the collision matrix against the channel radiusa, and the number of d + 58Ni basis functions N . These two parameters define the R-matrixframework. In Table 2, we give the elastic element of the collision matrix for different conditions.


Fig. 3. Elastic cross sections at Ed = 80 MeV (divided by the Rutherford cross section) with different lmax (upper panel)and Nr (lower panel) values. Larger values are superimposed to the lmax = 4 and Nr = 25 curves. Experimental data(crosses) are taken from Ref. [37].

Increasing a requires a simultaneous increase of N , since the basis should be able to reproducethe wave function over a wider range. Consequently, calculations performed with a = 18 fm havebeen complemented with N = 60, larger than our adopted value N = 50. The stability with thechannel radius represents a strong test of the accuracy and of the correctness of the numericalcalculation. This stability is better than 10−4 rad as soon as a is large enough to satisfy theR-matrix requirements (a � 15 fm, see Section 4.2). It is even better for large J values.

The elastic cross section for different R-matrix parameters is shown in Fig. 4. As expectedfrom the numerical values of Table 2, a = 9 fm is too small to apply the R-matrix formalism,and N = 25 is not large enough to provide an accurate description of the internal wave functions.Cross sections obtained in other conditions are hardly distinguishable. This conclusion is alsovalid when using the bin method.

In Table 3, we present a non-diagonal term of the collision matrix. The final state is char-acterized by (l′ = 0,L′ = 0, i′ = 3), i.e. it corresponds to the second pseudostate of the l = 0p + n partial wave (with Nr = 30). The convergence with respect to the R-matrix parameter a isslightly less good than for the elastic component. Obtaining a high accuracy in the breakup crosssection requires larger values for the channel radius. On the other hand, N = 50 and N = 60provide identical results with the channel radii considered here.

4.4. Breakup

Let us now investigate breakup matrix elements (36) which are used to compute breakupcross sections. Here we include odd partial waves of the p + n wave functions. In contrast withelastic scattering where they are negligible, they do play a role in breakup calculations. An otherspecificity of breakup calculations is the need for larger channel radii (see Table 3), in particular


Table 2Convergence of the elastic element UJπ

11 (multiplied by 102) as a function of the basis size N (a = 15 fm) and as afunction of the channel radius (N = 50). For a = 18 fm, the bracketed values correspond to N = 60. The conditions ofthe calculations are Emax = 60 MeV, Nr = 30, h = 0.4 fm, and lmax = 4.

102 × UJπ11 N = 30 N = 40 N = 50 N = 60

J = 0+ 9.84 5.06 5.07 5.04+6.08i +9.90i +9.87i +9.88i

J = 5− −8.02 −4.83 −4.83 −4.83−10.47i −10.36i −10.35i −10.35i

J = 20+ 86.54 86.50 86.50 86.50+18.40i +18.56i +18.56i +18.56i

a = 9 fm a = 12 fm a = 15 fm a = 18 fm

J = 0+ 5.12 5.07 5.07 5.06 (5.05)

+9.91i +9.87i +9.87i +9.89 (9.87)i

J = 5− −4.86 −4.84 −4.83 −4.83 (−4.83)

−10.08i −10.36i −10.35i −10.36 (−10.34)i

J = 20+ 87.32 86.53 86.50 86.49 (86.49)

+16.99i +18.50i +18.56i +18.56 (18.56)i

Fig. 4. Elastic cross sections at Ed = 80 MeV (divided by the Rutherford cross section) with different R-matrix condi-tions: channel radius (upper panel) and basis size (lower panel).

for odd l values and at large kr . We adopt here a = 27 fm and N = 100, which provide a highaccuracy for the breakup matrix elements up to kr = 1.2 fm−1.

In Fig. 5, we present the matrix element |U (kr )|2 for J = 0+ and a final deuteron state l′ =L′ = 0. We analyze here the dependence on the p +n basis. At low momenta, this dependence isnegligible, but becomes important for kr � 0.7 fm−1. Up to kr ≈ 1 fm−1, convergence is reachedwith Nr ≈ 50, a value larger than for elastic scattering.


Table 3See caption to Table 2 for the inelastic element (l′ = 0,L′ = 0, i′ = 3) of the collision matrix (Nr = 30).

102 × UJπ13 N = 30 N = 40 N = 50 N = 60

J = 0+ 0.36 0.99 0.98 0.99−2.35i −2.44i −2.44i −2.44i

J = 5− 0.09 −0.15 −0.16 −0.16+2.21i +2.16i +2.16i +2.16i

J = 20+ −5.74 −5.74 −5.74 −5.74+0.06i +0.06i +0.06i +0.06i

a = 9 fm a = 12 fm a = 15 fm a = 18 fm

J = 0+ 0.66 0.94 0.98 1.03−2.51i −2.47i −2.44i −2.46i

J = 5− −0.09 −0.11 −0.16 −0.18+2.27i +2.13i +2.16i +2.16i

J = 20+ −5.50 −5.61 −5.74 −5.72+0.09i +0.13i +0.06i +0.04i

Fig. 5. Breakup matrix element for J = 0+ , l′ = 0, L′ = 0, and for different Nr values.

In Fig. 6, we investigate the influence of the R-matrix parameters N and a, for J = 0+ andJ = 17−, which were illustrated in Ref. [15]. As suggested by Fig. 5, we take Nr = 70 to definethe pseudostate basis. We also use the bin method (Nb = 40). Both approaches are in very goodagreement with each other, except at large kr for l′ = L′ = 2 in the J = 0+ partial wave. Theresults are stable against variations of the R-matrix parameters (channel radius and number ofbasis functions). Strong changes are necessary to observe a deviation at the scale of the figure(we take a = 18 fm and N = 60 as a comparison).

The breakup cross section for zero-spin nuclei is determined as [23]

dσBU

dkr

= π

k2

∑Jπ

∑γ

(2J + 1)∣∣UJπ

c,γ (kr )∣∣2

, (40)

where c is the entrance channel and where UJπ

is calculated from Eq. (36). This cross section ispresented in Fig. 7 with the decomposition in deuteron partial waves. We show results obtainedwith the PS method, as the bin method provides almost indistinguishable cross sections. Themain components come from l = 0,1,2, but l = 4 is not negligible at large wave numbers. Themaximum moves from kr ≈ 0.2 fm−1 for l = 0 to kr ≈ 0.6 fm−1 for l = 4. This effect comesfrom the p + n centrifugal barrier, and means that the individual contributions for large l valuesneed special attention regarding the number of pseudostates.


Fig. 6. Breakup matrix element for J = 0+ (upper row) and J = 17− (lower row), and for some l′L′ values. The solidlines are obtained with the bin method (Nb = 40) and the other lines with the pseudostate method (Nr = 70). TheR-matrix parameters are a = 27 fm and N = 100 except when mentioned otherwise.

Fig. 7. Breakup cross section (40) with the individual components of the p + n partial waves l.

4.5. Comparison with previous works

In this subsection, we test the validity and precision of our method by a comparison withprevious works [7,10]. The conditions of the calculation are slightly different from those of Sec-tion 4, and are therefore adapted. In Refs. [7,10], only the deuteron angular momenta l = 0 andl = 2 are included, and the p–n interactions are different. Moro et al. [7] use a Poeschl–Tellerpotential, Egami et al. [10] use a single-Gaussian potential, whereas we adopt the Minnesota po-tential. We have tested that the elastic-scattering and breakup cross sections are weakly sensitiveto this choice.

Fig. 8 displays the elastic cross sections of Moro et al. [7] which are perfectly reproduced bythe present model if we use, as in Ref. [7], the Poeschl–Teller potential. We also show the resultswith the Minnesota potential, adopted throughout this paper. They confirm the weak sensitivityto this potential. Only large angles are affected.

In Fig. 9 we compare the breakup matrix element (36) with Refs. [7,10] for J = 17− andvarious (l′L′) values. As for elastic scattering, the R-matrix formalism is very close to the ref-erence calculations. Here, however, the role of additional l values, in particular odd values, is


Fig. 8. Elastic cross sections of Ref. [7] (squares) compared with the present method with the Poeschl-Teller (solid line)and Minnesota (dashed line) potentials.

Fig. 9. Breakup matrix elements of Ref. [7] (squares) and of Ref. [10] (circles) compared with the present method. Solidlines correspond to l = 0 and 2 adopted in Refs. [7,10], and dashed lines to all l values up to lmax = 4.

not negligible as shown by the dashed lines. In particular the second bump in the L′ = 15 andL′ = 19′ components disappears.

5. Conclusion

In the CDCC formalism, the Schrödinger equation is replaced by a set of coupled differentialequations. This set can be very large, and solving it with a good accuracy is one of the main issuesin CDCC calculations. We have shown that the Lagrange-mesh technique, associated with the R-matrix formalism, is well adapted to this requirement. It provides a simple and efficient methodto solve the CDCC equations, and allows a common treatment of the PS and bin methods. Theformalism depends on two parameters: the channel radius a and the number of basis functions N .As usual in R-matrix calculations, the channel radius stems from a compromise: it must be largerthan the range of the nuclear force, but should be chosen as small as possible to keep limitedN values. In practice, the main part of the computation time comes from the inversion of thecomplex matrix CJπ [Eq. (22)]. In the present work the maximum size is of the order of 104. Themethod can be made still faster by using propagation techniques [35,17], where the interval [0, a]


is split in subintervals. These techniques allow to deal with large number of coupled equations,since the number of Lagrange functions can be reduced in each subinterval. Consequently thesizes of the matrices to be inverted are smaller.

The Lagrange-mesh technique is not limited to elastic scattering. We have shown that themethod can be easily used for breakup problems. The present approach can be extended in vari-ous directions. As mentioned before, the Lagrange basis is well adapted to non-local potentials,since the matrix elements, at the Gauss approximation, do not need any numerical quadrature.Consequently, non-local potentials appearing in microscopic cluster descriptions of the projectile[38] can be treated by the R-matrix formalism [33]. Dealing with non-local potentials leads to anegligible increase of the computer times. On the other hand, CDCC calculations involving three-body projectiles are being developed [8,9]. These calculations are highly time-consuming sincethe size of the system strongly increases. The simplicity of the R-matrix formalism associatedwith the Lagrange-mesh technique provides very promising perspectives in this direction.

Acknowledgements

We are grateful to P. Capel for helpful discussions. This text presents research results of theIAP programme P6/23 initiated by the Belgian-state Federal Services for Scientific, Technicaland Cultural Affairs. T.D. is partly supported by the IAP programme.

References

[1] G.H. Rawitscher, Phys. Rev. C 9 (1974) 2210.[2] M. Yahiro, Y. Iseri, H. Kameyama, M. Kamimura, M. Kawai, Prog. Theor. Phys. Suppl. 89 (1986) 32.[3] N. Austern, Y. Iseri, M. Kamimura, M. Kawai, G. Rawitscher, M. Yahiro, Phys. Rep. 154 (1987) 125.[4] F.M. Nunes, I.J. Thompson, Phys. Rev. C 59 (1999) 2652.[5] C. Beck, N. Keeley, A. Diaz-Torres, Phys. Rev. C 75 (2007) 054605.[6] O.A. Rubtsova, V.I. Kukulin, A.M. Moro, Phys. Rev. C 78 (2008) 034603.[7] A.M. Moro, J.M. Arias, J. Gómez-Camacho, F. Pérez-Bernal, Phys. Rev. C 80 (2009) 054605.[8] T. Matsumoto, E. Hiyama, K. Ogata, Y. Iseri, M. Kamimura, S. Chiba, M. Yahiro, Phys. Rev. C 70 (2004) 061601.[9] M. Rodríguez-Gallardo, J.M. Arias, J. Gómez-Camacho, R.C. Johnson, A.M. Moro, I.J. Thompson, J.A. Tostevin,

Phys. Rev. C 77 (2008) 064609.[10] T. Egami, T. Matsumoto, K. Ogata, M. Yahiro, Prog. Theor. Phys. 121 (2009) 789.[11] N.C. Summers, F.M. Nunes, I.J. Thompson, Phys. Rev. C 74 (2006) 014606.[12] J.A. Tostevin, F.M. Nunes, I.J. Thompson, Phys. Rev. C 63 (2001) 024617.[13] P.C. Huu-Tai, Nucl. Phys. A 773 (2006) 56.[14] S. Hashimoto, K. Ogata, S. Chiba, M. Yahiro, Prog. Theor. Phys. 122 (2009) 1291.[15] M. Kamimura, M. Yahiro, Y. Iseri, S. Sakuragi, H. Kameyama, M. Kawai, Prog. Theor. Phys. Suppl. 89 (1986) 1.[16] A.M. Lane, R.G. Thomas, Rev. Mod. Phys. 30 (1958) 257.[17] P. Descouvemont, D. Baye, Rep. Prog. Phys. 73 (2010) 036301.[18] D. Baye, Phys. Status Solidi (b) 243 (2006) 1095.[19] P. Descouvemont, C. Daniel, D. Baye, Phys. Rev. C 67 (2003) 044309.[20] M. Hesse, J.M. Sparenberg, F. Van Raemdonck, D. Baye, Nucl. Phys. A 640 (1998) 37.[21] P. Descouvemont, E.M. Tursunov, D. Baye, Nucl. Phys. A 765 (2006) 370.[22] L.D. Tolsma, G.W. Veltkamp, Comput. Phys. Commun. 40 (1986) 233.[23] M. Yahiro, M. Nakano, Y. Iseri, M. Kamimura, Prog. Theor. Phys. 67 (1982) 1467.[24] R.A.D. Piyadasa, M. Kawai, M. Kamimura, M. Yahiro, Phys. Rev. C 60 (1999) 044611.[25] M. Ichimura, M. Igarashi, S. Landowne, C.H. Dasso, B.S. Nilsson, R.A. Broglia, A. Winther, Phys. Lett. B 67

(1977) 129.[26] I.J. Thompson, Comput. Phys. Rep. 7 (1988) 167.[27] H. Bachau, E. Cormier, P. Decleva, J.E. Hansen, F. Martin, Rep. Prog. Phys. 64 (2001) 1815.[28] R.F. Barrett, B.A. Robson, W. Tobocman, Rev. Mod. Phys. 55 (1983) 155.


[29] E.P. Wigner, Phys. Rev. 70 (1946) 606.[30] P.G. Burke, K.A. Berrington (Eds.), Atomic and Molecular Processes: An R-matrix Approach, Institute of Physics,

Bristol, 1993.[31] C. Bloch, Nucl. Phys. 4 (1957) 503.[32] D. Baye, J. Goldbeter, J.-M. Sparenberg, Phys. Rev. A 65 (2002) 052710.[33] M. Hesse, J. Roland, D. Baye, Nucl. Phys. A 709 (2002) 184.[34] T. Matsumoto, T. Kamizato, K. Ogata, Y. Iseri, E. Hiyama, M. Kamimura, M. Yahiro, Phys. Rev. C 68 (2003)

064607.[35] K.L. Baluja, P.G. Burke, L.A. Morgan, Comput. Phys. Commun. 27 (1982) 299.[36] D.R. Thompson, M. LeMere, Y.C. Tang, Nucl. Phys. A 286 (1977) 53.[37] E.J. Stephenson, C.C. Foster, P. Schwandt, D.A. Goldberg, Nucl. Phys. A 359 (1981) 316.[38] Y. Suzuki, H. Matsumura, M. Orabi, Y. Fujiwara, P. Descouvemont, M. Theeten, D. Baye, Phys. Lett. B 659 (2008)

160.

Documents

CDCC calculations with the Lagrange-mesh technique