Operation of Faddeev-Kernel in Configuration Space

0.34 MB·English

by  S. Ishikawa

Few-Body Systems 0, 1–20 (2008) Few-
Body
Systems
(cid:13)c by Springer-Verlag 2008
PrintedinAustria
Operation of Faddeev-Kernel in
7 Configuration Space
0
0
2
S. Ishikawa∗
n
a
J
Department of Physics, Science Research Center, Hosei University,2-17-1 Fujimi, Chiyoda,

6 Tokyo102-8160, Japan
1

1
v
4 Abstract. We present a practical method to solve Faddeev three-body equa-
4 tions at energies above three-body breakup threshold as integral equations in
0
coordinate space. This is an extension of previously used method for bound
1
0 states and scattering states below three-body breakup threshold energy. We
7 show that breakup components in three-body reactions produce long-range
0
effects on Faddeev integral kernels in coordinate space, and propose numeri-
/
h
calprocedurestotreat theseeffects. Usingthesetechniques, wesolve Faddeev
t
- equations for neutron-deuteron scattering to compare with benchmark solu-
l
c tions.
u
n
:
v
i 1 Introduction
X
r So far, a number of numerical methods to solve Faddeev three-body equations
a
for energies above three-body breakup threshold have been developed and then
applied to a system of nucleon-deuteron, which is considered as one of the most
basic quantum three-body systems [1, 2]. These methods are classified into two
groups: either to solve coupled integral equations for scattering amplitudes in
momentum space or to solve coupled partial differential equations for wave func-
tions in coordinate space.
In this paper, we will present a different approach for three-body scattering
problem above the breakup threshold energy, in which we solve the Faddeev
integral equations for wave functions in coordinate space. This approach has
been successfully applied to calculations of the three-nucleon bound states [3, 4]
andlow-energythree-nucleonscatteringbelowthebreakupthresholdenergywith
inclusionofthree-nucleonforces[4]andthelong-rangeCoulombinteraction[5,6].
Integral equations for scattering problems are generally written in the form
ofinhomogeneouslinearequations.Inthepreviousworks,weappliedaniterative
method called the method of continued fraction (MCF) to solve such equations,
whosedetails aregiven inRefs.[3,7]andreferencestherein.Abasicprocedureof
∗E-mail address: [email protected]
2 Faddeev-Kernelin Configuration Space
the algorithm in the MCF is to operate the integral kernel to a function made in
aprecedingstep,asarethosein mostiterative methods.Itis thusessential toes-
tablish precise operations of integral kernels for solving the equations accurately,
which is main subject of this paper.
The existence of three-body breakup channels causes some difficulties in
three-body calculations. In the momentum space approach, for example, the
effects appear as logarithmic singularities and discontinuities by a step function
intheintegralkernelsoftheequationssothatweneedtoperformtheintegration
very carefully [8]. In the differential equation approach, due to the breakup ef-
fectsoneneedstosetboundaryconditionsatverylongdistance,theorderoftens
or hundreds times larger than the range of interaction potentials [9, 10, 11, 12].
Since we treat the wave functions as solutions of the Faddeev integral equations,
the long-range behavior should appear in the integral kernel. In the present pa-
per, we will describe how this behavior appears in our kernel, and how to treat
it.
Basic notations and steps of the kernel operation in detail are explained in
Sec. 2 for a simple three-body system. In Sec. 3, we show numerical examples of
the kernel operation to a model function emphasizing some techniques to treat
breakup effects, and then compare our calculations with benchmark tests [1, 2].
Finally, we give a summary in Sec. 4.
2 Formulation
2.1 Notations
Let us consider a system of three identical particles (nucleons) 1, 2, and 3. We
use sets of Jacobi coordinates {x ,y } defined as
i i
x = r −r
i j k , (2.1)
y = r − 1(r +r )
( i i 2 j k
where (i,j,k) denotes (1,2,3) or its cyclic permutations and r is the position
i
vector of the nucleon i. For simplicity, we assume that the nucleons j and k
interact via a short range pair wise potential V = V(x ), where x = |x |, and
i i i i
that the potential supports a s-wave bound state (the deuteron) of energy E ,
d
whose radial part wave function is denoted by φd(x ).
i
We are going to obtain a wave function Ψ corresponding to a scattering pro-
cessinitiatedbyastatewithanucleonandadeuteronhavingrelativemomentum
p . Faddeev equations for the process in the form of integral equations are
0
Φ = Ξ +G V (Φ +Φ )= Ξ +G V PˆΦ , (2.2)
i i i i j k i i i i
where Φ ’s are Faddeev components to make Ψ as
i
Ψ = Φ +Φ +Φ , (2.3)
1 2 3
Ξ is an initial state consisting of the deuteron for the pair (j,k) and incoming
i
freenucleoni,andG isathree-bodychannelGreen’soperatorwiththeoutgoing
i
S.Ishikawa 3
boundary condition,
1
G = . (2.4)
i E+iε+ ¯h2∇2 + 3¯h2∇2 −V
m xi 4m yi i
The total energy in the three-body center of mass (c.m.) frame E is given as
3h¯2 3h¯2
E = p2+E = p2−|E |, (2.5)
4m 0 d 4m 0 d
where m denotes the nucleon mass. The operator Pˆ represents permutations of
the particle numbers,
PˆΦ = Φ +Φ . (2.6)
i j k
A partial wave decomposition is performed by introducing an angular func-
tion Y (xˆ ,yˆ),
α i i
Y (xˆ ,yˆ) = [Y (xˆ )⊗Y (yˆ )]J0 , (2.7)
α i i L i ℓ i M0
where L denotes the relative orbital angular momentum of the pair (j,k); ℓ the
orbital angular momentum of the spectator i with respect to the c.m. of the pair
(j,k); J the total angular momentum of the three-body system (J = L+ℓ);
0 0
M the third component of J . The set of the quantum numbers (L,ℓ,J ,M )
0 0 0 0
are represented by the index α. Furthermore, we use an index α to denote an
0
initial partial wave state specifically with L = 0.
2.2 Kernel Operation
Inthissubsection,wedescribehowtohandletheoperationoftheFaddeevkernel
GVPˆ on a given function Ξ,
hx,y|Ξi = Y (xˆ,yˆ)ξ (x,y) (2.8)
α α
α
X
to produce a new function Φ,
hx,y|Φi = hx,y|GVPˆ|Ξi
= Y (xˆ,yˆ)φ (x,y), (2.9)
α α
α
X
where we have dropped the particle number indices (i,j,k) for simplicity.
The kernel operation starts with the permutation operator Pˆ to define a
function χ (x,y),
α
χ (x,y) = (Y |Pˆ|Ξi. (2.10)
α α
Inthecaseofidenticalparticles,Pˆ isnothingbutacoordinateexchangeoperator,
whose operations are summarized in Appendix A:.
Next step is the operation of the Green’s operator G. In the case of the
scattering problem, where E > 0, the Green’s operator G possesses a pole cor-
responding to the deuteron bound state. In order to treat this pole, we apply a
standard subtraction method, in which we insert a trivial identity,
1= |Y φd)(φdY |+ 1− |Y φd)(φdY | , (2.11)
α0 α0 α0 α0
Xα0 " Xα0 #
4 Faddeev-Kernelin Configuration Space
between G and V in Eq. (2.9). This procedure extracts an elastic contribution
of the Green’s operator [13] and leads to an expression,
φ (x,y) = δ φd(x)F(e)(y)+φ(b,c)(x,y). (2.12)
α α,α0 α
Here, F(e)(y) represents an elastic component in the scattering given by
∞
F(e)(y)= y′2dy′G˘ (y,y′)ω(e)(y′), (2.13)
0,ℓ0
Z0
whereG˘ (y,y′) is a partial wave component of the free Green’s operator for the
0,ℓ
outgoing particle,
1
G˘ (y,y′) ≡ y y′
0,ℓ (cid:12)3¯h2p2+iε−T (y)(cid:12) !
(cid:12) 4m 0 ℓ (cid:12)
4(cid:12)m (cid:12)
(cid:12) (+) (cid:12)
= −3h¯(cid:12)2p0hℓ (p0y>)jℓ(p(cid:12)0y<) (2.14)
with
3h¯2 d2 2 d ℓ(ℓ+1)
T (y)= − + − . (2.15)
ℓ 4m dy2 ydy y2 !
(+)
In Eq. (2.14), j (p y) is the spherical Bessel function and h (p y) is the spher-
ℓ 0 ℓ 0
ical Hankel function with the outgoing wave, where the outgoing (+) and the
incoming(−)sphericalHankel functionsaredefinedwiththesphericalNeumann
function n (p y) as
ℓ 0
(±)
h (x) =−n (x)±ij (x). (2.16)
ℓ ℓ ℓ
The function ω(e)(y), which plays a role of the source for the elastic compo-
nent in Eq. (2.13), is given by
∞
ω(e)(y) = x2dxφd(x)V(x)χ (x,y). (2.17)
α0
Z0
The explicit expression of the Green’s function Eq. (2.14) gives the asymp-
totic form of F(e)(y) as
F(e)(y) → h(+)(p y) T(e), (2.18)
y→∞ ℓ0 0
where T(e) is the elastic T-matrix amplitude defined by
4m ∞
T(e) = −p y2dyj (p y)ω(e)(y). (2.19)
0(cid:18)3h¯2(cid:19)Z0 ℓ0 0
The second term in the right hand side of Eq. (2.12) expresses three-body
breakupandclosed-channelcomponentsinthescattering.Inourformalism,these
components are treated by expandingthe Faddeev kernel with respect to a spec-
tator particle state of momentum p,
2
u (y;p) ≡ pj (py), (2.20)
ℓ ℓ
π
r
S.Ishikawa 5
which satisfies a complete relation
δ(y −y′) ∞
= dpu (y;p)u (y′;p). (2.21)
yy′ ℓ ℓ
Z0
(b,c)
The function φ (x,y) thereby is written as a Fourier-Bessel transformation:
α
∞
φ(b,c)(x,y) = dpu (y;p) η (x;p)−δ φd(x)C (p) . (2.22)
α ℓ α α,α0 α
Z0 h i
Here, η (x;p) is defined as
α
η (x;p) = hx|G |ωˆ i, (2.23)
α L α
where G is a two-body Green’s operator
L
1
G = (2.24)
L
E +iε−T (x)−V(x)
q L
with
¯h2 d2 2 d L(L+1)
T (x) = − + − . (2.25)
L m dx2 xdx x2 !
The energy of the two-body subsystem E is given by
q
3h¯2 ¯h2
E = E − p2 = q2, (2.26)
q
4m m
and the p-dependence of the functions arises through this relation.
ThebreakupcomponentstemsfromtheintegralofthefirstterminEq.(2.22)
for the range of 0 ≤ p ≤ p = 4mE/3h¯2. In this range the energies of both
c
the spectator particle and the twqo-body subsystem are positive or zero, and thus
the integral survives at infinite values of x and y, see Eq. (2.33) below and Refs.
[14, 15]. The rest of the integral of the first term in Eq. (2.22), i.e., p < p < ∞,
c
as well as the second term in Eq. (2.22) damp for large values of x and y because
the energy of the two-body subsystem is negative. In this sense, we call these
components closed.
The source term in Eq. (2.23), ωˆ (x;p), is written as
α
ωˆ (x;p) = V(x)χˆ (x;p), (2.27)
α α
∞
χˆ (x;p) = y2dyu (y;p)χ (x,y). (2.28)
α ℓ α
Z0
The second term of the right hand side in Eq. (2.22) appears as a counter
part of the subtraction and C (p) is defined as
α
1 ∞
C (p)= x2dxφd(x)ωˆ (x;p). (2.29)
α α
E −E
q d Z0
Theapparentsingularity in C (p) cancels that of the two-body Green’s operator
α
G , which will be numerically shown in the following section, and thus, we can
L
6 Faddeev-Kernelin Configuration Space
apply a standard quadrature to perform the p-integration in Eq. (2.22) as far as
the both terms are treated together.
In calculating η (x;p), we transform Eq. (2.23) to an ordinary differential
α
equation:
[E −T (x)−V(x)]η (x;p) = ωˆ (x;p) (2.30)
q L α α
with boundary conditions
(+)
h (qx) (0 ≤ p ≤ p )
η (x;p) ∝ L c . (2.31)
α x→∞( h(L+)(i|q|x) (pc < p < ∞)
A treatment of the two-body Green’s operator at three-body breakup re-
gion, 0 ≤ p ≤ p will be described in Appendix B:. We here only note that the
c
asymptotic form of η (x;p) is given by
α
−qm
η (x;p) → h(+)(qx) ¯h2 hψˆ (q)|ωˆ i, (2.32)
α x→∞ L 1(cid:16)−iK ((cid:17)q) L α
L
where ψˆ (x;q) is a two-body scattering solution with the standing wave bound-
L
ary condition and K (q) is a scattering K-matrix for the two-body scattering
L
(See Appendix B:).
(b,c)
The asymptotic form of φ (x,y) is evaluated by the saddle-point approxi-
α
mation [14, 15] as
φα(b,c)(x,y)x→∞,→x/y fixed−eπ4ii−L−ℓ 4K30 3/2 eRiK5/0R2 Bα(Θ), (2.33)
(cid:18) (cid:19)
where we introduce a hyper radius R and a hyper angle Θ as
4
R = x2+ y2, (2.34)
3
r
3
x =RcosΘ, y = RsinΘ, (2.35)
4
r
and K is given by
0
m
K = E. (2.36)
0 ¯h2
r
B (Θ) is the breakup amplitude defined as
α
1 m 1
B (Θ) = − hψˆ (q¯)|ωˆ i. (2.37)
α p¯¯h21−iK (q¯) L α
L
Here, the momenta q¯and p¯are given as
4
q¯= K cosΘ, p¯= K sinΘ. (2.38)
0 0
3
r
S.Ishikawa 7
3 Numerical Analyses and Results
3.1 Model source function
Inthis section, we presenta numerical example of thekernel operation described
in the preceding section with a model source function that is restricted to L =
ℓ = J = 0 state but carries a feature of the presence of three-body breakup
0
channel similarly as the one used in Ref. [12]:
eiK0R
χ (x,y) = (3.1)
α (R+R )5/2
0
with R = 5 fm.
0
We choose the 3S -component of the Malfliet-Tjon model as presented in
1
Ref. [1] for the potential V(x) and the incident nucleon energy of E = 14.1
Lab
MeV, which gives K =0.416 fm−1, p =0.480 fm−1, and p =0.550 fm−1. In nu-
0 c 0
merical calculations below, mesh points for x- and y-variables in described in
Appendix C: are used.
3.2 Elastic part
In Fig. 1, we plot the real part of the elastic source function ω(e)(y), Eq. (2.17),
calculated with the model function Eq. (3.1). As shown in this figure, ω(e)(y)
reveals a long-range behavior, which is given by
eipcy
ω(e)(y) ∝ , (3.2)
y→∞ y5/2
whose oscillation length 2π is about 13 fm.
pc
In calculating the elastic component F(e)(y), we treat this long-range behav-
ior by rewriting Eq. (2.13) as
F(e)(y)= −n (p y)T(y)+ij (p y)T(e)+j (p y) S(y)−Sˆ , (3.3)
ℓ0 0 ℓ0 0 ℓ0 0
(cid:16) (cid:17)
where we have defined T(y), S(y), and Sˆ by
4m y
T(y) = −p y′2dy′j (p y′)ω(e)(y′) → T(e)
0(cid:18)3h¯2(cid:19)Z0 ℓ0 0 y→∞
4m y
S(y) = −p y′2dy′n (p y′)ω(e)(y′) → Sˆ. (3.4)
0(cid:18)3h¯2(cid:19)Z0 ℓ0 0 y→∞
In numerical integration in Eq. (3.4), we need to be careful for oscillational
behaviors of the spherical Bessel and spherical Neumann functions as well as
ω(e)(y). This is done by spline interpolation technique used in Ref. [5] taking
into account of oscillational behavior of the integrand carefully.
For the use of Eq. (3.3), one needs converged values of T(y) and S(y) for
y → ∞. In Fig. 2, we plot the real part of T(y) for an example. As is expected
from the long-range behavior of ω(e)(y), the convergence of T(y) becomes very
8 Faddeev-Kernelin Configuration Space
4x10-4
2x10-4
]
)
y
(
e) 0
(w

[
e
R
-2x10-4
-4x10-4
0 50 100
y (fm)
Figure 1. Real part of theelastic source function ω(e)(y).
slow.However, fromthefunctionalformof Eq.(3.2),weexpectthatthefunction
T(y) behaves asymptotically as
ei(p0−pc)y
T(y) → t +t , (3.5)
y→∞ 0 1 y3/2
wheret and t are expansion coefficients and thecoefficient t is considered as a
0 1 0
convergedvalueofT(e).Thewavelengthevaluatedbythisequationis 2π = 90
p0−pc
fm, which is actually observed in Fig. 2. The fitting coefficients in Eq. (3.5) are
evaluated by a least square fit. To do this, the calculated values of T(y) in a
range of 80 ≤ y ≤ ymax (fm) are used. In Table 1, the dependence of the result
fit
on ymax is displayed. From the table, we set ymax = 1000 fm to get a converged
fit fit
result in five digits of accuracy, which is denoted by the dashed line in Fig. 2.
We remark that this result is obtained in spite of the fact the deviation of T(y)
from the converged value is still about 0.5 % at y = 1000 fm, which is not shown
in Fig. 2.
Here, we consider a range of the variables {x,y} to be used in calculations.
In the Faddeev equation, the function χ (x,y) is always accompanied by the
α
potential V(x), which means that we need to calculate this function within the
range of potential, x , for the x-variable. On the other hand, there is no restric-
R
tion for the y-variable. In actual, Table 1 demonstrates that we need to calculate
χ (x,y) for a large value of y, i.e., 1000 fm.
α
Suppose that we calculate χ (x,y) by Eqs. (2.8) and (2.10) with a func-
α
tion ξ (x,y), which is obtained in a preceding iteration step, for a range of
α
0 ≤ x≤ x ,0 ≤ y ≤ymax . The formulae, Eqs. (A.4) and (A.5), show that we
R fit
n o
S.Ishikawa 9
Deviation
-0.060
+8 %
+6 %
+4 %
] +2 %
)
y
T( -0.065
e[ 0 %
R
-2 %
-4 %
-6 %
-0.070
0 200 400 600
y (fm)
Figure 2. The real part of the function T(y).The obtained converged value is shown by the
dashed line. Dotted lines with the indices in the right hand side axis denote the deviation of
thereal part of T(y) from theconverged value.
Table 1. Thereal part
of the fitting coefficient
t .
0
ymax (fm) Re[t ]
fit 0
532 -6.5454
622 -6.5452
983 -6.5449
1073 -6.5448
1163 -6.5448
1193 -6.5448
need to prepare the function ξ (x,y) for a range of
α
1
0 ≤ x ≤ x +ymax
2 R fit
3 1
0 ≤ y ≤ x + ymax (3.6)
4 R 2 fit
to perform the exchange operation to obtain χ (x,y) for the above range. If we
α
set x =10 fm and ymax=1000 fm, this turns to be {0 ≤ x ≤ 1005 (fm),0 ≤ y ≤
R fit
507.5 (fm)}, which is rather huge.
Tofacilitate numericalcalculations, welimittherangeof calculating χ (x,y)
α
to {0 ≤ x ≤ x ,0 ≤ y ≤ y } by choosing the value of y adequately, and
R M M
10 Faddeev-Kernelin Configuration Space
approximate the value of χ (x,y) for y > y using a form of
α M
ei 43K0y a1(x) a2(x)
χ (x,y) = a (x)+ + , (3.7)
α 0≤x≤xR,y≥yM py2/5 (cid:18) 0 y y2 (cid:19)
where the coefficients a (x) are determined by a least square fit to χ (x,y) for
n α
y < y and for each value of x.
M
Withachoiceofx = 10fmandy = 80fm,bywhichtherangeforξ (x,y)
M M α
becomes {0 ≤ x ≤ 85 (fm),0 ≤ y ≤ 47.5 (fm)}, we obtain the equivalent results
for T(y) and its asymptotic value T(e) to the previously shown. This procedure
reduces the amount of calculations considerably without loss of accuracy, and
will be applied in the following analyses.
Together with the function S(y) and its asymptotic value Sˆ calculated sim-
ilarly, the elastic component F(e)(y) is constructed using Eq. (3.3), whose real
partis plotted in Fig. 3. Note that the effect of the slow convergence in T(y) and
S(y) functions appears as a small oscillation of the amplitude of F(e)(y) with
the wave length of 2π = 90 fm, which exists up to a large distance where the
p0−pc
convergences of T(y) and S(y) are achieved.

0.2
]
)
y
(
e) 0.0
((cid:1)

[
e
R
-0.2
0 100 200 300
y (fm)
Figure 3. The real part of the elastic function (e)(y).
F
3.3 Breakup and closed channel parts
First step in the calculation of the three-body breakup and closed channel con-
tributions is the Fourier-Bessel transformation of χ (x,y) with respect to the
α
coordinate y, Eq. (2.28). We again face the problem of slow convergence in the