Three-body continuum states on a Lagrange mesh

h

esh. Thealism

tinuumy using

tion

ry of aclei,tivelynoughadrupole

Nuclear Physics A 765 (2006) 370–389

Three-body continuum states on a Lagrange mes

P. Descouvemonta,∗,1, E. Tursunova,2, D. Bayeb,a

a Physique Nucléaire Théorique et Physique Mathématique, CP229 Université Libre de Bruxelles,B-1050 Brussels, Belgium

b Physique Quantique, CP165/82 Université Libre de Bruxelles, B-1050 Brussels, Belgium

Received 19 October 2005; received in revised form 10 November 2005; accepted 14 November 2005

Available online 5 December 2005

Abstract

Three-body continuum states are investigated with the hyperspherical method on a Lagrange mR-matrix theory is used to treat the asymptotic behaviour of scattering wave functions. The formis developed for neutral as well as for charged systems. We point out some specificities of constates in the hyperspherical method. The collision matrix can be determined with a good accuracy bpropagation techniques. The method is applied to the6He (= α + n + n) and6Be (= α + p + p) systems,as well as to14Be (= 12Be+ n + n). For 6He, we essentially recover results of the literature. Applicato 14Be suggests the existence of an excited 2+ state below threshold. The calculatedB(E2) value shouldmake this state observable with Coulomb excitation experiments. 2005 Elsevier B.V. All rights reserved.

1. Introduction

Three-body systems present a large variety of interesting features [1,2]. The discovehalo structure in6He [3] triggered many experimental and theoretical works on exotic nusuch as6He,11Li or 14Be. The bound-state spectroscopy of Borromean systems is now relawell known. On the experimental side, current intensities of radioactive beams are high efor precise measurements of spectroscopic properties, such as energies, r.m.s. radii or qu

* Corresponding author.E-mail address: [email protected] (P. Descouvemont).

1 Directeur de Recherches FNRS.2 Permanent address: Institute of Nuclear Physics, 702132 Ulugbek, Tashkent, Uzbekistan.

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

P. Descouvemont et al. / Nuclear Physics A 765 (2006) 370–389 371

accurate

y sys-nsionalinto a

o many

range-ion thatximationments-mesho-body

tates.e spec-exper-been

uationsonantbe

anydy sys-ribed over

ee-bodys, how-

resentappli-

, due to

d.-ks

moments. On the theoretical side, several methods have been developed, and providesolutions of the three-body Schrödinger equation.

The hyperspherical harmonic method (HHM) is known to be well adapted to three-bodtems [4,5]. The six Jacobi coordinates are replaced by five angles, and a single-dimecoordinate, the hyperradius. The HHM transforms the three-body Schrödinger equationset of coupled differential equations depending on the hyperradius. It has been applied texotic nuclei.

Recently, we have combined the HHM with the Lagrange-mesh technique [6]. The Lagmesh method (see Ref. [7] and references therein) is an approximate variational calculatresembles a mesh calculation. The matrix elements are calculated at the Gauss approassociated with the mesh. They become very simple. In particular, the potential matrix eleare replaced by their values at the mesh points. In spite of its simplicity, the Lagrangemethod is as accurate as the corresponding variational calculation. This was shown for tw[7] as well as for three-body [6] systems.

In the present work, we extend the formalism of Ref. [6] to three-body continuum sThe information provided by continuum states is a natural complement to the bound-stattroscopy. Experimentally, three-body continuum states are investigated through breakupiments (see for example Ref. [8]). On the theoretical point of view, various methods havedeveloped. Some of them, such as the complex scaling method [9], or the analytic continin the coupling constant [10] deal with resonances only; they cannot be applied to non-restates. Other methods, such as theR-matrix theory [11] are more difficult to apply, but canused for non-resonant, as well as for resonant states.

Applications of theR-matrix method to two-body systems have been performed for myears in nuclear as well as in atomic physics. In nuclear physics, applications to three-botems are more recent [12]. TheR-matrix theory allows the use of a variational basis to descunbound states. It is based on an internal region, where the wave function is expandethe basis, and on an external region, where the asymptotic behaviour can be used. In thrsystems, the hyperspherical formalism is very efficient for bound states. For unbound stateever, it raises problems owing to the long range of the coupling potentials [12]. In theR-matrixframework this can be solved by using propagation techniques [13].

In two-body systems, the Lagrange-mesh technique associated with theR-matrix formalismhas been applied in single- [14] and multi-channel [15] calculations. The purpose of the pwork is to extend the method to three-body systems. Another development concerns thecation to charged systems. Many exotic nuclei are unbound, even in their ground statesthe Coulomb force. We show applications to theα + n + n and12Be+ n + n systems, for whichtwo-body potentials are available in the literature. The mirror systems are also investigate

In Section 2, we summarize the three-body formalism, and present theR-matrix method. Section 3 is devoted to applications to6He and14Be, with the mirror systems. Concluding remarare given in Section 4.

2. Three-body continuum states

2.1. Hamiltonian and wave functions

Let us consider three particles with mass numbersAi (in units of the nucleon massmN ), andspace coordinatesr i . A three-body Hamiltonian is given by

372 P. Descouvemont et al. / Nuclear Physics A 765 (2006) 370–389

ect

lism

he

6]

ith

nian

l

slow,alistic

H =3∑

i=1

Ti +3∑

i>j=1

Vij (rj − r i ), (1)

whereTi is the kinetic energy of nucleoni, andVij a nucleus–nucleus potential. We neglthree-body forces in this presentation.

The HHM is known to be an efficient tool to deal with three-body systems. This formais well known, and we refer to Refs. [2,5] for details. Starting from coordinatesr i , one definesthe Jacobi coordinatesxk and yk (k = 1,2,3). We adopt here the notations of Ref. [6]. Thyperradiusρ and hyperangleαk are then defined as

ρ2 = x2k + y2

k , αk = arctanyk

xk

. (2)

The hyperangleαk and the orientationsΩxkand Ωyk

provide a set of anglesΩ5k . In thisnotation the kinetic energy reads

Tρ =3∑

i=1

Ti − Tcm = − h2

2mN

(∂2

∂ρ2+ 5

ρ

∂

∂ρ− K2(Ω5k)

ρ2

). (3)

