aa r X i v : . [ nu c l - t h ] O c t Constrained correlated-Gaussians for hyperspherical calculations
Y. Suzuki
1, 2 and K. Varga Department of Physics, Niigata University, Niigata 950-2181, Japan RIKEN Nishina Center, Wako 351-0198, Japan Department of Physics and Astronomy, Vanderbilt University, Nashville, Tennessee 37235, USA (Dated: October 23, 2018)We formulate a hyperspherical approach within standard configuration interaction calculationsaiming at a description of large-scale dynamics of N -particle system. The channel wave functionand the adiabatic channel energy are determined by solving a hyperradius-constrained eigenvalueproblem of the adiabatic Hamiltonian. The needed matrix elements are analytically evaluatedusing correlated Gaussians with good orbital angular momentum and parity. The feasibility of theapproach is tested in three- α system. A spectrum of the adiabatic channel energies is determineddepending on the degree of localization of the basis functions. I. INTRODUCTION
The hyperspherical coordinate system is a natural ex-tension of the three-dimensional spherical polar coor-dinates to a set of N -particle coordinates. The hy-perspherical approach attempts to solve an N -particleSchr¨odinger equation by expressing the total wave func-tion as a product of the hyperradial and hyperangularparts and can be used to solve bound and continuumstate problems.The main advantage of the hyperspherical method isthat it provides a unified framework to describe quan-tum dynamics of complex reactions such as decay, fusionor fission. In other methods the choice of relevant co-ordinates is not trivial. For example, in nuclear fusioninitially the relative distance between the nuclei mightbe the most important coordinate, but later other coor-dinates will be more suitable and necessary. In the hyper-spherical approach the hyperradius captures all featuresof the complicated dynamical processes and describe dy-namical properties of the system emerging at differenthyperradial distances.Two realizations of the hyperspherical approach arewidely used (see Refs. [1–5] for reviews on the hyper-spherical approach and its applications). In one approachthe hyperangular part is expanded in terms of the hyper-spherical harmonics that are eigenfunctions of the hy-perangular kinetic-energy operator, and a coupled hy-perradial equation is solved by including the interactionof the particles. In another approach, often called theadiabatic hyperspherical approach, the adiabatic Hamil-tonian consisting of the hyperangular kinetic energy andthe interaction potential is diagonalized first to obtainthe adiabatic channel energies and channel wave func-tions. The total wave function is then expanded in termsof the basis set of the adiabatic channel wave functions.The adiabatic channel energies give hints on how the sys-tem responds as a function of hyperradial distances.The advantage of the first approach is that the hy-perspherical harmonics are known, but two difficultiesmay hinder the application. One problem is that con-vergence of the hyperspherical harmonics expansion is slow even for short-range potentials [6], and it becomesprohibitively slow when a long-range potential, like theCoulomb coupling potential, acts at large hyperradialdistances. This slowness is related to the fact that thehyperangular kinetic energy and the interaction poten-tial do not commute [7, 8]. The slow convergence causeshuge discrepancies, e.g. in the triple- α reaction rate atlow temperatures [9–12]. Another problem is that solv-ing the coupled differential equation in the hyperradialcoordinate may become hard when a number of avoidedcrossings occur in the potential energy curves.Although the basic idea of the hyperpsherical methodis not limited to three-body systems, its extension tomore-particle system is impeded by the lack of appropri-ate basis functions that can be flexibly used in the hyper-spherical calculation. References [13–19] discuss recentdevelopments in going beyond the three-body problems.In both realization of the hyperspherical approach, onecalculates the matrix element of an operator O , h Ψ ′ |O| Ψ i ρ = R , (1)where h . . . i ρ = R indicates that the matrix element is tobe evaluated by integrating in all the coordinates butthe hyperradius ρ , which is fixed to R . The integral oftype (1) is hard to evaluate in general because specify-ing the hyperangle coordinates for the N -particle systemis considerably involved and integrating in those coordi-nates requires many-dimensional integrations. Althoughsome progress has recently been made with correlatedGaussian (CG) basis functions [20–23], the total orbitalangular momentum is limited to L = 0 and 1.In this paper we will examine the possibility of usingthe CG as suitable basis functions in hyperspherical cal-culations. The CG proposed many years ago [24, 25] isextended to describe motion with non-zero total orbitalangular momentum, especially with the help of the globalvector representation [26–28]. Together with the stochas-tic variational method [26, 27, 29] to select efficiently theparameters of the CG, many problems have accuratelybeen solved with the CG. See, e.g. Refs. [30–32] for somerecent applications of CG.We attempt to formulate the hyperspherical approachin standard configuration interaction calculations. Fol-lowing the spirit of the second realization of the hyper-spherical approach, we set up a number of basis func-tions that are expected to be important at ρ ≈ R ,calculate the matrix elements of the adiabatic Hamil-tonian using the full coordinate integration instead ofEq. (1), and determine both the channel wave functionand the adiabatic channel energy at ρ ≈ R by solving ahyperradius-constrained eigenvalue equation. We showthat this scheme can be achieved using the CG basisfunctions. The emphasis of this paper is not on solv-ing a specific problem with the hyperspherical approachbut on carefully examining its feasibility and discussingproblems that may occur.We present our formulation in Sec. II, and show inSec. III how to evaluate the needed matrix elements. InSec. IV we test our method in three- α system that is thesimplest possible system but contains all the complexitiesnevertheless. Section V is a summary and discussions. II. SCHR ¨ODINGER EQUATION INHYPERSPHERICAL APPROACHA. Hyperspherical coordinates
Let r i ( i = 1 , . . . , N ) denote the position coordinate ofthe i th particle. The mass m of all particles is assumedto be the same, although the case of unequal mass can betreated by defining mass-scaled coordinates. We define aset of relative coordinates, x i ( i = 1 , . . . , N − x i = r ii + 1 (cid:16) r i +1 − i i X j =1 r j (cid:17) . (2)The set x i together with the center of mass (c.m.) coor-dinate, x N ≡ R cm = P Ni =1 r i /N , defines a transforma-tion matrix U from the single-particle coordinates to therelative and c.m. coordinates: x i = N X j =1 U ij r j ( i = 1 , . . . , N ) . (3)Conversely, r i is expressed as r i = P Nj =1 U − ij x j .The square of the hyperradius ρ is defined by ρ = N X i =1 ( r i − R cm ) = 1 N N X j>i =1 ( r i − r j ) , (4)which is equal to ρ = N − X i =1 x i . (5)Let Ω denote a set of the hyperangle coordinates con-structed from dimensionless coordinates, ξ i = x i /ρ ( i = 1 , . . . , N − P N − i =1 ξ i = 1.The volume element for integration excluding R cm is d x ≡ d x . . . d x N − = ρ d − dρd Ω , (6)where d = 3( N −
1) (7)is the degree of freedom excluding the c.m. motion.Since ρ /N is the mean-square-radius operator, ρ mea-sures the global size of the system. Or ρ is a kind of col-lective coordinate responding to a large-scale change ofthe system [33]. Suppose that the system develops into f subsystems or clusters, each of which consists of N i particles ( P fi =1 N i = N ). ρ is rewritten as ρ = f X i =1 ρ i + ρ (8)with ρ i = N i X j =1 ( r ( i − j − R i ) , ρ = f X i =1 N i ( R i − R cm ) , (9)where ( i −
1) = P i − k =1 N k with (0) = 0, and R i is thec.m. coordinate of the i th cluster. ρ i is the hyperradiusof the i th cluster, and ρ rel stands for the hyperradiusthat measures the spatial extension of the relative mo-tion of the clusters. In such phenomena that include theformation of f subsystems, the contribution of P fi =1 ρ i to ρ remains finite, whereas ρ takes an increasinglylarge value as ρ increases. Moreover, since ρ is invariantwith respect to the number of clusters f , the hyperspher-ical coordinates have the unique advantage that they cantreat any decomposition of the system in a unified way. B. Equation of motion in hyperspherical approach
The Hamiltonian H of the system consists of the ki-netic energy T and the interaction potential V : H = T + V. (10)With the c.m. kinetic energy being subtracted, T reads T = − ~ m N − X i =1 ∂ ∂ x i , (11)and it is separated into hyperradial ( T ρ ) and hyperangu-lar ( T Ω ) parts, T = T ρ + T Ω : T ρ = − ~ m (cid:16) ∂ ∂ρ + d − ρ ∂∂ρ (cid:17) , T Ω = ~ Λ mρ , (12)where Λ is the squared grand angular momentum oper-ator that can in principle be expressed in terms of thehyperangle coordinates and their derivatives. The adia-batic Hamiltonian H Ω = T Ω + V (13)is a kind of a generalized potential. As usual, V andconsequently H Ω is assumed to contain no derivative op-erator with respect to ρ .Let the total wave function Ψ JMπ of the system belabeled by the total angular momentum J , its z com-ponent M , the parity π . The Schr¨odinger equation, H Ψ JMπ = E Jπ Ψ JMπ , reads as( T ρ + H Ω )Ψ JMπ = E Jπ Ψ JMπ . (14)The channel wave function Φ JMπν ( ρ, Ω) and the adiabaticchannel energy (or adiabatic potential) U Jπν ( ρ ) are de-fined by solving the eigenvalue problem of H Ω , H Ω Φ JMπν ( ρ, Ω) = U Jπν ( ρ )Φ JMπν ( ρ, Ω) . (15)In order to exploit the fact that ρ can be treated asa c -number in Eq. (15), we need the matrix element h Φ JMπν | H Ω | Φ JMπν ′ i ρ = R . Since its evaluation is, however,hard as already mentioned, we take a different route.We set up a number of independent basis functionsΦ JMπl ( l = 1 , . . . , M ) that satisfy h Φ JMπl | ρ | Φ JMπl i = R , (16)and assume that the ν th ‘channel wave function’ Φ JMπRν at h ρ i = R is given as a combination of the basis functionsΦ JMπRν = M X l =1 c JπRν,l Φ JMπl . (17)The coefficients c JπRν,l are determined by solving the con-strained eigenvalue problem h Φ JMπl | H Ω − U JπRν | Φ JMπRν i = 0 ( l = 1 , . . . , M ) , (18)subject to h Φ JMπRν | ρ | Φ JMπRν i = R . (19)Both Φ JMπl and Φ
JMπRν are normalized, and the matrixelements in Eqs. (16), (18), and (19) are evaluated byintegrating in all the coordinates. Appendix A shows howto determine the adiabatic channel energies U JπRν at R andmutually orthogonal channel wave functions Φ JMπRν .We calculate the channel wave functions at a numberof mesh points R i and assume the total wave function tobe approximated by their combinationsΨ JMπ = X iν χ Jπiν Φ JMπR i ν . (20)Equation (14) reduces to the following equation for χ Jπiν : X jν ′ h Φ JMπR i ν | T ρ + H Ω − E Jπ | Φ JMπR j ν ′ i χ Jπjν ′ = 0for all i and ν. (21) The condition (16) is necessary to look for a suitablebasis set at R because each piece of H Ω shows different ρ -dependence and hence such set may change dependingon ρ . Short-ranged interactions in V become importantin the region of small ρ , while long-ranged interactionslike the Coulomb potential contribute at large ρ as well. T Ω is also long-ranged. Moreover, since the Coulombpotential and T Ω do not commute each other, one has totake account of both terms simultaneously [7].As a measure of the localization of a wave function Φ,we introduce the standard deviation σ of ρ : σ = h Φ | ( ρ − h Φ | ρ | Φ i ) | Φ ih Φ | ρ | Φ i = h Φ | ρ | Φ ih Φ | ρ | Φ i − . (22)Φ JMπRν is obtained as a combination of Φ
JMπl s. Therefore,even though Φ
JMπl s are all set to have σ ’s within a certainrange, it may happen that the σ value of Φ JMπRν is farbeyond its range. To obtain U JπRν around ρ ≈ R , it isuseful to check the σ value of Φ JMπRν . We will discuss thisproblem later.Once Φ
JMπR i ν s are determined at R i ’s, Eq. (21) can besolved in a standard linear algebra. Note that the ma-trix element of T ρ is already available at the stage ofsolving the eigenvalue problem of H Ω . This is in sharpcontrast to the standard hyperspherical method where nohyperradial function is employed and thus one has to usenumerical differentiations with respect to ρ or e.g. slowvariable discretization method [34, 35].In what follows we omit the superscripts JM π . III. CORRELATED GAUSSIAN ASHYPERSPHERICAL BASIS FUNCTIONA. Correlated Gaussian and its generating function
We adopt the CG as the basis function. We use matrixnotations to make equations compact. For example, x denotes a column vector of dimension ( N −
1) whose i thelement is x i . A tilde symbol e indicates a transposeof a column vector or a matrix, e.g. e x is the row vectorand ρ may be written as e xx , where the scalar productof 3-dimensional vectors is implicitly understood: e xx = P N − i =1 x i · x i = P N − i =1 x i .The CG with the total orbital angular momentum L and its z component M reads f uAKLM ( x ) = N uAKL | e u x | K + L Y LM ( ce u x ) e − e x A x , (23)where a column vector u = ( u i ) of dimension ( N −
1) anda symmetric, positive-definite ( N − × ( N −
1) matrix A = ( A ij ) are both (variational) parameters to character-ize the CG. Both A and u are assumed to be real in thispaper. The exponential part, e − e x A x , is invariant underthe coordinate rotation, whereas the spherical harmonics Y LM describes the rotational motion through the globalvector, e u x = P N − i =1 u i x i [26–28, 36, 37]. ce u x stands forthe polar and azimuthal angles of e u x . N uAKL is the nor-malization constant determined from h f uAKLM | f uAKLM i = 1. K is a non-negative integer parameter related to thelocalization in ρ motion of the CG [38]. It should benoted that the CG has simple hyperadial dependence f uAKLM ( x ) ∼ ρ κ e − ρ e ξ A ξ ( κ = 2 K + L ) , (24)where ξ = ( ξ i ) is a column vector of dimension ( N − T ρ .Let us introduce the generating function for the CG, g ( s , A, x ) = e − e x A x + e sx , (25)where s = ( s i ) is a column vector of dimension ( N − s i . With a choice of s i = αu i e , where α is an auxiliary real parameter and e is a three-dimensional unit vector ( e = e · e = 1), theCG is generated as follows[26, 27]: f uAKLM ( x ) = N uAKL B KL Z d e Y LM (ˆ e ) × (cid:16) d K + L dα K + L g ( αu e , A, x ) (cid:17) α =0 (26)with B KL = 4 π (2 K + L )!2 K K ! (2 K + 2 L + 1)!! . (27)Here ( ) α =0 indicates that α is set to zero after thedifferentiation. B. Basic matrix elements
The CG matrix elements for various operators areavailable in the literature [26–28, 38]. We recapitulatesthe basic procedure to derive them with emphasis onthe relationship to the Gauss hypergeometric function(GHF) [39, 40], which has hitherto never been recognized.Applying Eq. (26) leads to the CG matrix element: h f u ′ A ′ K ′ LM | ˆ O| f uAKLM i = N u ′ A ′ K ′ L B K ′ L N uAKL B KL Z d e ′ Z d e Y ∗ LM ( ˆ e ′ ) Y LM (ˆ e ) × (cid:16) d K ′ + L +2 K + L dα ′ K ′ + L dα K + L Z d x e − e x B x + e vx O ( x ) (cid:17) α =0 α ′ =0 . (28)Here O ( x ) is determined by acting ˆ O on g ( s ; A, x ) or f uAKLM ( x ). The matrix B and the vector v are defined by B = 12 ( A + A ′ ) , v = s + s ′ , (29)where s = αu e and s ′ = α ′ u ′ e ′ . For a class of operators, the integral in Eq. (28) overthe whole region of x takes the form Z d x e − e x B x + e vx O ( x ) = P O (cid:16) π N − det B (cid:17) e e v B − v . (30)Appendix B lists some examples of O ( x ) and P O . In allthose cases, P O consists of terms with the form T kk ′ l ( u ′ A ′ , uA ) α k α ′ k ′ ( αα ′ e · e ′ ) l , (31)each of which is characterized by non-negative integers, k , k ′ , l , and the coefficient T kk ′ l ( u ′ A ′ , uA ). The exponentin Eq. (30) is14 e v B − v = pα + p ′ α ′ + qαα ′ e · e ′ , (32)where p = 14 e uB − u, p ′ = 14 e u ′ B − u ′ , q = 12 e uB − u ′ . (33) e uAu ′ or ( e uAu ′ ) stands for the inner product, P N − i,j =1 u i A ij u ′ j . Expanding e e v B − v in a power seriesof α and α ′ , and combining it with the term of Eq. (31),we perform the operation in Eq. (28), obtaining the con-tribution of term (31) to the matrix element as follows: h f u ′ A ′ K ′ LM |O| f uAKLM i ∼ N u ′ A ′ K ′ L B K ′ L N uAKL B KL (cid:16) π N − det B (cid:17) × T kk ′ l ( u ′ A ′ , uA )(2 K + L )!(2 K ′ + L )! × n X n = n p K − k − n p ′ K ′ − k ′ − n q n + L − l B nL ( K − k − n )!( K ′ − k ′ − n )!(2 n + L − l )! , (34)where n and n are given by n = min( K − k, K ′ − k ′ ) ,n = n L ≧ l (cid:2) l − L +12 (cid:3) for l > L. (35)Here Gauss’s symbol [ x ] stands for the greatest integerthat is less than or equal to x .The sum in Eq. (34) can be expressed with the GHFas follows. By using B nL (27), the sum reduces to n X n = n p K − k − n p ′ K ′ − k ′ − n q n + L − l B nL ( K − k − n )!( K ′ − k ′ − n )!(2 n + L − l )!= 4 π p K − k p ′ K ′ − k ′ q L − l (2 z ) n ( K − k − n )!( K ′ − k ′ − n )!(2 L + 2 n + 1)!! × n − n X m =0 ( − K + k + n ) m ( − K ′ + k ′ + n ) m m !( L + n + ) m P L,n l ( m ) z m , (36)where ( a ) m is Pochhammer’s symbol( a ) m = Γ( a + m )Γ( a ) (37)expressed with the Gamma function Γ. If a is negative,( a ) m = ( − m Γ( − a +1) / Γ( − a − m +1). If a is a negativeinteger, a = − k , ( − k ) m = 0 for m > k . z in Eq. (36) isdefined by z = q pp ′ = ( e uB − u ′ ) ( e uB − u )( e u ′ B − u ′ ) , (38)and takes a value in the interval [0 , P L,n l ( m ) inEq. (36) is a polynomial of m with the order l − n , P L,n l ( m ) = m !(2 m + L + 2 n )!( m + n )!(2 m + L − l + 2 n )! . (39)Because of m i z m = ( z ddz ) i z m for any non-negative inte-ger i , P L,n l ( m ) z m may be replaced by P L,n l (cid:0) z ddz (cid:1) z m ,which makes Eq. (36) further compact: n X n = n p K − k − n p ′ K ′ − k ′ − n q n + L − l B nL ( K − k − n )!( K ′ − k ′ − n )!(2 n + L − l )!= 4 π p K − k p ′ K ′ − k ′ q L − l (2 z ) n ( K − k − n )!( K ′ − k ′ − n )!(2 L + 2 n + 1)!! × P L,n l (cid:0) z ddz (cid:1) γ K − k − n ,K ′ − k ′ − n ,L + n ( z ) . (40)Here γ K,K ′ ,L ( z ), introduced in Ref. [38], is nothing butthe GHF γ K,K ′ ,L ( z ) = F ( − K, − K ′ ; L + ; z ) , (41)which is actually a polynomial of z with the ordermin( K, K ′ ) because K and K ′ are both non-negative in-tegers in the present case.Equations (34) and (40) constitute a basic formula tocalculate the matrix element. Let us consider the overlapmatrix element, for which O ( x ) = 1, P O = 1, k = k ′ = l = 0, T ( u ′ A ′ , uA ) = 1, leading to h f u ′ A ′ K ′ LM | f uAKLM i = N u ′ A ′ K ′ L B K ′ L N uAKL B KL (cid:16) π N − det B (cid:17) (2 K + L )!(2 K ′ + L )! × π p K p ′ K ′ q L K ! K ′ !(2 L + 1)!! γ K,K ′ ,L ( z ) . (42)In the diagonal case of u ′ = u, A ′ = A, K ′ = K , z isunity and γ K,K,L (1) is easily obtained by using F ( a, b ; c ; 1) = Γ( c )Γ( c − a − b )Γ( c − a )Γ( c − b ) , (43)which is valid provided that Re ( c − a − b ) >
0. Thenormalization constant is then given by N uAKL = vuut A ) √ π N − Γ(2 K + L + )( e uA − u ) K + L . (44) TABLE I: F KKLkk ′ l (1) for some sets of ( k, k ′ , l ). SeeEq. (47). Note that F KKLk ′ kl (1) = F KKLkk ′ l (1). M and M stand for M = K + L + and M = 2 K + L + ,respectively. k k ′ l F KKLkk ′ l (1)0 0 0 11 0 0 K M M M M − L −
12 0 0 K ( K − M ( M − M ( M − K M M ( M − K M ( M − M ( M − − K ( L + 1) M M M ( M − M ( M − − L + 1) M M + ( L + 1)( L + 2) Substitution of Eqs. (27), (33), and (44) into Eq. (42) andthe use of Eq. (43) completes the overlap matrix element: h f u ′ A ′ K ′ LM | f uAKLM i = (cid:16) det AA ′ (det B ) (cid:17) (cid:16) e uB − u e uA − u (cid:17) (2 K + L ) (cid:16) e u ′ B − u ′ e u ′ A ′− u ′ (cid:17) (2 K ′ + L ) × (cid:16) e uB − u ′ | e uB − u ′ | √ z (cid:17) L γ K,K ′ ,L ( z ) p γ K,K,L (1) γ K ′ ,K ′ ,L (1) . (45)Combining Eqs. (34), (40), and (42) enables us to ex-press the contribution of term (31) in relation to the over-lap matrix element: h f u ′ A ′ K ′ LM |O| f uAKLM i ∼ h f u ′ A ′ K ′ LM | f uAKLM i× T kk ′ l ( u ′ A ′ , uA ) p − k p ′− k ′ q − l F KK ′ Lkk ′ l ( z ) , (46)where F KK ′ Lkk ′ l ( z ) = K ! K ′ !( K − k − n )! ( K ′ − k ′ − n )! ( L + ) n × z n γ K,K ′ ,L ( z ) P L,n l (cid:0) z ddz (cid:1) γ K − k − n ,K ′ − k ′ − n ,L + n ( z ) . (47)Equations (45), (46), and (47) give a powerful formulafor the matrix element. We only need to determine T kk ′ l ( u ′ A ′ , uA ), which contributes to the matrix elementprovided that both K − k − n and K ′ − k ′ − n arenon-negative.Small values of k , k ′ , and l are usually needed. Forexample, in all the classes of Eq. (B2), possible sets of( k, k ′ , l ) are (0 , , , (1 , , , (0 , , , (0 , , , (2 , , , , , (1 , , , (0 , , , (0 , , , , P L,n l ( m ) turns out to be simple. For l = 0, P L, ( m ) = 1. For l = 1, P L, ( m ) = 2 m + L ( L ≧
1) and P L, ( m ) = 2 ( L = 0). For l = 2, P L, ( m ) = (2 m + L − m + L ) ( L ≧
2) and P L, ( m ) =2(2 m + 2 L + 1) ( L = 0 , P L,n l ( z ddz ) is given, itsaction on γ K,K ′ ,L ( z ) is performed by using z ddz γ K,K ′ ,L ( z ) = (cid:0) L + 12 (cid:1)(cid:2) γ K,K ′ ,L − ( z ) − γ K,K ′ ,L ( z ) (cid:3) , (48)which is derived from the well-known formulas involvingthe GHF. Table I tabulates F KKLkk ′ l (1) for the above cases.With C set to the unit matrix in Eqs. (B5) and (B6)and using Table I, the expectation values of ρ and ( ρ −h ρ i ) are given by ( κ = 2 K + L ) h f uAKLM | ρ | f uAKLM i ≡ h ρ i = 32 Tr A − + κ e uA − u e uA − u , (49) h f uAKLM | ( ρ − h ρ i ) | f uAKLM i = 32 Tr A − + 2 κ e uA − u e uA − u − κ (cid:16) e uA − u e uA − u (cid:17) . (50)The σ value of Eq. (22) is readily obtained for f uAKLM . C. Hamiltonian matrix element
We show how to calculate the matrix element of H Ω .First we note that the relative distance vector, r i − r j , isexpressed as a combination of x k , r i − r j = N − X k =1 ( U − ik − U − jk ) x k ≡ g ω ( ij ) x , (51)where ω ( ij ) is a column vector of dimension ( N − r i − r j ) = e x T ( ij ) x , where T ( ij ) = ω ( ij ) g ω ( ij ) is a symmetric ( N − × ( N −
1) matrix. A Gaussianpotential e − a ( r i − r j ) is expressed as e − a e x T ( ij ) x , and itsmatrix element reduces to the overlap (45): h f u ′ A ′ K ′ LM | e − a e x T ( ij ) x | f uAKLM i = G u ′ A ′ : uAK ′ L : KL ( aT ( ij ) ) , (52)where G u ′ A ′ : uAK ′ L : KL ( T )= R u ′ A ′ : A ′ + TK ′ L R uA : A + TKL h f u ′ A ′ + TK ′ LM | f u A + TKLM i (53)with R uA : A ′ KL = N uAKL N u A ′ KL = (cid:16) det A det A ′ (cid:17) (cid:16) e uA ′− u e uA − u (cid:17) K + L . (54)The matrix element of three-body force of Gaussian formfactor can be obtained in a similar way.The matrix elements of Coulomb and Yukawa poten-tials are obtained by applying the above result [36]. Forexample, by expressing the Yukawa potential as1 r e − µr = 2 √ π Z ∞ dt exp (cid:0) − t r − µ t (cid:1) , (55) its matrix element is obtained by a numerical integrationof Eq. (52) with an appropriate change of the range pa-rameter a . Equation (52) is valid for not only T ( ij ) butany positive-definite symmetric matrix. For example, us-ing the unit matrix I we obtain h f u ′ A ′ K ′ LM | ρ | f uAKLM i = Z ∞ dt G u ′ A ′ : uAK ′ L : KL ( tI ) , (56)which is computed with e.g. Gauss-Laguerre quadrature.We turn to the hyperangular kinetic energy T Ω . Weobtain its matrix element without expressing Λ in termsof Ω, but in the indirect way [20] that utilizes the identity, T Ω = T − T ρ . The matrix elements of T and T ρ arerespectively obtained as follows. As for T , we start from T g ( s , A, x )= ~ m (3Tr A − e ss + 2 e s A x − e x A x ) g ( s , A, x ) . (57) T kk ′ l ( u ′ A ′ , uA ) contributed by each term of Eq. (57) isread from Appendix B. The use of Eq. (46) leads to h f u ′ A ′ K ′ LM | T | f uAKLM i = ~ m h f u ′ A ′ K ′ LM | f uAKLM i h
32 Tr B −
32 Tr C C − ( e uu − e uC u + e uC u ) 1 e uB − u F KK ′ L ( z ) − ( e u ′ u ′ + 2 e u ′ C u ′ + e u ′ C u ′ ) 1 e u ′ B − u ′ F KK ′ L ( z )+ ( e uu ′ + e uC u ′ − e u ′ C u − e uC u ′ ) 1 e uB − u ′ F KK ′ L ( z ) i , (58)where C, C , and C are the matrices defined by C = 12 ( A − A ′ ) , C = CB − , C = B − C B − . (59)As for T ρ , we use Eq. (24) to obtain the relation T ρ f uAKLM = − ~ mρ h κ + ( d − κ − (2 κ + d ) e x A x + ( e x A x ) i f uAKLM . (60)As in Eq. (56), the matrix element of T ρ is obtained byperforming the following integration: h f u ′ A ′ K ′ LM | T ρ | f uAKLM i = − ~ m Z ∞ dt R u ′ A ′ : A ′ + tIK ′ L R uA : A + tIKL × h f u ′ A ′ + tIK ′ LM | h κ + ( d − κ − (2 κ + d ) e x A x + ( e x A x ) i | f u A + tIKLM i , (61)where the matrix elements of e x A x and ( e x A x ) arereadily available from Eqs. (B5) and (B6). FollowingRefs. [20–22], the matrix element of T Ω is given as h f u ′ A ′ K ′ LM | T Ω | f uAKLM i = 12 (cid:16) h f u ′ A ′ K ′ LM | T − T ρ | f uAKLM i + h f uAKLM | T − T ρ | f u ′ A ′ K ′ LM i (cid:17) . (62)All the matrix elements needed to solve Eq. (21) arethus available in the CG basis functions. The presentapproach thus reduces the whole task to a standard linearalgebra of matrices in place of the coupled differentialequation commonly used in the hyperspherical approach. D. Permutation symmetry
The permutation symmetry for identical particles hasto be imposed on the wave function. Its incorporation inthe CG is very easy [27, 38, 41].The permutation P induces the coordinate transfor-mation: x → T P x , where the ( N − × ( N −
1) matrix T P is easily determined. Since P just rearranges the la-bels of r i , ρ remains unchanged (see Eqs. (4) and (5)): ρ = e xx → g T P x T P x = e x f T P T P x = e xx , concluding f T P T P = I. (63)The CG acted by P transforms to P f uAKLM ( x ) = N uAKL | f u P x | K + L Y LM ( df u P x ) e − e x A P x = N uAKL N u P A P KL f u P A P KLM ( x ) , (64)where u P = f T P u, A P = f T P AT P . (65)Since det A P = det A and f u P A − P u P = e uA − u , Eq. (44)confirms N u P A P KL = N uAKL , establishing P f uAKLM ( x ) = f u P A P KLM ( x ) . (66)The CG keeps its functional form under the permu-tation, and its effect results in simply changing the CGparameters, u and A , as in Eq. (65). The basis functionΦ uAKLM that is constructed from f uAKLM ( x ) and satisfiesthe symmetry requirement is given byΦ uAKLM = X P ǫ P f u P A P KLM ( x ) , (67)where ǫ P is the phase of P . IV. TEST OF THREE- α SYSTEM
In order to learn similarity to and dissimilarity fromthe usual adiabatic channel energy, we use the same Hamiltonian as that of Ref. [11]. The mass of the α par-ticle is ~ /m =10.5254408 MeV fm , and the charge con-stant is e = 1.4399644 MeV fm. The two-body potential V αα ( r ) consists of a modified Ali-Bodmer potential [42]and the Coulomb potential: V αα ( r ) = 125 e − r / . − . e − r / . + 4 e r erf (0 . r ) , (68)where the length and energy are given in units of fm andMeV, respectively. The three-body potential is chosen tobe a hyperscalar potential, V ααα = v e − a e xx , (69)where the range parameter a is √ /R with R =2.58 fm,and the potential strength v is L -dependent: It is − .
737 MeV for L = 0 to reproduce the Hoyle reso-nance energy, and − .
463 MeV for L = 2 to fit thelowest 2 + state energy of C. A. Specification of correlated-Gaussian parameters
We use Φ uAKLM , Eq. (67), as the basis functions Φ R i ,l .The label l stands for K, u , and A . u contains just 1 pa-rameter, assuming that u is normalized: u = sin ζ, u =cos ζ (0 ≦ ζ < π ). ζ is discretized by M ζ meshes. Thematrix A for three-body system contains 3 parameters, A , A (= A ) , A . It may be prescribed with threeparameters ( d , d , d ) as e x A x = X j>i =1 d ij ( r i − r j ) . (70)Roughly speaking, d ij controls the distance between par-ticles i and j . In analogy to the prescription used inRefs. [11, 12], we specify d ij by two angles θ (0 ≦ θ < π/ φ (0 ≦ φ ≦ π ) that define the ‘shape’ of three parti-cles: d = ¯ d h θ cos (cid:0) φ + 23 π (cid:1)i ≡ ¯ d λ + ,d = ¯ d [1 + sin θ cos φ ] ≡ ¯ d λ ,d = ¯ d h θ cos (cid:0) φ − π (cid:1)i ≡ ¯ d λ − . (71) θ = 0 and θ = π/ φ can be restricted to[0 , π/ θ and φ by M θ and M φ meshes.The matrix A reads as A = A / ¯ d , where A = (cid:16) λ + + ( λ + λ − ) − √ ( λ − λ − ) − √ ( λ − λ − ) ( λ + λ − ) (cid:17) , (72)and ¯ d is determined from the constraint (16). B. Results
As defined in Eq. (16), R is a c -number representing p h ρ i . In Refs. [11, 12], R stands for both the hyper-radius operator and its value, although the hyperradiusthere corresponds to 3 / ρ of the present paper. To avoidconfusion, we employ a point- α root-mean-square (rms)radius R rms as a length scale, R rms = r h ρ i , (73)which is computed as 3 − R from R in [11, 12].Figure 1 displays K -dependence of the minimum ex-pectation value of H Ω calculated by a single configura-tion with L π = 0 + . For each K , u and A are varied onthe meshes discretized with M ζ , M θ , and M φ , subject to R rms = 1 .
54 fm. The minimum of the adiabatic channelenergy occurs around that rms value [11]. The mini-mum of the curve is − .
62 MeV at K = 4 and graduallyincreases with K . The contribution of T Ω to the min-imum expectation value increases as K increases, whilethe sum of the potentials, V + V + V C , shows a moder-ate change for K ≧
2, probably because it is determinedmainly by the global size of the system. The curve la-beled H is the sum of the minimum energy of H Ω andthe expectation value of T ρ , that is, the total energy, cal-culated by the optimal configuration. The expectationvalue of T ρ increases from 7 to 27 MeV as K increasesfrom 0 to 20. Figure 2 is the same as Fig. 1 but for L π = 2 + . The configuration is again constrained to sat-isfy R rms = 1 .
54 fm. The minimum energy of H Ω is 6.79MeV at K = 2.Table II lists some properties of the single configura-tion used in Figs. 1 and 2. It is noted that the standarddeviation σ decreases as K increases. The configurationwith K = 4 giving the energy minimum for L = 0 has σ = 0 . σ value corresponds tothe degree of localization, ( √ σ − R rms ≈ .
064 fmaround R rms . If we want to use more localized configu-rations, we have to increase K . Since the overlap withthe K = 4 configuration decreases very slowly as listed inOverlap column, the energy loss may not be very large.In L = 2 case, the K dependence of σ and the overlapintegral appears to decrease faster than the L = 0 case.Figure 3(a) plots the minimum expectation value of H Ω for L π = 0 + as a function of R rms . The minimum energyis obtained by a single configuration determined similarlyto the case of Fig. 1, but with slightly finer meshes. Welearn how each term of H Ω responds to the expansionof the system as R rms increases. The H Ω curve showsa minimum around R rms =1.6 fm, and reaches a broadtiny peak at 14.6-14.7 fm, where the contribution of eachpiece of H Ω displays a sudden change as magnified inFig. 3(b). Before the peak, V C , T Ω , and V are maincontributors to the H Ω curve, whereas after the peakboth contributions of T Ω and V get small and V C playsa dominant role. Note, however, that the contribution of E ( M e V ) H Ω T Ω V V V C H FIG. 1: (Color online) Minimum expectation value of H Ω calculated by a single configuration, Φ uAKL =0 M =0 ,Eq. (67), as a function of K . The minimum is searchedfor by varying θ, φ, ζ on the meshes discretized with M θ = 30, M φ = 20, and M ζ = 30 under the constraintthat R rms is kept to 1.54 fm. The contributions of T Ω ,the nuclear potentials (two-body V , and three-body V ) as well as the Coulomb potential ( V C ) to theminimum energy are also drawn. The curve denoted H is the variation of the total energy. E ( M e V ) H Ω T Ω V V V C H FIG. 2: (Color online) The same as Fig. 1 but for L π = 2 + . V persists up to large distances beyond 14 fm, despitethe fact that the range of V αα is much shorter than thatvalue. This long-range effect is due to the αα resonance.Although the minimum expectation value of H Ω changes smoothly with R rms as seen in Fig. 3(a), thecontributions of T Ω and V show some kinks, especiallywhen R rms changes from 1.6 to 1.7, 2.4 to 2.5, and 4.8 to4.9 fm. At these points the optimal K value also changesas follows: 4 →
3, 3 →
2, and 2 →
1, respectively. However,the minimum expectation value of H Ω is often not verysensitive to the change of K but several K configurationsgive almost equal results, whereas the contribution of T Ω seems to be more sensitive to K . This is understoodfrom the degree of localization of the CG. In fact, the σ TABLE II: The properties of the single configurationused in Figs. 1 and 2. σ is the standard deviation andOverlap is the overlap integral with K = 4 ( L = 0) or K = 2 ( L = 2) basis function. L = 0 L = 2 K σ
Overlap σ Overlap0 0.577 0.923 0.477 0.9771 0.480 0.973 0.431 0.9952 0.450 0.991 0.400 1.0003 0.435 0.998 0.374 0.9964 0.419 1.000 0.378 0.9855 0.413 0.999 0.381 0.9786 0.403 0.996 0.373 0.9527 0.399 0.992 0.365 0.9418 0.394 0.987 0.338 0.9299 0.389 0.981 0.338 0.91910 0.382 0.975 0.332 0.90911 0.370 0.96712 0.370 0.96213 0.358 0.95214 0.364 0.95015 0.354 0.94016 0.356 0.93617 0.346 0.92718 0.342 0.91919 0.334 0.90920 0.335 0.905 value of Φ uAKLM decreases with increasing K , and hencethe contribution of T Ω tends to increase.Now we mix various configurations to solve Eq. (18).The constrained equation is solved at the following fourpoints: R rms =1.6, 2.5, 5.0, 14.5 fm. The lowest adiabaticchannel energy of Ref. [11] exhibits different character atthese points, a steep slope close to the minimum, and abroad plateau close to the three- α threshold. The CG ba-sis functions are generated by including different K, θ, φ ,and ζ parameters. K is tested up to 20. The mesh pointsare discretized with M θ = 30, M φ = 21, and M ζ = 45.To avoid possible linear-dependence of the generated ba-sis functions, we exclude any basis function that has over-lap of more than 0.95 with other basis functions. We alsoexclude any configuration whose expectation value of H Ω is larger than a cut-off energy, E c . The value of E c is abit arbitrary, and it is taken fairly large compared to theexpected lowest adiabatic channel energy. The actual ba-sis size is around 250. Note that the basis functions allhave h ρ i = R but they have different σ values within σ ≦ σ ≦ .
5, case B those with σ ≦ . rms (fm)-20-15-10-505101520 E ( M e V ) H Ω T Ω V V V C
14 14.5 15 15.5R rms (fm)-3-2-101234 E ( M e V ) H Ω T Ω V V V C FIG. 3: (Color online) (a) Minimum expectation valueof H Ω calculated by a single configuration, Φ uAK L =0 M =0 ,as a function of R rms . The minimum is searched for byvarying K as well as u and A that are discretized on themeshes with M θ = 30, M φ = 21, and M ζ = 45. Thecontributions of T Ω , V , V , and V C to the minimumenergy are also drawn. (b) An enlarged figure of (a) in R rms = 14 . − . σ value of Φ JMπRν and if it is not larger than σ that characterizes each case, we accept that Φ JMπRν asa solution, otherwise it is discarded. Figure 4 plots theadiabatic channel energies in each case at four R rms radii.The solution of the constrained equations, (18) and (19),appears to be obtained stably. With the increase of thebasis size from case A to case C, the density of the adi-abatic channel energies considerably increases. Note thedifferent energy scale in Fig. 4(a) to 4(d).It is interesting to compare the present adiabatic chan-nel energies with those of Ref. [11]. The latter uses ba-sis functions quite different from ours: At each ρ , thechannel wave function is expanded in terms of a com-bination of the product of the Wigner D function andfifth-order basis splines for the two hyperangles. It in-cludes no ρ -dependence. In contrast to this, our chan-nel wave function has finite ρ -dependence, and there-fore receives influence from the adiabatic Hamiltonian0 U ( M e V ) (a)A B C -100102030405060 U ( M e V ) (b)A B C U ( M e V ) (c)A B C U ( M e V ) (d)A B C FIG. 4: (a) Adiabatic channel energies U at R rms = 1 . σ ≦ . σ ≦ .
75 in case B, while all the basis functions selected are allowed in case C. Seethe text for detail. Figures (b), (c), and (d) are the same as (a) but for R rms = 2 . , .
0, and 14 . ρ values. Thus both energies need not be nec-essarily the same but a comparison may indicate char-acteristics of different types of calculations. Three low-est adiabatic channel energies in MeV obtained in [11]are − . , . , . R rms = 1 . − . , . , . . , . , .
38 at 5.0 fm, and 0 . , . , . − . , − . − . − . , . , .
12 at 2.5 fm,1 . , . , .
68 at 5.0 fm, and 2 . , . , .
52 at 14.5 fm.The energy spacing of our calculation is much narrowerthan that of Ref. [11]. Which of the two calculationsgives lower value for the lowest channel energy seems todepend on R rms . At R rms = 2 . − .
75 MeV. Since their σ values are larger than 1, theyare not drawn in Fig. 4(b). Note, however, that the high-est one among the four is predicted to be − .
09 MeV with σ = 1 . − . σ valuesare larger than 1. One of them is located at 1.79 MeVand has σ = 1 .
08. It is still considerably higher than thelowest energy, 0 .
46 MeV, of Ref. [11].
V. SUMMARY
We have formulated hyperspherical calculations usingthe flexibility of the correlated Gaussians. Differentlyfrom conventional hyperspherical methods, the channelwave function and the adiabatic channel energy are de-fined by solving the hyperradius-constrained eigenvalueequation of the adiabatic Hamiltonian. This approachenables us to perform standard configuration interactioncalculations.This work takes a non-conventional venue by allowingthe spread of the value of the hyperradius for a givenbasis function. While in previous hyperspherical calcula-tions, e.g. in [11], the basis functions belong to a givenhyperradius, in the present work “the basis functions arelocalized”, that is, the hyperradii of the basis functionsreside in a narrow region around a predefined hyperra-dius. This approach has its advantages and disadvan-tages. A slight disadvantage is that one can not definea sharp hyperradius, so that the direct comparison toconventional calculations is not simple, whereas an ad-vantage is that the basis functions directly couple theneighboring regions which may help to resolve compli-cated dynamical processes. A further advantage is theeasier access to larger systems.The present formulation is expected to have many ap-1plications. As an example, we just mention one prob-lem, the fragmentation or decay of a nucleus into several α particles at large distances, as discussed in [43–45].The approaches employed there have limitations in tak-ing into account important effects such as couplings withother configurations, the angular momentum dependenceof the adiabatic potential, and the removal of spuriouscenter-of-mass excitations. Since the issue is exactly con-cerned with how the system evolves as it expands, it isworthwhile attempting at resolving those problems in thehyperspherical approach. ACKNOWLEDGMENTS
We are grateful to H. Suno for his interest in this workand for sending us his calculated results. We thank D.Blume, Q. Guan, T. Morishita, and I. Shimamura foruseful discussions at the early stage of the work.
Appendix A: Solving a constrained eigenvalueproblem
The aim of this appendix is to solve a problem of ob-taining the eigenvalues and corresponding eigenfunctionsof a Hermitian operator H with a constraint. It is for-mulated as follows: Let Q be a positive-definite Hermi-tian operator, and let ( φ , φ , . . . , φ n ) be a given set ofnormalized, independent basis functions. Obtain, in thespace spanned by the set, as many Φ’s possible that makethe expectation value of H h Φ | H | Φ ih Φ | Φ i (A1)stationary under the constraint h Φ | Q | Φ ih Φ | Φ i = q, (A2)where q is a positive constant. The condition h φ i | Q | φ i i = q is assumed for each φ i in the text. It may not be ab-solutely necessary, however, although the number of so-lutions may depend on how many basis functions satisfythe condition.This type of problem appears in several cases. See, forexample, Ref. [46] for the optimization of Φ and Ref. [47]for the determination of Φ free from some configurations.The present problem has distinct differences from thosecases in that the available configuration space is presetand several solutions are requested if possible.We construct an orthonormal set, ( ψ , ψ , . . . , ψ n ),that makes Q diagonal. To do this, we first diagonal-ize the overlap matrix ( h φ i | φ j i ): n X j =1 h φ i | φ j i u ( k ) j = b k u ( k ) i , (A3) where P ni =1 u ( k ) i u ( l ) i = δ kl . The basis set u k defined by u k = 1 √ b k n X i =1 u ( k ) i φ i ( k = 1 , . . . , n ) (A4)is orthonormal, h u k | u l i = δ kl . Next, diagonalizing Q inthe set u k , n X j =1 h u i | Q | u j i ψ ( k ) j = q k ψ ( k ) i , (A5)with P ni =1 ψ ( k ) i ψ ( l ) i = δ kl , we construct the set ψ k as ψ k = n X i =1 ψ ( k ) i u i , (A6)which has the desired property, h ψ k | ψ l i = δ kl , and h ψ k | Q | ψ l i = q k δ kl . It is easy to express ψ i in terms ofthe original set φ i ’s.We attempt to obtain Φ’s step by step. Defining aHermitian operator H ′ with a Lagrange multiplier λ , H ′ ( λ ) = H − λ ( Q − q ) , (A7)we solve the eigenvalue problem, H ′ ( λ )Φ( λ ) = E ′ ( λ )Φ( λ ) , (A8)using the set ψ i , and calculate the expectation value, F ( λ ) = h Φ( λ ) | Q − q | Φ( λ ) i = h Φ( λ ) | Q | Φ( λ ) i − q, (A9)where Φ( λ ) is normalized. Focusing always on the lowest-energy solution for any λ , we vary λ to find a zero of F ( λ ): F ( λ ) = 0. Then Φ( λ ) satisfies the constraint (A2) andthat is the solution to be found: Φ = Φ( λ ) with theenergy E = E ′ ( λ ).To determine the next solution, we define a configu-ration space of dimension ( n −
1) by removing Φ fromthe set ( ψ , ψ , . . . , ψ n ), and follow the above procedureto find a successful solution Φ . This process continuesuntil no new solution is found. Clearly Φ i ’s determinedin this way are orthogonal to each other.We can show that Φ i and Φ j for i = j have no couplingmatrix element of H if λ i = λ j . Since both functions areorthogonal, the matrix element of H reduces to that of Q as follows: h Φ j | H | Φ i i = h Φ j | H ′ ( λ i ) + λ i ( Q − q ) | Φ i i = E ′ ( λ i ) h Φ j | Φ i i + λ i h Φ j | Q − q | Φ i i = λ i h Φ j | Q | Φ i i . (A10)Because of h Φ j | H | Φ i i = h Φ i | H | Φ j i ∗ = λ j h Φ i | Q | Φ j i ∗ = λ j h Φ j | Q | Φ i i , it follows that( λ i − λ j ) h Φ j | Q | Φ i i = 0 . (A11)If λ i = λ j , h Φ j | Q | Φ i i vanishes and consequently h Φ j | H | Φ i i must vanish.If λ i and λ j are accidentally equal, the above argumentdoes not apply and it is not clear whether or not H hasthe coupling matrix element.2 Appendix B: Examples of correlated-Gaussianmatrix element
We show the examples of Eq. (30) that appear fre-quently. For O ( x ), let us consider the following terms:(i) e wx , (ii) ( e wx )( f w ′ x ) , (iii) e x C x , (iv) ( e wx )( e x C x ) , (v) ( e x C x )( e x C ′ x ) . (B1)Here w and w ′ are column vectors of dimension ( N − C and C ′ are ( N − × ( N −
1) symmetric matrices, and they areall independent of x . The integral of Eq. (30) is easilyobtained. Corresponding to (i)-(v) classes, P O reads(i) 12 e w B − v , (ii) 12 e w B − w ′ + 14 e w B − v f w ′ B − v , (iii) 32 Tr B − C + 14 e v B − CB − v , (iv) 12 e w B − CB − v + 34 ( e w B − v ) Tr B − C + 18 ( e w B − v )( e v B − CB − v ) , (v) 32 Tr B − CB − C ′ + 12 e v B − CB − C ′ B − v + (cid:16)
32 Tr B − C + 14 e v B − CB − v (cid:17) × (cid:16)
32 Tr B − C ′ + 14 e v B − C ′ B − v (cid:17) . (B2) As an example, we show how to obtain the matrix ele-ment of e x C x belonging to class (iii). The corresponding P O reads 32 Tr B − C + 14 e v G v (B3)with G = B − CB − , and it comprises four terms: T = 32 Tr B − C, T = 14 e uG u,T = 14 e u ′ G u ′ , T = 12 e uG u ′ , (B4)where T kk ′ l is an abbreviation of T kk ′ l ( u ′ A ′ , uA ). Equa-tion (46) immediately gives us the following result: h f u ′ A ′ K ′ LM | e x C x | f uAKLM i = h f u ′ A ′ K ′ LM | f uAKLM i h
32 Tr B − C + e uG u e uB − u F KK ′ L ( z )+ e u ′ G u ′ e u ′ B − u ′ F KK ′ L ( z ) + e uG u ′ e uB − u ′ F KK ′ L ( z ) i . (B5)The matrix element of ( e x C x ) is obtained similarly: h f u ′ A ′ K ′ LM | ( e x C x ) | f uAKLM i = h f u ′ A ′ K ′ LM | f uAKLM i h
32 Tr G C + 94 (Tr B − C ) + (2 e uG u + 3Tr B − C e uG u ) 1 e uB − u F KK ′ L ( z )+ (2 e u ′ G u ′ + 3Tr B − C e u ′ G u ′ ) 1 e u ′ B − u ′ F KK ′ L ( z ) + (2 e uG u ′ + 3Tr B − C e uG u ′ ) 1 e uB − u ′ F KK ′ L ( z )+ ( e uG u ) ( e uB − u ) F KK ′ L ( z ) + ( e u ′ G u ′ ) ( e u ′ B − u ′ ) F KK ′ L ( z ) + ( e uG u ′ ) ( e uB − u ′ ) F KK ′ L ( z )+ 2 ( e uG u )( e u ′ G u ′ )( e uB − u )( e u ′ B − u ′ ) F KK ′ L ( z ) + 2 ( e uG u )( e uG u ′ )( e uB − u )( e uB − u ′ ) F KK ′ L ( z ) + 2 ( e u ′ G u ′ )( e u ′ G u )( e u ′ B − u ′ )( e u ′ B − u ) F KK ′ L ( z ) i , (B6)where G = B − CB − CB − . [1] M. V. Zhukov, B. V. Danilin, D. V. Fedorov, J. M. Bang,I. J. Thompson, and J. S. Vaagen, Phys. Rep. , 151 (1993).[2] C. D. Lin, Phys. Rep. , 1 (1995). [3] R. Krivec, Few-Body Syst. , 199 (1998).[4] E. Nielsen, D. V. Fedorov, A. S. Jensen, and E. Garrido,Phys. Rep. , 373 (2001).[5] C. H. Greene, P. Giannakeas, and J. P´erez-R´ıos, Rev.Mod. Phys. , 035006 (2017).[6] I. J. Thompson, B. V. Danilin, V. D. Efros, J. S. Vaagen,J. M. Bang, and M. V. Zhukov, Phys. Rev. C , 024318(2000).[7] J. Macek, J. Phys. B , 831 (1968).[8] A. A. Kvitsinsky and V. V. Kostrykin, J. Math. Phys. , 2802 (1991).[9] N. B. Nguyen, F. M. Nunes, and I. J. Thompson, Phys.Rev. C , 054615 (2013).[10] S. Ishikawa, Phys. Rev. C , 055804 (2013).[11] H. Suno, Y. Suzuki, and P. Descouvemont, Phys. Rev. C , 054607 (2016).[12] H. Suno, Y. Suzuki, and P. Descouvemont, Phys. Rev. C , 014004 (2015).[13] N. Barnea, W. Leidemann, and G. Orlandini, Phys. Rev.C , 054001 (2000).[14] N. Barnea, W. Leidemann, and G. Orlandini, Phys. Rev.C , 054003 (2003).[15] N. Barnea, W. Leidemann, and G. Orlandini, Phys. Rev.C , 064001 (2010).[16] S. Bacca, N. Barnea, and A. Schwenk, Phys. Rev. C ,034321 (2012).[17] N. K. Timofeyuk, Phys. Rev. C , 064306 (2002).[18] N. K. Timofeyuk, Phys. Rev. C , 054314 (2008).[19] M. Gattobigio, A. Kievsky, and M. Viviani, Phys. Rev.C , 024001 (2011).[20] J. von Stecher and C. H. Greene, Phys. Rev. A , 022504(2009).[21] S. T. Rittenhouse, J. von Stecher, J. P. D’Incao, N. P.Mehta, and C. H. Greene, J. Phys. B: At. Mol. Opt. Phys. , 172001 (2011).[22] D. Rakshit and D. Blume, Phys. Rev. A , 062513(2012).[23] K. M. Daily and C. H. Greene, Phys. Rev. A , 012503(2014).[24] S. F. Boys, Proc. R. Soc. London, Ser. A , 402 (1960).[25] K. Singer, Proc. R. Soc. London, Ser. A , 412 (1960).[26] K. Varga and Y. Suzuki, Phys. Rev. C , 2885 (1995).[27] Y. Suzuki and K. Varga, Stochastic Variational Ap- proach to Quantum-Mechanical Few-Body Problems , Lec-ture Notes in Physics, m , (Springer, Berlin, 1998).[28] Y. Suzuki, J. Usukura and K. Varga, J. Phys. B: At. Mol.Opt. Phys. , 31 (1998).[29] K. Varga, Y. Suzuki, and R. G. Lovas, Nucl. Phys. A571 , 447 (1994).[30] J. Mitroy, S. Bubin, W. Horiuchi, Y. Suzuki, L. Adamow-icz, W. Cencek, K. Szalewicz, J. Komasa, D. Blume, andK. Varga, Rev. Mod. Phys. , 693 (2013).[31] W. Horiuchi and Y. Suzuki, Phys. Rev. C , 011304(2014).[32] D. Mikami, W. Horiuchi, and Y. Suzuki, Phys. Rev. C , 064302 (2014).[33] Y. Suzuki, Prog. Theor. Exp. Phys. , 043D05.[34] O. I. Tolstikhin, S. Watanabe, and M. Matsuzawa, J.Phys. B: At. Mol. Opt. Phys. , L 389, (1996).[35] H. Suno, J. Chem. Phys. , 064318 (2011).[36] Y. Suzuki, W. Horiuchi, M. Orabi and K. Arai, Few-BodySyst. , 33 (2008).[37] S. Aoyama, K. Arai, Y. Suzuki, P. Descouvemont, andD. Baye, Few-Body Syst. , 97 (2012).[38] Y. Suzuki and W. Horiuchi, Phys. Rev. C , 044320(2017).[39] A. Erd´elyi, W. Magnus, F. Oberhettinger, and F. G.Tricomi, Higher Transcendental Functions , Vol. I,(MacGraw-Hill, New York, 1953).[40] M. Abramowitz and I. A. Stegun,