In Eq. (3),Tcm is the c.m. kinetic energy, andK2 is a five-dimensional angular momentum [1whose eigenfunctions (with eigenvaluesK(K + 4)) are given by

Yxy

KLML(Ω5) = φ

xy

K (α)[Yx (Ωx) ⊗ Yy (Ωy)

]LML,

φxy

K (α) = N xy

K (cosα)x (sinα)y P(y+ 1

2 ,x+ 12 )

n (cos2α), (4)

whereP(α,β)n (x) is a Jacobi polynomial andN xy

K a normalization factor [5] (herek is implied).In these definitions,K is the hypermomentum, (x, y ) are the orbital momenta associated w(x,y), andn is a positive integer defined by

n = (K − x − y)/2. (5)

Introducing the spin componentχSMS yields the hyperspherical function with total spinJ

YJMγK (Ω5) = [

Yxy

KL (Ω5) ⊗ χS]JM

, (6)

where indexγ stands for (x, y,L,S).A wave functionΨ JMπ , solution of the Schrödinger equation associated with Hamilto

(1), is expanded over basis functions (6) as

Ψ JMπ(ρ,Ω5) = ρ−5/2∑γK

χJπγK(ρ) YJM

γK (Ω5), (7)

whereπ = (−1)K is the parity of the three-body relative motion, andχJπγK(ρ) are hyperradia

wave functions which should be determined. Rigorously, the summation over (γK) should con-tain an infinite number of terms. In practice, this expansion is limited by a maximumK value,denoted asKmax. For weakly bound states, it is well known that the convergence is ratherand that largeKmax values must be used. Typically 100–200 terms are necessary for reKmax values.

The radial functionsχJπ (ρ) are derived from a set of coupled differential equations
γK


eus–

ad-on to,the

0)

umfe can

(10).-After

[− h2

2mN

(d2

dρ2− LK(LK + 1)

ρ2

)− E

]χJπ

γK(ρ) +∑K ′γ ′

V JπKγ,K ′γ ′(ρ)χJπ

γ ′K ′(ρ) = 0, (8)

with LK = K + 3/2. The potential terms are given by the contribution of the three nuclnucleus interactions

V JπKγ,K ′γ ′(ρ) =

3∑i=1

(V

Jπ(Ni)

Kγ,K ′γ ′(ρ) + VJπ(Ci)

Kγ,K ′γ ′(ρ)), (9)

where we have explicitly written the nuclear(N) and Coulomb (C) terms.Assuming the use of (x1,y1) for the coordinate set, the contributioni = 1 is directly deter-

mined from

VJπ(1)

Kγ,K ′γ ′(ρ) =∫

YJM∗γK (Ω5)V23

(ρ cosα√

µ23

)YJM

γ ′K ′(Ω5) dΩ5, (10)

whereµij = AiAj/(Ai + Aj). The termsi = 2,3 are computed in the same way, with anditional transformation using the Raynal–Revai coefficients [16]. Definition (10) is commthe nuclear and Coulomb contributions. Integrations overΩx andΩy are performed analyticallywhereas integration over the hyperangleα is treated numerically. For the Coulomb potential,ρ dependence is trivial; we have

3∑i=1

VJπ(Ci)

Kγ,K ′γ ′(ρ) = zJπKγ,K ′γ ′

e2

ρ, (11)

wherezJπKγ,K ′γ ′ is an effective charge, independent ofρ, and calculated numerically from Eq. (1

and from Raynal–Revai coefficients [17]. Examples of matriceszJπ are given in Ref. [17] for theα+p+p system. Knowing the analyticalρ-dependence of the potential is crucial for continustates (see below). Notice that, to derive Eq. (11), one assumes the 1/|rj − r i | dependence othe Coulomb potential. Using a point-sphere definition is straightforward, as the differencbe included in the nuclear part.

2.2. Asymptotic behaviour of the potential

For smallρ values the potential must be determined by numerical integration of Eq.However, analytical approximations can be derived for largeρ values. For the Coulomb interaction, definition (11) is always valid. Let us now consider the nuclear contribution.integration overΩx andΩy , a matrix element between basis states (4) is written as

Vxy,′

x′y

KL,K ′L′ (ρ) = δLL′δy′y

π/2∫0

φxy

K (α)VN

(ρ cosα√

µ23

)φ

′xy

K ′ (α)sin2 α cos2 α dα

= N xy

K N ′xy

K ′ δLL′δy′y

1

ρ3

ρ∫0

P(y+ 1

2 ,x+ 12 )

n

(2u2

ρ2− 1

)VN

(u√µ23

)

× P(y+ 1

2 ,′x+ 1

2 )

n′

(2u2

2− 1

)(1− u2

2

)y+ 12(

u)x+′

x

u2 du. (12)
ρ ρ ρ


totaltum

n re-f)forma-tained

ns.ae

-ust be

metryblem.

To deal with the spin, the coupling order in Eq. (6) is modified in order to introduce thespin of the interacting particlesjx = x + S. This is achieved with standard angular-momenalgebra, involving 6j coefficients. If the tensor force is not included, we also havex = ′

x .For largeρ values, and if the potential goes to zero faster than 1/u2, we can use the followingexpansions [18]

P (α,β)n (2x − 1) =

n∑m=0

c(α,β)m xm,

c(α,β)m = (−1)n+m

m!(n − m)! (β + n + 1) (α + β + n + m + 1)

(β + m + 1) (α + β + n + 1),

(1− x)α =∞∑

m=0

(α

m

)(−x)m, (13)

and we end up with the asymptotic expansion of the potential

Vxy,′

x′y

KL,K ′L′ (ρ) ≈ δLL′δy′y

1

ρx+′x+3

∞∑k=0

vk

ρ2k, (14)

where

vk = N xy

K N ′xy

K ′

∞∫0

ux+′x+2k+2V

(u√µ23

)du

×∑

m1,m2

(−1)k−m1−m2

(y + 1

2k − m1 − m2

)c(y+ 1

2 ,x+ 12 )

m1 c(y+ 1

2 ,′x+ 1

2 )m2 . (15)

Owing to the finite range of the potential, the upper limit in the integrals (12) has beeplaced by infinity. Up to a normalization factor, the contribution of eachk value is a moment othe potential. As it is well known [12], the leading term isv0/ρ

3 for x = ′x = 0. Expansion (14

is carried out for the three nucleus–nucleus potentials with additional Raynal–Revai transtions for the second and third terms. Analytic expansions of potentials (10) are finally obwith

3∑i=1

VJπ(Ni)

Kγ,K ′γ ′(ρ) ≈ 1

ρ3

∞∑k=0

vk

ρ2k, (16)

where coefficientsvk are obtained fromvk after Raynal–Revai and spin coupling transformatioLet us evaluate coefficientsvk for 6He= α +n+n, with theα–n potential taken from Kanad

et al. [19]. Coefficientsv0 to v4 are given in Table 1 forJ = 0+. We also provide the amplitudof the centrifugal term

vcent= h2

2mN

(K + 3/2)(K + 5/2), (17)

which depends onρ as 1/ρ2. It is clear from Table 1 that coefficientsvk are large and increasing with k. Integrals in (15) must be computed with a high accuracy. Special attention mpaid to partial waves involving two-body forbidden states. In this case, we use a supersymtransform of the potential [20], in order to remove forbidden states in the three-body pro


and

ularity

the

waveslarger

endousn 2.3.3.

nalation.dyan beyknown.

Table 1Coefficientsv0 to v4 in 6He for J = 0+,L = S = 0, and for typical partial waves (energies are expressed in MeVlengths in fm). The bracketed values represent the power of 10, andγ = x, y

K, γ K ′, γ ′ v0 v1 v2 v3 v4 vcent ρmax

0,0,0 0,0,0 3.40(3) −7.46(3) −2.02(4) −1.53(5) −1.78(6) 78 43704,0,0 4,0,0 1.18(3) −1.20(5) 7.31(6) −2.13(8) 2.87(9) 741 1608,0,0 8,0,0 −2.59(3) −1.19(5) 5.46(7) −6.66(9) 4.98(11) 2068 1254,2,2 4,2,2 2.61(4) −1.27(6) 5.40(7) −1.39(9) 1.81(10) 741 35208,2,2 8,2,2 5.49(4) −7.82(6) 1.06(9) −1.02(11) 6.78(12) 2068 2660

0,0,0 4,0,0 −3.41(3) 8.04(4) −1.09(6) 4.27(6) 1.43(7)

0,0,0 8,0,0 1.19(3) −1.08(5) 6.21(6) −1.75(8) 2.33(9)

0,0,0 4,2,2 9.62(3) −2.41(5) 3.47(6) −1.37(7) −4.61(7)

0,0,0 8,2,2 1.40(4) −9.90(5) 4.80(7) −1.30(9) 1.71(10)

This transformation is carried out numerically, and the resulting potential presents a singat short distances.

From Table 1, we evaluate theρ value where the nuclear part is negligible with respect tocentrifugal term. In other words,ρmax is defined as

|v0|ρ3

max= ε × vcent

ρ2max

. (18)

Values ofρmax are given in Table 1 by assumingε = 0.01. In general they are larger for loK values for two reasons: (i) the centrifugal term is of course lower, and (ii) low partial wgenerally involve forbidden states which lead to singularities in the potential, and hence tovalues ofv0.

From theρmax values displayed in Table 1, it is clear that the channel radiusa of theR matrixmust be very large. Using basis functions valid up to these distances would require trembasis sizes. This is solved by using a propagation technique which is presented in Sectio

In the analytical expansion of the potential, the maximum valuekmax is determined from therequirement

vkmax+1

a2kmax+2 vkmax

a2kmax. (19)

This yields typical valueskmax≈ 3–4, depending on the system and on the partial wave.

2.3. Three-body R matrix

2.3.1. Principle of the R matrixThe R-matrix theory is well known for many years [11]. It allows matching a variatio

function over a finite interval with the correct asymptotic solutions of the Schrödinger equWe summarize here the main ingredients of theR-matrix theory and emphasize its three-boaspects. TheR-matrix method is based on the assumption that the configuration space cdivided into two regions: an internal region, with radiusa, where the solution of (8) is given bsome variational expansion, and an external region where the exact solutions of (8) areThis is formulated as

χJπγK,int(ρ) =

N∑cJπγKi ui(ρ), (20)

i=1


gentials

e

21].

t wetion asingoing

ase is

tore

where theN functionsui(ρ) represent the variational basis, andcJπγKi are the correspondin

coefficients. In the external region, it is assumed that only the Coulomb and centrifugal potdo not vanish; we have, for an entrance channelγ ′K ′,

χJπγK,ext(ρ) = AJπ

γK

[H−

γK(kρ)δγ γ ′δKK ′ − UJπγK,γ ′K ′H+

γK(kρ)], (21)

where the amplitude is chosen as

AJπγK = iK+1(2π/k)5/2, (22)

and whereUJπ is the collision matrix, andk =√

2mNE/h2 is the wave number [12]. If the threparticles do not interact, Eq. (21) is a partial wave of a 6-dimension plane wave [16]

exp[i(kx . x + ky . y)

] = (2π)3

(kρ)2

∑xyLMLK

iKJK+2(kρ)Yxy

KLML(Ω5ρ)Yxy∗

KLML(Ω5k). (23)

For charged systems, we have

H±γK(x) = G

K+ 32(ηγK, x) ± iF

K+ 32(ηγK, x), (24)

whereGK+3/2 andFK+3/2 are the irregular and regular Coulomb functions, respectively [The Sommerfeld parametersηγK are given by

ηγK = zJπγK,γK

mNe2

h2k, (25)

wherezJπ is the effective-charge matrix (11);η therefore depends on the channel. Notice thaneglect non-diagonal terms of the Coulomb potential. This is in general a good approximathese terms are significantly smaller than diagonal terms [17]. For neutral systems, theand outgoing functionsH±

γK(x) do not depend onγ and are defined as

H±γK(x) = ±i

(πx

2

)1/2[JK+2(x) ± iYK+2(x)

], (26)

whereJn(x) andYn(x) are Bessel functions of first, and second kind, respectively. The phchosen to recover the plane wave in absence of interaction(U = I ).

For bound states (E < 0), the external wave function is written as

χJπγK,ext(ρ) = BJπ

γK W−ηγK,K+2(2κρ), (27)

whereWab(x) is a Whittaker function, andBJπγK the amplitude (κ2 = −2mNE/h2). For neutral

systems, we have

χJπγK,ext(ρ) = CJπ

γK (κρ)1/2KK+2(κρ), (28)

whereKn(x) is a modified Bessel function.

2.3.2. The Bloch–Schrödinger equationThe basic idea of theR-matrix theory is to solve Eq. (8) over the internal region. To res

the hermiticity of the kinetic energy, one solves the Bloch–Schrödinger equation(H +L(L) − E

)Ψ JMπ = L(L)Ψ JMπ, (29)


e. for

r short-ccuratehyper-and

e size40 fm)

l radiusysicso

rm.

with the Bloch operatorL(L) defined as

L(L) = h2

2mN

∑γK

∣∣YJMγK

⟩δ(ρ − a0)

1

ρ5/2

(∂

∂ρ− LγK

ρ

)ρ5/2⟨YJM

γK

∣∣, (30)

whereL is a set of arbitrary constantsLγK . In the following, we assumeLγK = 0 for positiveenergies. Formulas presented in this subsection are given for any channel radiusa0, which canbe different froma, defined in Section 2.3.1.

Let us define matrixCJπ as

CJπγKi,γ ′K ′i′ = ⟨

uiYJMγK

∣∣H +L(L) − E∣∣ui′YJM

γ ′K ′⟩I, (31)

where subscriptI means that the matrix element is evaluated in the internal region only, i.ρ a0. Using the partial-wave expansion (7) and the continuity of the wave function atρ = a0,we obtain theR-matrix ata0 from

RJπγK,γ ′K ′(a0) = h2

2mNa0

∑i,i′

ui(a0)(CJπ

)−1γKi,γ ′K ′i′ui′(a0). (32)

2.3.3. R-matrix propagation and collision matrixAs shown in Section 2.2, the nuclear potential extends to very large distances, even fo

range nucleus–nucleus interactions. In other words, the asymptotic behaviour (21) is not abelow distances which may be as large as 1000 fm or more. This is a drawback of thespherical method, where even for largeρ values, two particles can still be close to each othercontribute to the three-body matrix elements.

It is clear that using basis functions valid up to distances of 1000 fm is not realistic, as thof the basis would be huge. On the other hand, using a low channel radius (typically 30–would keep the basis size in reasonable limits, but would not satisfy the key point of theR-matrixtheory, namely that the wave function has reached its asymptotic behaviour at the channea0. This problem can be solved with propagation techniques, well known in atomic ph[13]. The idea is to usea0 as a starting point for theR matrix; its value is small enough tallow reasonable basis sizes. TheR matrix is then propagated froma0 to a, where the Coulombasymptotic behaviour (21) is valid. Betweena0 anda, the wave functionsχJπ (ρ) are still givenby Eq. (8), but with the potential replaced by its (analytical) asymptotic expansion.

More precisely, the internal wave functions in the different intervals are given by

χJπγK,int(ρ) =

∑Ni=1 cJπ

γKi ui(ρ) for ρ a0,

χJπγK(ρ) for a0 ρ a,

(33)

whereχγK(ρ) are solutions of Eq. (8) with the analytical expansion (16) of the potential teTheR matrix is first computed ata0 with Eq. (32) (typical values area0 ≈ 20–40 fm). Then

we considerN0 sets of initial conditions forχ(ρ), whereN0 is the number ofγK values (fromnow on we drop theJπ index for clarity). We combine these sets as matrixχ0(ρ), and choose

χ0(a0) = I , (34)

whereI is the unit matrix.According to the definition of theR matrix [11], we immediately find the derivative ata0

χ ′0(a0) = 1

R−1(a0)χ0(a0) = 1R−1(a0). (35)

a0 a0


(16)

n, ands not relylations

3].horter(typicalcase,ssumes

msn, the

botht

n the

Knowing functionsχ0γK and their derivatives ata0, they are then propagated untila by usingthe Numerov algorithm [22], well adapted to the Schrödinger equation. The analytical formof the potential is used, with a summation limited to a fewk values. TheR matrix ata is thendetermined by using Eq. (35) withχ0(a) andχ ′

0(a). We have

R(a) = 1

aχ0(a)

(χ ′

0(a))−1

. (36)

Notice that the propagatedR matrix (36) does not depend on the choice ofχ0(a0). In Ref. [13],the propagation is performed through the Green function defined in the intermediate regioexpanded over a basis. The method presented here uses the Numerov algorithm, and doeon the choice of a basis. The analytical form of the potential in this region makes calcufast and accurate.

Finally the collision matrix is obtained from theR matrix at the channel radiusa with

UJπ = (ZJπ

)−1ZJπ , (37)

and

ZJπγK,γ ′K ′ = H−

γK(ka)δγ γ ′δKK ′ − ka(H−

γK(ka))′RJπ

γK,γ ′K ′(a), (38)

where the derivation is performed with respect toka.Lower values of the channel radiusa can be used by employing the Gailitis method [2

In this method the asymptotic forms (24) are generalized with the aim of using them at sdistances. This means that the propagation should be performed in a more limited rangevalues fora area ∼ 200–400 fm). However this does not avoid propagation which, in anyis very fast. In addition, the Gailitis method cannot be applied to charged systems, as it afrom the very beginning that the coupling potentials decrease faster than 1/ρ.

The extension of theR-matrix formalism to bound states is well known for two-body syste[24]. Basically, theLγK constants are defined so as to cancel the r.h.s. of Eq. (29). Theproblem is reduced to a matrix diagonalization with iteration on the energy [24,25].

2.3.4. Wave functionsOnce the collision matrix is known, the internal wave function (33) can be determined in

intervals. Although the choice ofχ0(a0) is arbitrary, functionsχ(ρ) entering Eq. (33) do nodepend on that choice. In the intermediate regiona0 ρ a, functionsχ(ρ) and χ0(ρ) arerelated to each other by a linear transformation

χ(ρ) = χ0(ρ)M. (39)

Matrix M is deduced by using the asymptotic behaviour (21) atρ = a,

χ(a) = χ0(a)M = χext(a), (40)

whereχext(a) is the matrix involving all entrance channels [see Eq. (21)]. It depends ocollision matrix.

CoefficientscJπγKi defining the internal wave function in the intervalρ a0 are finally obtained

by

cJπγKi = h2

2mN

∑γ ′K ′i′

(C−1)Jπ

γKi,γ ′K ′i′

(dχJπ

γ ′K ′

dρ

)ρ=a0

ui′(a0). (41)


thodotice

tion ofnergy

o-h

stems.n thatcan be

odyd a su-aches

stems,xpan-ently,nucleus

showhases.

is int be as-hases.

2.3.5. The Lagrange-mesh methodUp to now, the basis functionsui(ρ) are not specified. We use here the Lagrange-mesh me

which has been proved to be quite efficient in two-body [26] and three-body [6] systems. Nhowever that its application to three-body continuum states is new.

When dealing with a finite interval, theN basis functionsui(ρ) are defined as [14]

ui(ρ) = (−1)N−i

(1− xi

a0xi

)1/2ρPN(2ρ/a0 − 1)

ρ − a0xi

, (42)

where thexi are the zeros of a shifted Legendre polynomial given by

PN(2xi − 1) = 0. (43)

The basis functions satisfy the Lagrange condition

ui(a0xj ) = (a0λi)−1/2δij , (44)

where theλi are the weights of the Gauss–Legendre quadrature corresponding to the[0,1] inter-val, i.e. half of the weights corresponding to the traditional interval[−1,1].

The main advantage of the Lagrange-mesh technique is to strongly simplify the calculamatrix elements (31) if the Gauss approximation is used. Matrix elements of the kinetic e(T + L) are obtained analytically [14]. Integration overρ provides matrix elements of the ptential by a simple evaluation of the potential atρ = a0xi . The potential matrix is diagonal witrespect toi andi′.

In Ref. [6], we applied the Lagrange-mesh technique to bound states of three-body syAs the natural interval ranges from zero to infinity, we used a Laguerre mesh. It was showthe Gauss quadrature is quite accurate for the matrix elements, and that computer timesstrongly reduced.

3. Applications

3.1. Conditions of the calculations

Here we apply the method to the6He and14Be nuclei. Theα–n and 12Be–n interactionsare chosen as local potentials. They contain Pauli forbidden states (one in = 0 for α − n, andone in = 0,1 for 12Be–n) which should be removed for a correct description of three-bstates [6,12]. For bound states, two methods are available: the use of a projector [27], anpersymmetric transformation of the nucleus–nucleus potential [20]. Although both approprovide different wave functions, spectroscopic properties are similar [6]. For unbound syit turns out that the projector technique is quite difficult to apply with a good accuracy. Esions similar to Eq. (16) for the projection operator provide non-local potentials. Consequall applications presented here are obtained with supersymmetric partners of the nucleus–potentials.

As collision matrices can be quite large, it is impossible to analyze all elements. Tothe essential information derived from the collision matrix, we rather present some eigenpThose presenting the largest variation in the considered energy range are shown.

Analyzing the collision matrix in terms of eigenphases raises two problems. First, itgeneral not obvious to link the eigenphases at different energies. As eigenphases cannosociated with given quantum numbers, there is no direct way to draw continuous eigenp


giveneen the

As ma-appearontri-nolli-

produce

lect to

somee thehtion.lr of basisnd thise of theyently

onalith

wOne ofour isig. 2,nphases

.conver-

he

The procedure can be strongly improved by analyzing the eigenvectors. Starting from aenergy, eigenphases for the next energy are chosen by minimizing the differences betwcorresponding eigenvectors.

A second problem associated with eigenphases arises from the Coulomb interaction.trix elements of the Coulomb force are not diagonal, the corresponding phase shifts do notin a simple way, as in two-body collisions. Consequently, in order to extract the nuclear cbution UN from the total collision matrixU , we perform two calculations: a full calculatioproviding U , and a calculation without the nuclear contribution, providing the Coulomb csion matrixUC . Then we define the nuclear collision matrixUN by

U = U1/2C UNU

1/2C . (45)

AsU andUC are symmetric and unitary, the same properties hold forUN . Examples of Coulombphase shifts will be given in the next sections.

3.2. Application to 6He and 6Be

The conditions of the calculation are those of Ref. [6]. Theα–n potentialVα–n has beenderived by Kanada et al. [19]. It contains spin–orbit and parity terms. Then–n potential is theMinnesota interaction [28]. As bare nucleus–nucleus potentials cannot be expected to rethe 6He ground-state energy with a high accuracy, we renormalizeVα–n by a factorλ = 1.051(note that this value was misprinted in Ref. [6]). This value reproduces the6He experimentaenergy−0.97 MeV and provides 2.44 fm for the r.m.s. radius. The convergence with respKmax and to the Lagrange-mesh parameters has already been discussed in Ref. [6].

Let us first illustrate the importance of the propagation technique. In Fig. 1, we plotelements of theJ = 0+ collision matrix under different conditions. In each case, we comparphase shifts for two channel radii:a0 = 20 fm anda0 = 30 fm. The calculation is performed witand without propagation. ForK = 0, reasonable values can be obtained without propagaHowever, for largerK values (K = 8 is displayed withx = y = 0 andx = y = 4), the channeradius should be quite large to reach convergence. To keep the same accuracy, the numbefunctions should be increased. However, one basis function per fm is a good estimate, aleads to unrealistically large basis sizes. This convergence problem is due to the long rangpotential. The propagation technique (performed here up toa = 250 fm) allows us to get a verhigh stability (better than 0.1 at all energies) even for rather small channel radii. Consequcalculations with highK values are still feasible.

To illustrate the diagonalization of the collision matrix, we compare in Fig. 2 the diagphase shifts with the corresponding eigenphases. We have selected a particular case, wJ =2+, andKmax = 2. With these conditions the collision matrix is 4× 4, and presents a narroresonance near 2 MeV. In the upper part of Fig. 2, we plot the diagonal phase shifts.them presents a 180 jump, characteristical of narrow resonances. This resonant behavialso observable in two other partial waves. After diagonalization of the collision matrix (Flower part) the resonant behaviour shows up in one eigenphase only. The three other eigesmoothly depend on energy.

The convergence with respect toKmax is illustrated in Fig. 3 with theJ = 0+ eigenphasesIt turns out that, at low energies, high hypermomenta are necessary to achieve a precisegence. However, above 4 MeV,Kmax= 20 is sufficient to obtain an accuracy of 2.

Fig. 4 gives the eigenphases forJ = 0+,1−,2+ in 6He and6Be (Kmax is taken as 24, 19and 16, respectively). As expected, the 2+ phase shift of6He presents a narrow resonance. T


g

[29] isd, for

nilin et

s

Fig. 1. α + n + n phase shifts (J = 0+) for channel radiia0 = 20 fm (N = 20) anda0 = 30 fm (N = 30), and fordifferent partial waves. Solid lines are obtained with propagation up toa = 250 fm of theR matrix (curves correspondinto differenta0 are undistinguishable), and dashed lines without propagation.

theoretical energy (about 0.2 MeV) is however underestimated as the experimental valueE = 0.82 MeV. In order to provide meaningful properties for this state, we have readjusteJ = 2+, the scaling factor toλ = 1.020, which provides the correct energy. The 0+ and 1− phaseshifts show broad structures near 1.5 MeV. Similar phase shifts have been obtained by Daal. [30,31] and by Thompson et al. [12] with other potentials.

In 6Be, the ground state is found atE = 1.26 MeV with a widthΓ = 65 keV. These valueare in reasonable agreement with experiment [29] (E = 1.37 MeV,Γ = 92± 6 keV), the width


ithtomb

Fig. 2. Diagonal phase shifts (upper panel) and eigenphases (lower panel) for theα +n+n system (J = 2+,Kmax= 2).

Fig. 3. Energy dependence ofα + n + n eigenphases (J = 0+) for differentKmax values.

being underestimated by the model due to the lower energy. Experimentally, a 2+ state is knownnearE = 3.0 MeV with a width ofΓ = 1.16± 0.06 MeV. These properties are consistent wthe theoretical 2+ eigenphase, which presents a broad structure nearE ≈ 4 MeV. The largesCoulomb eigenphases (J = 0+) are shown as dotted lines in Fig. 4. As expected, the Coul


st

d evenare not

,

he E2

ticd iniptions

Fig. 4. Eigenphases of6He and6Be for differentJ values (solid lines). For6Be, dotted lines represent the largeCoulomb eigenphases forJ = 0+.

interaction plays a dominant role at low energies, but it cannot be completely neglectenear 10 MeV. Coulomb eigenphases for other spin values are very similar and thereforepresented. Energies and widths are given in Table 2.

In Table 2, we also present the E2 transition probability in6He. For the narrow 2+ resonancewe use the bound-state approximation. Without effective charge, theB(E2) value for the 0+ →2+ transition is underestimated with respect to the experimental value [8]. However, tmatrix element is very sensitive to the effective charge. A small correction (δe = 0.05e) providesaB(E2) within the experimental error bars.

3.3. Application to 14Be

As shown in previous works [34–36], a12Be+n+n three-body model can provide a realisdescription of14Be. The spectroscopy of the14Be ground state has already been investigatenon-microscopic [34–37] and microscopic [38] models. Here we extend three-body descrto 14Be excited and continuum states.


ingll

po-

ta

sncom-factor

n-

,yvidell

the

eter-

Table 26He and6Be properties. Unless specified, experimental data are taken from Ref. [29]

6He 6Be

present exp. present exp.

E(0+) (MeV) −0.97 −0.97 1.26 1.37Γ (0+) (keV) 65 92± 6E(2+) (MeV) 0.8 0.82 ≈ 4.0 3.04Γ (2+) (MeV) 0.04 0.113± 0.020 ≈ 1.0 1.16± 0.06√

〈r2〉 (fm) 2.44 2.33± 0.04a

2.57± 0.10b

2.45± 0.10c

B(E2,0+ → 2+) (e2 fm4) 1.23 (δe = 0) 3.2± 0.6d

2.69 (δe = 0.05e)

a Ref. [3], b Ref. [32], c Ref. [33], d Ref. [8].

The13Be ground state is expected to be a virtuals wave, with a large and negative scatterlength (as < −10 fm) [39]. In addition, the existence of a 5/2+ d-state near 2 MeV is weestablished. These properties can be reproduced by a12Be–n potential

V (r) = − V0 + Vs · s1+ exp((r − r0)/a)

, (46)

where is the relative angular momentum ands the neutron spin. In Eq. (46),r0 = 2.908 fm,a = 0.67 fm, V0 = 43 MeV, Vs = 6 MeV. The range and diffuseness of the Woods–Saxontential are taken from Ref. [36]. The amplitudesV0 andVs provideE(5/2+) = 2.1 MeV, andas = −47 fm, which are consistent with the data. For then–n potential, we use the Minnesointeraction, as for the6He study.

With these potentials, the14Be ground state is found atE = −0.16 MeV, which representan underbinding with respect to experiment (−1.34± 0.11 MeV [40]). This calculation has beeperformed withKmax = 24, which ensures the convergence. The underbinding problem ismon to all three-body approaches, and can be solved in two ways. (i) A renormalizationλ = 1.08 provides a ground-state energy at−1.30 MeV, i.e. within the experimental uncertaities. This procedure leads to a slightly bound13Be ground state, which might influence the14Beproperties. (ii) A three-body phenomenological termV (123), taken as in Ref. [12], i.e.,

V(123)Kγ,K ′γ ′(ρ) = −δKK ′δγ γ ′ V3/

[1+ (ρ/ρ3)

2], (47)

reproduces the experimental energy with an amplitudeV3 = 4.7 MeV (according to Ref. [12]we takeρ3 = 5 fm). This potential is diagonal in(K,γ ), and is simply added to the two-bodterm [see Eqs. (10), (11)]. In6He, it was shown that both readjustments of the interaction prosimilar results [6]. However the renormalization factor is larger for14Be, and both methods wibe considered in the following.

The convergence with respect toKmax is illustrated in Fig. 5. ForJ = 0+, the calculationshave been done withKmax up to 24. The energies obtained with renormalization or withthree-body potential are very similar. This confirms the conclusion drawn for the6He nucleus[6].

Spectroscopic properties of14Be are given in Tables 3 and 4. The r.m.s. radii have been dmined with 2.57 fm as12Be radius. For the ground state, we have

√〈r2〉 = 3.10 fm or 3.14 fm,


a

ble

al. 5w

re of the

havein

ifferent

Fig. 5. Energy of the14Be 0+ (squares) and 2+ (circles) states as a function ofKmax. Full symbols correspond torenormalized12Be–n potential, and open symbols to a phenomenological three-body term (see text).

Table 3Properties of the14Be 0+ and 2+ states.λ is the renormalization factor of the12Be–n potentialandV3 is this amplitude of the three-body potential

λ = 1.08,V3 = 0 λ = 1, V3 = 4.7

0+ E (MeV) −1.34 −1.34√〈r2〉 (fm) 3.10 3.14

PS=1 0.046 0.033

2+ E (MeV) −0.15 −0.03√〈r2〉 (fm) 2.99 3.04

PS=1 0.192 0.165

0+ → 2+ B(E2) (e2 fm4) 0.48 (δe = 0) 0.64 (δe = 0)

3.18 (δe = 0.05e) 4.05 (δe = 0.05e)

in nice agreement with experiment (3.16± 0.38 fm, see Ref. [41]). In all cases, theS = 1 com-ponent (denoted asPS=1) is small (< 5%). The decomposition in shell-model orbitals (see Ta4) shows that the 0+ state is essentially (≈ 70%) (2s1/2)

2, with small (2d3/2)2 and (2d5/2)

2

admixtures.RegardingJ = 2+, we have considered values up toKmax = 16, where the number of parti

waves is 172. Going beyondKmax = 16 would require too large computer memories. Figshows the energy convergence with respect toKmax. For both potentials, the energy is belothreshold, and the r.m.s. radius is close to 3 fm. A partial-wave analysis provides 19% ofS = 1admixture, a value much larger than in the ground state. Table 4 suggests that the structu2+ state is spread over many components. The(s1/2d5/2) component is dominant (≈ 23%) butother(sd) and(pf ) orbitals also play a role.

E2 transition probabilities are also given in Table 3. Without effective charge, weB(E2,0+ → 2+) = 0.48 and 0.64e2 fm4, which is lower than for the corresponding transition6He. However, the amplitudes of the proton and neutron E2 operators being even more din 14Be than in6He, theB(E2) values strongly depend on the effective charge. Forδe = 0.05e,


of

ties

tr”es

. As

Table 4Components (in %) in14Be wave functions

0+ (λ = 1.08,V3 = 0) 0+ (λ = 1,V3 = 4.7)

(p3/2)2 2.0 2.2(p1/2)2 1.0 1.1(s1/2)2 70.4 73.1(d5/2)2 14.6 13.0(d3/2)2 11.2 9.8

2+ (λ = 1.08,V3 = 0) 2+ (λ = 1,V3 = 4.7)

p3/2p3/2 7.7 8.2p3/2f7/2 9.6 9.4p1/2p3/2 18.0 18.6p1/2f5/2 5.8 5.7s1/2d5/2 23.2 23.5s1/2d3/2 19.2 18.5d5/2d5/2 5.6 5.3d3/2d5/2 4.0 3.6d3/2d3/2 3.0 2.9

Fig. 6. Left panel: radial functionsχ(ρ) for the 0+ state. The curves are labeled byx, y , n. Right panel: probabilityP(rnn, rBe–nn), deduced from Eq. (48) withrnn = √

2x andrBe–nn = √6/7y. Contour levels are plotted by steps

0.005.

we find B(E2) = 3.18 or 4.05e2 fm4 according to the potential. Such transition probabilishould be measurable through Coulomb excitation experiments.

In Figs. 6–7, we present the 0+ and 2+ radial wave functions and probabilitiesP(x, y) definedas

P Jπ(x, y) =∫

dΩx dΩy x2y2∣∣Ψ JMπ(x,y)

∣∣2. (48)

The dominantS = 0 components are plotted. The 0+ probability shows two well distincmaxima, which resemble the maxima found in6He, corresponding to “dineutron” and “cigaconfigurations. Partial wavesχJπ

γK(ρ) have maxima forρ > 5 fm. This corresponds to distanc

larger than in6He [6] where the maxima of the main components are located near 4 fmexpected, the 2+ probability is similar to the 0+ probability, with two maxima.


asthe 2

A very

s. As

paredlinghyper-

eiques.is

opertiesitha

ions.mody

to the

f partialhbachuse a

entials

Fig. 7. See Fig. 6 for the 2+ state.

Three-body eigenphases are displayed in Fig. 8. As for the14Be spectroscopy the use ofthree-body potential does not qualitatively change the phase shifts. The 1− phase shift presenttwo jumps but they cannot be directly assigned to physical resonances. On the contrary,+phase shift shows a narrow resonance near 2 MeV. For the sake of completeness,12O + p + p

mirror phase shifts are also shown in Fig. 8. As expected, no narrow structure is found.broad 0+ resonance shows up near 8 MeV, and should correspond to the14Ne ground state.

4. Conclusion

In this work, we have extended the three-body formalism of Ref. [6] to unbound statefor two-body systems, the Lagrange-mesh technique, associated with theR-matrix method, pro-vides an efficient and accurate way to derive collision matrices and wave functions. Comwith two-body systems, three-bodyR-matrix approaches are more tedious, owing to the couppotentials which extend to very large distances. This behaviour is inherent to the use ofspherical coordinates which provide three-body potentials behaving as 1/ρ3, even for short-rangtwo-body interactions. This problem can be efficiently solved by using propagation technHere, we propagate the wave function and theR matrix by using the Numerov algorithm. Thformalism has been extended to charged systems.

The 6He system has essentially been used as a test of the method, as most of its prare available in the literature. TheB(E2,0+ → 2+) experimental value can be reproduced wa small effective chargeδe = 0.05e. We have determinedα + p + p phase shifts, and foundgood agreement with experiment for the6Be ground-state properties.

Application to three-body12Be+n+n states is new, and has been developed in two directThe bound-state description of14Be provides evidence for a 2+ bound state, as expected frothe shell model. The study of the12Be+ n + n system has been complemented by three-bphase shifts, which suggest the existence of a second narrow 2+ resonance atEx ≈ 3.4 MeV.

A limitation of the method is the slow convergence of the phase shifts with respectmaximum hypermomentumKmax. To achieve a full convergence, values up toKmax = 20 ormore are necessary. This problem is even stronger for high spins, where the number owaves increases rapidly. A possible solution to this problem would be to apply the Fesreduction method [42] to scattering states. Another possible development would be toprojection technique to remove Pauli forbidden states [10]. In that case, asymptotic pot(15) are non local, which makes the calculation still heavier.


poten-

exoticditions.sting

prob-tractionffairs.

Fig. 8. Three-body12Be–n–n and12O–p–p eigenphases. Solid lines correspond to a renormalized core–nucleontial, and dotted lines to a phenomenological three-body term.

The present model offers an efficient way to investigate bound and unbound states. Innuclei, most low-lying states are unbound, and a rigorous analysis requires scattering conThe inclusion of the Coulomb interaction still extends the application field, and is intereeven for non-exotic nuclei. In this context, an accurate analysis of unboundα + α + α statesseems desirable in view of its strong interest in the triple-α reaction rate [43].

Acknowledgements

We are grateful to Prof. F. Arickx for useful discussions about the three-body Coulomblem. This text presents research results of the Belgian program P5/07 on interuniversity atpoles initiated by the Belgian-state Federal Services for Scientific, Technical and Cultural AOne of the authors (E.M.T.) is supported by the SSTC.


1.yashi,

18.

Nucl.

. Sato,

References

[1] B. Jonson, Phys. Rep. 389 (2004) 1.[2] M.V. Zhukov, B.V. Danilin, D.V. Fedorov, J.M. Bang, I.J. Thompson, J.S. Vaagen, Phys. Rep. 231 (1993) 15[3] I. Tanihata, H. Hamagaki, O. Hashimoto, Y. Shida, N. Yoshikawa, K. Sugimoto, O. Yamakawa, T. Koba

N. Takahashi, Phys. Rev. Lett. 55 (1985) 2676.[4] P.M. Morse, H. Feshbach, Methods in Theoretical Physics, vol. II, McGraw–Hill, New York, 1953.[5] C.D. Lin, Phys. Rep. 257 (1995) 1.[6] P. Descouvemont, C. Daniel, D. Baye, Phys. Rev. C 67 (2003) 044309.[7] D. Baye, M. Hesse, M. Vincke, Phys. Rev. E 65 (2002) 026701.[8] T. Aumann, et al., Phys. Rev. C 59 (1999) 1252.[9] Y.K. Ho, Phys. Rep. 99 (1983) 1.

[10] V.I. Kukulin, V.M. Krasnopol’sky, J. Phys. A 10 (1977) 33.[11] A.M. Lane, R.G. Thomas, Rev. Mod. Phys. 30 (1958) 257.[12] I.J. Thompson, B.V. Danilin, V.D. Efros, J.S. Vaagen, J.M. Bang, M.V. Zhukov, Phys. Rev. C 61 (2000) 0243[13] V.M. Burke, C.J. Noble, Comput. Phys. Commun. 85 (1995) 471.[14] D. Baye, M. Hesse, J.-M. Sparenberg, M. Vincke, J. Phys. B 31 (1998) 3439.[15] M. Hesse, J.-M. Sparenberg, F. Van Raemdonck, D. Baye, Nucl. Phys. A 640 (1998) 37.[16] J. Raynal, J. Revai, Nuovo Cimento A 39 (1970) 612.[17] V. Vasilevsky, A.V. Nesterov, F. Arickx, J. Broeckhove, Phys. Rev. C 63 (2001) 034606.[18] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover, London, 1972.[19] H. Kanada, T. Kaneko, S. Nagata, M. Nomoto, Prog. Theor. Phys. 61 (1979) 1327.[20] D. Baye, Phys. Rev. Lett. 58 (1987) 2738.[21] I.J. Thompson, A.R. Barnett, J. Comput. Phys. 64 (1986) 490.[22] J. Raynal, in: Computing as a Language of Physics, Trieste 1971, IAEA, Vienna, 1972, p. 281.[23] M. Gailitis, J. Phys. B 9 (1976) 843.[24] D. Baye, P. Descouvemont, Nucl. Phys. A 407 (1983) 77.[25] P. Descouvemont, M. Vincke, Phys. Rev. A 42 (1990) 3835.[26] M. Hesse, J. Roland, D. Baye, Nucl. Phys. A 709 (2002) 184.[27] V.I. Kukulin, V.N. Pomerantsev, Ann. Phys. 111 (1978) 330.[28] D.R. Thompson, M. LeMere, Y.C. Tang, Nucl. Phys. A 286 (1977) 53.[29] D.R. Tilley, C.M. Cheves, J.L. Godwin, G.M. Hale, H.M. Hofmann, J.H. Kelley, C.G. Sheua, H.R. Weller,

Phys. A 708 (2002) 3.[30] B.V. Danilin, I.J. Thompson, J.S. Vaagen, M.V. Zhukov, Nucl. Phys. A 632 (1998) 383.[31] B.V. Danilin, T. Rogde, J.S. Vaagen, I.J. Thompson, M.V. Zhukov, Phys. Rev. C 69 (2004) 024609.[32] L.V. Chulkov, et al., Europhys. Lett. 8 (1989) 245.[33] G.D. Alkhazov, A.V. Dobrovolsky, A.A. Lobodenko, Nucl. Phys. A 734 (2004) 361.[34] D. Baye, Nucl. Phys. A 627 (1997) 305.[35] A. Adahchour, D. Baye, P. Descouvemont, Phys. Lett. B 356 (1995) 445.[36] I.J. Thompson, M.V. Zhukov, Phys. Rev. C 53 (1996) 708.[37] T. Tarutina, I.J. Thompson, J.A. Tostevin, Nucl. Phys. A 733 (2004) 53.[38] P. Descouvemont, Phys. Rev. C 52 (1995) 704.[39] M. Thoennessen, S. Yokoyama, P.G. Hansen, Phys. Rev. C 63 (2001) 014308.[40] G. Audi, A.H. Wapstra, Nucl. Phys. A 565 (1993) 1.[41] I. Tanihata, T. Kobayashi, O. Yamakawa, S. Shimoura, K. Ekuni, K. Sugimoto, N. Takahashi, T. Shimoda, H

Phys. Lett. B 206 (1988) 592.[42] H. Feshbach, Ann. Phys. 19 (1962) 287.[43] H. Fynbo, et al., Nature 433 (2005) 136.

Documents

Three-body continuum states on a Lagrange mesh