JM-ECS: A hybrid method combining the J -matrix and ECS methods for scattering calculations
Y. Bidasyuk, W. Vanroose, J. Broeckhove, F. Arickx, V. Vasilevsky
JJM-ECS: A hybrid method combining the J -matrix and ECS methods for scatteringcalculations Y. Bidasyuk ∗ Departement Wiskunde-Informatica, Universiteit Antwerpen, Antwerpen, Belgium andBogolyubov Institute for Theoretical Physics, Kyiv, Ukraine
W. Vanroose, † J. Broeckhove, ‡ and F. Arickx § Departement Wiskunde-Informatica, Universiteit Antwerpen, Antwerpen, Belgium
V. Vasilevsky
Bogolyubov Institute for Theoretical Physics, Kyiv, Ukraine (Dated: October 23, 2018)The paper proposes a hybrid method for calculating scattering processes. It combines the J -matrixmethod with exterior complex scaling and an absorbing boundary condition. The wave function isrepresented as a finite sum of oscillator eigenstates in the inner region and it is discretized on a gridin the outer region. The method is validated for a one- and two-dimensional model with partialwave equations, and a calculation of p -shell nuclear scattering with semi-realistic interactions. PACS numbers: 02.60.Cb, 21.60.Gx, 25.55.Ci, 03.65.Nk
I. INTRODUCTION
Numerical simulations of breakup reactions are amongthe most complicated problems in the theoretical physicsof nuclear reactions. The main complication arises fromthe necessity of modeling the asymptotic behavior of thewave function in a many-body continuum. To simplifythe asymptotic behavior, the wave function is usuallyexpressed in terms of hyperspherical harmonics [1–3].But, in breakup problems, the hyperpotentials decay veryslowly and convergence with respect to the hyper-angularmomentum is also very slow. This results in very largesystems of equations, even for the simplest breakup prob-lems.The Complex Scaling approach [4, 5] can be used toavoid the problems associated with the asymptotic be-havior of the many-body wave function. However, thismethod is suitable only to analyze resonances in the sys-tem and it does not provide any breakup information,such as cross sections, that can be compared with exper-iments.In the numerical solution of breakup problems de-scribed by the Schr¨odinger equation it is important todescribe each of the breakup channels in a correct way.Once the asymptotic form in each of the channels is un-derstood, it is possible to formulate a set of equationswhose solution yields both the wave function in the in-teraction region, and the asymptotic amplitudes in eachbreakup channel. This is the approach taken both by the R -matrix [6] and the J -matrix method [7]. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected]
However, for some of the breakup channels the asymp-totic wave function is only known at very large distancefrom the target. As a consequence very large numeri-cal domains are needed to cover the interaction and thenear-field regions. In other cases, such as the three bodybreakup in the J -matrix representation, the explicit ex-pression of the asymptotic wave function is not known.The need for an explicit asymptotic wave fuction canbe avoided by the introduction of absorbing boundaryconditions. These enforce an “outgoing wave” bound-ary condition on the numerical solution, and have beenused successfully in atomic and molecular physics [8] andacoustic and electromagnetic scattering problems [9]. In-stead of solving in a single system both the wave functionin the interaction region and the reaction rates, thesemethods calculate the cross section as a post-processingstep. First the equations are solved with the absorbingboundary conditions and then, in a second step, the am-plitudes are extracted from the numerical solution.Absorbing boundary conditions are easy to implementin a grid based discretization of a partial differential equa-tion. They have been implemented in calculations basedon finite difference, B-splines, finite elements etc [8, 10].In nuclear physics, however, one prefers to use a L -basis, such as the oscillator eigenstates. Such basis iswell-adapted to the fully microscopic description of thecompound nucleus. This paper reports on the initial ef-forts to introduce this approach in the context of nuclearfew-body systems. We propose a hybrid approach withan L representation of the wave function in the inter-action region, and with a discretized grid representationin the outer region. Our aim is to combine the strengthsand benefits of both approaches.The paper is structured as follows. In section II, af-ter a short review of the ECS and J -matrix method, weintroduce the hybrid method where the wave function isrepresented in the inner region by oscillator states and in a r X i v : . [ nu c l - t h ] N ov the outer region by a grid. We also discuss how the ob-servables are extracted from the numerical representationof the scattered wave. In section III we present numeri-cal results that validate the proposed method with modelproblems for one- and two-dimensional partial wave equa-tions. In section IV we show results for nuclear p -shellscattering. We use m = 1 and (cid:126) = 1 throughout thepaper, except where mentioned. II. THE HYBRID J -MATRIX WITH ECSMETHODA. Exterior Complex Scaling as an outgoing waveboundary condition Exterior Complex Scaling (ECS) was introduced by B.Simon [11] and has been initially used for the calculationof resonance positions and widths [5]. It has also foundwidespread application by providing an absorbing bound-ary condition in atomic and molecular breakup problemswith charged particles [12, 13]. A review is given in [8].In these problems the outgoing wave is very complicated,and its form depends on the effective interaction of theoutgoing particle. It is determined by the position of theother charged particles involved in the breakup problem.The derivation of ECS is based on an analytical con-tinuation in the complex plane. It differs from the well-known
Complex Scaling (CS) by starting the scaling pro-cedure well into the asymptotic region.We briefly explain why it leads to outgoing waveboundary conditions and apply these to the Schr¨odingerequation.
1. Enforcing outgoing wave boundary conditions
Consider a one-dimensional Helmholtz equation (cid:18) − d dρ − k (cid:19) u ( ρ ) = f ( ρ ) (1)on a domain [0 , L ], with constant wave number k > ∈ R , homogeneous Dirichlet conditions on the left bound-ary, and outgoing wave boundary conditions on the rightboundary. The right hand side, f ( ρ ), is such that it iszero outside the interval [0 , a ] with a < L .Let us focus this discussion on the right boundary in L . For a < ρ < L the equation is a second order homoge-neous differential equation with constant coefficients. Inthis region the general solution can then be written as u ( ρ ) = Ae ikρ + Be − ikρ , (2)a linear combination of two fundamental solutions. En-forcing outgoing wave boundary conditions in ρ = L means that coefficient B should be zero. This can be realized by the following mixed type boundary conditionon the solution u (cid:48) ( L ) = i k u ( L ) . (3)Note that this boundary condition requires the explicitknowledge of the wave number k .ECS is an alternative way to enforce the same outgo-ing wave condition. When we analytically continue theequation (1) to complex ρ ∈ C , the general solution ofthe equation, where it is homogeneous, remains a linearcombination of the same fundamental modes as in (2).When we impose, instead of the condition (3), a homo-geneous Dirichlet boundary condition in the point L (cid:48) ∈ C that lies inside the region where the equation is homoge-neous, we find that u ( L (cid:48) ) = 0 leads to B/A = exp(2 ikL (cid:48) )or | B || A | = exp( − k Im( L (cid:48) )) . (4)So if the point L (cid:48) ∈ C is chosen such that Im( L (cid:48) ) > k Im( L (cid:48) ) (cid:29) | B | (cid:28) | A | . As a result, we have,effectively, enforced outgoing wave boundary conditionsin the point L (cid:48) .The linear combination of fundamental modes withcoefficients A and B describes the solution everywherewhere the equation (1) is homogeneous. The equation ishomogeneous between the point L and L (cid:48) . So the coef-ficient B of the solution is still much smaller than A inthe point L and we can conclude that we also have anoutgoing wave boundary condition in L .It is important to note that, in contrast to the mixedboundary condition (3), enforcing Dirichlet boundaryconditions in L (cid:48) does not require the knowledge of k . Thismakes it possible to describe both inelastic processes andoutgoing boundary conditions in higher dimensions.
2. ECS contour
We implement this by extending the interval [0 , L ] witha contour that connects L with L (cid:48) . The point L (cid:48) is typi-cally chosen such that | L | < | L (cid:48) | and k | L (cid:48) | (cid:29) , Re( L (cid:48) )]. z ( ρ ) = (cid:26) ρ if ρ ≤ L.ρ + if ( ρ ) if L < ρ <
Re( L (cid:48) ) . (5)where f is an increasing function — linear, quadratic,. . . — with f ( Re ( L (cid:48) )) = Im ( L (cid:48) ) and f ( L ) = 0.In particular, ECS considers a linear function, and thecoordinate transformation is therefore written as z ( ρ ) = (cid:40) ρ for ρ ≤ L.L + ( ρ − L )e i θ for ρ > L. (6)
3. Application to the Schr¨odinger equation
The application of these outgoing wave boundary con-ditions to the Schr¨odinger equation is straightforward.Consider, as an example, the two-particle model prob-lem for which the radial equation is( T + V ( ρ ) + l ( l + 1)2 ρ − E ) ψ sc ( ρ ) = − V ˆ j l ( kρ ) , (7)where T = ( − / d /dρ denotes the kinetic energy op-erator, k = √ E is the wave number, l = 0 , , , . . . the angular momentum, V the potential that only de-pends on the radial coordinate ρ , and ˆ j l ( kρ ) the Ricatti-Bessel function, which is the radial part of the incom-ing wave. The radial solution of this equation is ψ ( ρ ) =ˆ j l ( ρ ) + ψ sc ( ρ ), a sum of the incoming wave, ˆ j l ( ρ ), and thescattered wave, ψ sc ( ρ ). The latter should fit the outgoingwave boundary conditions.For problems with short range potentials theSchr¨odinger equation reduces to the Helmholtz equationoutside the range of the potentials with a wave number k .We can then apply the ECS transformation in the regionwhere it reduces to a Helmholtz equation.
4. Discretization
ECS has been implemented in finite differences, B-splines and spectral elements [8]. In this section wediscretize the differential operator that appears in theSchr¨odinger equation with finite differences using theShortley-Weller formula for non-uniform grids [14]. Itenables the discretization of the operator on the complexcontour and uses complex valued mesh widths. Such amesh mesh is illustrated on Fig. 1
Re(z) I m ( z ) h he i θ L Re(L’)
FIG. 1: The choice of grid points for the finite difference rep-resentation of the Helmholtz equation on a exterior complexscaled (ECS) domain. The original problem was stated on[0 , L ] with outgoing wave boundary conditions in L . The do-main is extended with [ L, L (cid:48) ] where a complex grid distance h e iθ is used. In L (cid:48) homogeneous Dirichlet boundary condi-tions are enforced. Consider the differential operator d /dz ( ρ ) on thecontour defined by the coordinate transformation (6).We define a uniform grid( z i ) ≤ i ≤ n on [0 , L ] with z = 0 and z n = L and mesh width h = 1 /n ∈ R ,and a second uniform grid on the complex contour( z i ) n ≤ i ≤ n + m on [ L, L (cid:48) ]with z n + m = L (cid:48) and complex mesh width h γ = ( L (cid:48) − L ) /m . The union of these two grids is the grid( z i ) ≤ i ≤ n + m on [0 , L ] ∪ [ L, L (cid:48) ] (8)in the entire ECS domain. To approximate the secondderivative in each grid point z i we use d u ( z i ) dz ≈ h i − + h i × (cid:18) h i − u ( z i − ) − (cid:18) h i − + 1 h i (cid:19) u ( z i ) + 1 h i u ( z i +1 ) (cid:19) , (9)where h i − and h i are the left and right mesh widths re-spectively, which may be complex valued. The formulareduces to regular second order central differences when h i − = h i , i.e., in the interior real region [0 , L ], and in theinterior of the complex contour [ L, L (cid:48) ], since the scalingfunction f is taken to be linear. The only exception isthe point z n where we lose at most an order of accuracy.However, with ample discretization steps, the overall ac-curacy is anticipated to match up to second order. Notethat a higher order discretization at the hinge is proposedin [15] to maintain a second order accuracy throughoutthe domain.The Hamiltonian of the one-dimensional model (7) isthen a ( n + m ) by ( n + m ) matrix H l ∈ C n + m . H l = T + l ( l + 1)2 diag (1 /z i ) + diag ( V ( z i )) , (10)where T is the finite difference representation of the sec-ond derivative on the complex contour, and diag is a di-agonal matrix with the potential evaluated in each gridpoint z i .For two-dimensional problems, the Hamiltonian is dis-cretized starting from the one-dimensional discretizedHamiltonian using the Kronecker product ⊗ H D = H l ⊗ I + I ⊗ H l + diag ( V ( z i , z j )) , (11)where I is a ( n + m ) by ( n + m ) unit matrix, and diag ( V )is the two body potential evaluated at the grid points.The matrix H D is now a complex valued ( n + m ) by( n + m ) matrix.The numerical solution of the PDE (7) is found bysolving the linear system( H l − E ) x = b , (12)where b is a numerical representation of the right handside of the equation.Some examples of one- and two-dimensional wave func-tions on the ECS domain are presented in Fig. 2 and 3. Ψ L ρ FIG. 2: The real and imaginary part of the scattered wavefor an Helmholtz problem (1) on the ECS grid. Because of thehomogeneous Dirichlet boundary condition the wave is onlyoutgoing in L .FIG. 3: The real part of a 2D partial wave of a three-particle s -wave scattering problem on an ECS domain, (49). Theproblem has a homogeneous Dirichlet boundary condition onthe domain [0 , L (cid:48) ] and because of the ECS this leads to aproblem on the domain [0 , L ] with homogeneous Dirichletboundary condition on the south and west boundary, andoutgoing wave boundary conditions on the north and eastboundary where x = L or y = L . We show the numericalsolution of the problem discretized with finite differences.
5. Extraction of the observables
Once we have the numerical solution of the equation onan ECS domain, we need to extract physical observablessuch as cross sections.In one dimension the amplitude A of the outgoing wavecan be extracted from the numerical solution with thehelp of the Wronskian. Indeed, outside the range of the potential V the solution can be written as a linear com-bination ψ sc = A ˆ h + l ( kr )+ B ˆ h − l ( kr ), where ˆ h ± l are the in-and outgoing Ricatti-Hankel functions. The coefficient A is then A = W (cid:16) ψ sc ( ρ ) , ˆ h − l ( kρ ) (cid:17) /W (cid:16) ˆ h + l ( kρ ) , ˆ h − l ( kρ ) (cid:17) , (13)where ρ ∈ [ a, L ] is outside the range of the potential butstill on the real part of the ECS domain. The Wronskianis calculated as W ( u, v ) = u (cid:48) v − v (cid:48) u . Since only the valuesof ψ sc are known on the grid points, its first derivative isapproximated by central differences.In two dimensions the problem is more complicated.The solution inside the numerical box has not reached itsfar-field form and a projection on the asymptotic stateis inaccurate since the single particle and double particlebreakup wave functions live in the same regions of space,especial near the edges of the domain. A procedure toextract the observables with the help of surface integralswas proposed by McCurdy, Horner and Rescigno in [16].The amplitude can be written as a surface integral f l ,l ( k , f ) = 12 (cid:90) S ( τ k ,l ( ρ ) τ k ,l ( ρ ) ∇ ψ sc ( ρ , ρ ) − ψ sc ( ρ , ρ ) ∇ ( τ k ,l ( ρ ) τ k ,l ( ρ ))) · d S , (14)where τ k ,l ( ρ ) is the solution of the radial equation (cid:18) T + l ( l + 1)2 ρ + V − k (cid:19) τ k ,l ( ρ ) = 0 . (15)In a similar way τ k ,l ( ρ ) fits the same equation with T , V , l and k such that k + k = 2 E . The contour of thesurface integral, (14), should lie in the region where theSchr¨odinger equation is homogeneous. In the literaturethe functions τ k,l are often denoted as φ k,l . In this paperwe have chosen a different notation to avoid confusionwith the basis functions of the oscillator representationthat will be considered in the next section.The single differential cross section (SDCS) is then dσdE = 8 π k k k | f ( k , k ) | , (16)where k is the momentum of the incoming wave [16].From these amplitude formulas it is also clear why weuse exterior complex scaling rather than complex scal-ing. In the latter the complete domain is rotated intothe complex plane as ρ → ρe iθ and homogeneous Dirich-let boundary conditions are enforced at the end of thegrid. However, after this transformation it is very hardto extract scattering information. In contrast, the solepurpose of ECS is to provide the correct outgoing bound-ary conditions. It leaves the solution on the real part ofthe grid unchanged and the scattering cross sections canbe extracted. B. The J -Matrix Method In the original J -Matrix method for quantum scatter-ing [7, 17] the wave function is represented in terms ofsome L -basis set that leads to a tri-diagonal structureof the Hamiltonian matrix for the free-particle problem(Jacobi shape of the matrix). This tri-diagonal structureallows the use of the complete basis set without actu-ally working with matrices of infinite size. A review ofthe applications of the J -matrix method can be found in[18].A popular L -basis set for calculations with short-range potentials is the set of eigenstates of the radialharmonic oscillator. The basis states are then generalizedLaguerre polynomials multiplied by a weight function φ i,l ( ρ ) = ( − i N i,l b − / (cid:16) ρb (cid:17) l × ρ exp (cid:18) − ρ b (cid:19) L l +1 / i (cid:18) ρ b (cid:19) (17)with a normalization N i,l = (cid:115) i !Γ( i + l + 3 / , (18)where i = 0 , , . . . and oscillator length b = (cid:112) (cid:126) /mω that is related to the oscillator frequency ω . Note thatwe have incorporated the weight ρ of integration comingfrom the radial coordinates in the function (17) such that (cid:82) ∞ φ i,l ( ρ ) φ j,l ( ρ ) dρ = δ ij Each basis function has a classical turning point R i,l = b √ i + 2 l + 3 . (19)The oscillator basis set is complete over L and the solu-tions of (7) can always be represented as a linear combi-nation of all oscillator states ψ l ( ρ ) = ∞ (cid:88) i =0 c i,l φ i,l ( ρ ) . (20)This representation reduces the Schr¨odinger equation toan infinite system of linear equations for c i,l .As already mentioned above, the kinetic energy oper-ator T (the free particle Hamiltonian) is a tri-diagonalmatrix in the J -matrix method. For the oscillator basisthe non-zero elements are J i,j = (cid:0) i + l + (cid:1) (cid:126) ω for j = i, − (cid:113) i (cid:0) i + l + (cid:1) (cid:126) ω for j = i − , − (cid:113) ( i + 1) (cid:0) i + l + (cid:1) (cid:126) ω for j = i + 1 . (21)However, the potential energy matrix in this represen-tation will be dense. Thus for an accurate treatment ofthe problem we need to deal with an infinitely-sized denseHamiltonian matrix. As long as we deal with short-range potentials only, this dense potential matrix can be safelytruncated to a finite matrix.This relies on the fact that highly excited oscillatorstates (with large index i ) oscillate rapidly between theorigin and the corresponding classical turning point R i,l .So, if the range of the potential is less then R i,l , its matrixelements will average to negligibly small value because ofannihilating oscillatory contributions. The potential ma-trix can then be truncated at some i = N determined bythe desired accuracy. Beyond this point, the Hamiltonianmatrix is approximated by the tri-diagonal (asymptotic)form. We therefore refer to the dense part of the Hamilto-nian matrix corresponding to i ≤ N as the “interactionregion”, and to the tridiagonal part for i > N as the“asymptotic region”.In the asymptotic region, the tri-diagonal structure ofthe matrix leads to a simple three-term recurrence rela-tion for the oscillator expansion coefficients of the solu-tion: J i,i − c i − + ( J i,i − E ) c i + J i,i +1 c i +1 = 0 ∀ i ≥ N (22)Since this is a second order recurrence relation the { c i } can be obtained as a linear combination of two linearlyindependent fundamental solutions. In the J -matrixmethod the regular b i,l and the irregular n i,l are cho-sen as fundamental solutions that can be easily obtainedfrom the explicit form of the matrix J , eq. (21), c i,l = b i,l + t n i,l ∀ i ≥ N. (23)These fundamental solutions of the recurrence relationhave a corresponding coordinate-space solution ψ l ( ρ ) = B l ( ρ ) + t N l ( ρ ) , (24)where B l ( ρ ) = (cid:80) ∞ i =0 b i,l φ i,l ( ρ ), N l ( ρ ) = (cid:80) ∞ i =0 n i,l φ i,l ( ρ )are usually referred to as Bessel-like and Neumann-likerespectively, due to their asymptotic behavior. In thiscontext the coefficient t in this linear combination corre-sponds to the scattering t -matrix.The full solution is then ψ l ( ρ ) = χ l ( ρ ) + B l ( ρ ) + t N l ( ρ ) , (25)where the function χ l ( ρ ) = (cid:80) ∞ i =0 c i,l φ i,l ( ρ ) is nonzeroonly in the interaction region. This leads to c i,l = (cid:26) c i,l + b i,l + t n i,l when i < Nb i,l + t n i,l when i ≥ N . (26)With this form of the solution we reduce the set of un-knowns to { c i,l , t } and the linear system to be solved has N + 1 dimensions. Solving the linear system we obtainsimultaneously the wave function of the system, { c i,l } ,and the scattering information, t . C. Introduction of the JM-ECS method
For multi-particle scattering and reaction (breakup)problems it is very hard to obtain an explicit form ofthe asymptotic solutions, also in the oscillator represen-tation. It is then a natural approach to avoid such ex-plicit forms by introducing the ECS transformation inthe asymptotic JM region by enforcing outgoing waveboundary conditions.A formal introduction of the ECS within the oscillatorbasis representation is hard for several reasons. First, ap-plying the ECS coordinate transformation (6) to oscilla-tor functions destroys the orthogonality as well as the tri-diagonal structure of the kinetic energy matrix. Second,the convergence of the oscillator basis will be stronglyaffected. Most probably, the wave function within theabsorbing layer will be poorly reproduced even by a hugetruncated basis set. Finally, the coordinate transforma-tion makes an analytical calculation of matrix elementsdifficult.An alternative approach is to combine the grid and theoscillator representation and represent the wave functionin the asymptotic region with finite differences. This be-comes possible due to a specific property of the oscilla-tor basis. For highly excited oscillator states the oscilla-tor expansion coefficients are related to the values of thewave function on the grid of classical turning points [19]: c n,l = (cid:90) ∞ φ i,l ( ρ ) ψ l ( ρ ) dρ ≈ b (cid:113) /R n,l (cid:20) Ψ l ( R n,l ) + O (cid:18) R n,l (cid:19)(cid:21) . (27)This allows to couple the oscillator representation ofthe wave function with the coordinate-space representa-tion in the high- n region. A similar argument has beenused for building the so-called modified J -matrix method(MJM) [19].
1. A hybrid representation of the wave function
In the hybrid JM-ECS method we represent our one-dimensional wave function as a vector Ψ in C n + m , where Ψ = ( c , c , . . . , c n − ,ψ ( R n ) , ψ ( R n + h ) , . . . , ψ ( R n +( m − h )) . (28)The first n elements represent the wave function in the interaction region in the oscillator representation, whilethe remaining m elements represent the wave function inthe asymptotic region on an equidistant grid that startsat R n , the n -th classical turning point, and runs up to R n + ( m − h with a grid distance h = R n − R n − . (29)We assume that the matching point, which connects theoscillator to finite difference representation, corresponds to an large index n such that the asymptotic formula forthe expansion coefficient, (27), can be applied.Again, the kinetic energy operator in this hybrid rep-resentation is tridiagonal since it is tridiagonal bothin finite-differences and oscillator representations. Oneshould only be careful near the matching point betweenboth representations so that the asymptotic formula (27)is a sufficiently good approximation for representing thesolution.To obtain the kinetic energy in the final point of theoscillator representation, we use the tridiagonal kineticenergy formula (21). It involves a recurrence relationconnecting the three terms c n − , c n − and c n . The lat-ter, the coefficient c n , is unknown. Only ψ ( R n ) is avail-able. Using the asymptotic relation (27), however, wecan calculate the required matrix element as follows:( Jc ) n − = J n − ,n − c n − + J n − ,n − c n − + J n − ,n b (cid:112) /R n ψ l ( R n ) . To calculate the kinetic energy in the first point of thefinite difference grid, the second derivative of the wavefunction has to be known. To approximate the latter witha finite difference formula, one needs the wave functionin the grid points R n − , R n and R n + h . We again apply(27) to obtain ψ ( R n − ) in terms of c n − : ψ (cid:48)(cid:48) ( R n ) = c n − / ( b (cid:112) /R n − ) − ψ ( R n ) + ψ ( R n + h ) h . The coupling between both representations around thematching point is sketched in Fig. 4, together with theterms involved to determine the correct matching.
Oscillator Finite Differences Tc n − (cid:122) (cid:125)(cid:124) (cid:123) c n − c n − Ψ( R n ) Ψ( R n + h ) (cid:124) (cid:123)(cid:122) (cid:125) Ψ (cid:48)(cid:48) ( R n ) (cid:116)(cid:100) (cid:100) (cid:116) FIG. 4: This figure illustrates how the kinetic energy matrixelements are calculated in the last point R n − of the oscillatorrepresentation and in the first point R n of the finite differencerepresentation. To calculate T applied on a solution vectorwe need to translate the oscillator representation to the gridand vice versa. The structure of the matrix representation of the ki-netic operator around the matching point is then . . . . . . J n − ,n − J n − ,n − J n − ,n b (cid:112) /R n / ( h b (cid:112) /R n − ) − /h /h . . . . . . ... c n − ψ ( R n )... . (30)The representation of the potential operator in the hy-brid JM-ECS method is more complex. Similar to theJM method the potential matrix in the interaction re-gion, covered by the oscillator representation, is dense.Therefore, the full hybrid potential matrix will be densein the interaction region and diagonal in the asymptoticfinite-differences region. It is clear that the computa-tional complexity of the problem is determined by thesize of the interaction region because of the large numberof potential matrix elements that needs to be calculated.Shifting the matching point further outward will there-fore strongly affect the computation time in a negativeway. Increasing the finite-differences asymptotic part will have almost no effect on the latter because of the diago-nal potential.The introduction of the outgoing wave boundary con-dition is now straightforward by extending the finite dif-ference grid with an ECS contour as explained in sectionII A. We use m grid points to cover both the real partand the ECS part of the finite difference grid.Similar to the finite difference representation (11)one can easily construct a representation for a two-dimensional wave function. Instead of the vector form(28), the solution can be represented by a matrix Ψ ∈ C ( n + m ) × ( n + m ) that has the following structure Ψ = c . . . c n − d n . . . d n + m c . . . c n − d n . . . d n + m c n − . . . c n − n − d n − n . . .d n . . . d n n − ψ ( R n , R n ) . . . ψ ( R n , R n + m ) d n +1 0 . . . d n +1 n − ψ ( R n +1 , R n ) . . . ψ ( R n +1 , R n + m ) . . . . . .d n + m . . . d n + m n − ψ ( R n + m , R n ) . . . ψ ( R n + m , R n + m ) . (31)An example of this wave function is presented in Fig. 5.The spatial wave function can then, depending on theregion in the two-dimensional domain, be written as fol-lows. For ρ , ρ < R n one has ψ ( ρ , ρ ) = n − (cid:88) i,j =0 Ψ ij φ i ( ρ ) φ j ( ρ ) , (32)for ρ < R n and k ≥ n one writes ψ ( ρ , R l ) = n − (cid:88) i Ψ il φ i ( ρ ) (33)and for ρ < R n and l ≥ n , ψ ( R k , ρ ) = n − (cid:88) j Ψ kj φ j ( ρ ) (34)and, finally, for k, l ≥ nψ ( R k , R k ) = Ψ kl . (35) The Hamiltonian of a 2D problem is again constructedas a Kronecker product of the 1D Hamiltonians and atwo-body potential. The locations of non-zero matrixelements are the same as in the Kronecker product ofthe two two-particle potential matrices. This leads toa number of dense blocks distributed over a large sparsematrix. In this case the computational complexity is alsodetermined by the size of the finite-differences part, asit not only extends the diagonal, but also increases thenumber of dense blocks.
2. Extracting Scattering information from a wave functionin a hybrid representation
As already mentioned at the beginning of this section,the extraction of observables is not always straightfor-ward in the JM method. This mainly comes from thefact that in the original formulations one has to solve thescattering problem and define the scattering parameterssimultaneously.
FIG. 5: The real part of the 2D wave function of the Three-particle s-wave scattering problem of (49) in the hybrid ap-proach with 40 oscillator states and the oscillator parameter b = 0 .
8. The values of X and Y in the oscillator region arechosen as the corresponding classical turning points to showthe similar behavior of the wave function in different repre-sentations. Compare with Fig. (3). Note the jump in thewave function near the J-matrix and finite difference bound-ary. This jump clearly indicates the different regions in therepresentation. In the hybrid JM-ECS approach in one dimension, oncea solution is obtained, we still need to extract the ob-servables. Since the wave function is represented in theasymptotic region with a finite-difference representation,we can use a standard Wronskian technique to extractthe scattering amplitude (see section II A 5).Obtaining the scattering information from the solutionin two dimensions requires the surface integral of sectionII A 5. This approach still needs to be elaborated in thehybrid representation. We consider for this discussion acontour S of the surface integral (14) that is piecewiseparallel to one of the axes. This is illustrated in Fig. 6.We first integrate ρ from 0 to a value R k ∈ [ a, L ] whilefixing ρ = R k . The point R k is located on the finitedifference grid in the region where the equation is homo-geneous. Typically R k is chosen a few grid points beforethe end of the real part of the grid. This procedure leadsto a line integral parallel to the ρ axis. We denote it as I .The second part of the surface integral integrates ρ from R k to 0 while keeping ρ equal to R k . This is a lineintegral parallel to the ρ axis and is denoted as I .The amplitude, (14), now consists of two parts f l ,l ( k , k ) = I + I , (36) FIG. 6: There are four region in the hybrid representation.In region a) ρ , ρ < R n and the wave function is a sum overproduct of two oscillator functions. In b) and d) one onecoordinate, respectively ρ < R n and ρ < R n , is representedin the oscillator state. In region c) we use finite differencefor both coordinates. The amplitude is calculated with thehelp of a surface integral along the path that is indicated bythe solid line a distance R k away from the axes. Note thatbetween L and Re( L (cid:48) ) the grid is complex valued because ofECS. where I = 12 (cid:90) R k (cid:18) τ k ,l ( ρ ) τ k ,l ( R k ) ∂ψ sc ( ρ , R k ) ∂ρ − ψ sc ( ρ , R k ) τ k ,l ( ρ ) ∂τ k ,l ( R k ) ∂ρ (cid:19) dρ (37)and I = 12 (cid:90) R k (cid:18) τ k ,l ( R k ) τ k ,l ( ρ ) ∂ψ sc ( R k , ρ ) ∂ρ − ψ sc ( R k , ρ ) τ k ,l ( ρ ) ∂τ k ,l ( R k ) ∂ρ (cid:19) dρ . (38)This is in fact a surface integral of an in-product of twovector-valued functions. We can reverse the bounds on I without switching the sign of the argument, since thesurface vector also changes direction.It is important to note that I probes the 2D wavefunction ψ sc ( ρ , ρ ) where it is either represented as asum over oscillator states in the first coordinate and afinite difference in the second or as a finite difference inboth coordinates. The integral never probes the regionwhere both coordinates are represented in the oscillatorrepresentation (see region (a) in Fig. 6).To arrive at a practical expression for the amplitude,we split I into (cid:82) R n ( . . . ) + (cid:82) R k R n ( . . . ) into an integrationfrom 0 to R n , the part covered by the oscillator represen-tation, and an integration from R n to R k , covering thegrid.In order to calculate the integrals we require the so-lutions τ k ,l ( ρ ) and τ k ,l ( ρ ) of equation (15). We solvethese equations also in the hybrid representation. The so-lutions are vectors t ( k ) and t ( k ) ∈ C ( n + m ) representingthe numerical solution of τ k ,l ( ρ ) and τ k ,l ( ρ ), respec-tively, with k + k = E . For ρ ∈ [0 , R n ] the solution isa finite sum τ k ,l ( ρ ) = n − (cid:88) i t ( k ) i φ i ( ρ )and, similarly, the 2D wave function in region (b) of Fig.(6) equals ψ sc ( ρ , R k ) = (cid:80) n − i =0 Ψ ik φ i ( ρ ).The calculation of the amplitude requires the firstorder derivatives ∂τ k ,l ( R k ) /∂ρ and ∂ψ sc ( ρ , R k ) /∂ρ .These will be approximated by central differences from t ( k ) and Ψ . We define t (cid:48) k ( k ) = t ( k ) k +1 − t ( k ) k − h , (39)and similarly for t (cid:48) ( k ) k . We also define Ψ (cid:48) , ik = Ψ i +1 k − Ψ i − k h (40)and Ψ (cid:48) , ik = Ψ i k +1 − Ψ i k − h . (41)So the first term in the integral of I becomes (cid:90) R n (cid:18) − τ k ,l ( ρ ) τ k ,l ( R k ) ∂ψ sc ( ρ , R k ) ∂ρ + ψ sc ( ρ , R k ) τ k ,l ( ρ ) ∂τ k ,l ( R k ) ∂ρ (cid:19) dρ = (cid:90) R n − (cid:32) n − (cid:88) i =0 t i ( k ) φ i ( ρ ) (cid:33) t k ( k ) n − (cid:88) j =0 Ψ (cid:48) , jk φ j ( ρ ) + n − (cid:88) j =0 Ψ jk φ j ( ρ ) (cid:32) n − (cid:88) i =0 t i ( k ) φ i ( ρ ) (cid:33) t (cid:48) k ( k ) dρ , (42)where we have used that τ k ,l ( R k ) = t k ( k ) since thevector t ( k ) represents the solution in the hybrid repre-sentation. In the next step we reorder the sum and the integral and find= − t k ( k ) n − (cid:88) i =0 n − (cid:88) j =0 t i ( k ) Ψ (cid:48) , jk (cid:90) R n φ i ( ρ ) φ j ( ρ ) dρ + t (cid:48) k ( k ) n − (cid:88) j =0 n − (cid:88) i =0 t i ( k ) Ψ jk (cid:90) R n φ j ( ρ ) φ i ( ρ ) dρ ≈ − t k ( k ) n − (cid:88) i =0 t i ( k ) Ψ (cid:48) , ik + t (cid:48) k ( k ) n − (cid:88) i =0 t i ( k ) Ψ ik + O (max i,j To illustrate and benchmark the proposed hybrid JM-ECS method, we consider two model problems that arederived from a two-body and a three-body problem.0When these systems are expressed in spherical coordi-nates around the center of mass and expanded in spheri-cal harmonics one typically ends up with a large systemof coupled radial equations, known as the partial waveequations. These partial waves can either be coupled 1Dor coupled 2D problems. We choose to benchmark ourmethods to model problems formulated as an uncoupledpartial wave problem.The two-body problem reduces to a one-dimensionalradial problem with the form given by (7).Similarly a two-dimensional radial problem is derivedfrom a three-body problem and is, expressed as: (cid:18) T + V ( ρ ) + l ( l + 1)2 ρ + T + V ( ρ ) + l ( l + 1)2 ρ + V ( ρ , ρ ) − k (cid:19) ψ sc ( ρ , ρ )= − √ k (cid:16) [ V ( ρ ) + V ( ρ , ρ )] ϕ l , ( ρ )ˆ j l ( k ρ )+ [ V ( ρ ) + V ( ρ , ρ )] ϕ l , ( ρ )ˆ j l ( k ρ ) (cid:17) ,ψ sc (0 , ρ ) = 0 ∀ ρ > ψ sc ( ρ , 0) = 0 ∀ ρ with outgoing wave boundary conditions for ρ → ∞ or ρ → ∞ , (47)where ρ and ρ are two radial coordinates represent-ing the distances of the first and second particle to thecenter of the coordinate system. The two body poten-tial is V ( ρ , ρ ). The angular momenta l and l arenon-negative integers.The total wave function is the sum of an incomingwave and a scattered wave. If we model an impact ion-ization problem, for example, the incoming wave will bea product of the target in the ground state, ϕ l , and theincoming wave, j l ( k ρ ), of the second particle. Takinginto account symmetrization this is1 √ k (cid:16) ϕ l , ( ρ )ˆ j l ( k ρ ) + ϕ l , ( ρ )ˆ j l ( k ρ ) (cid:17) , where ϕ l , ( ρ ) is an eigenstate of the operator ( T + V + l ( l + 1) / ρ ) with energy E . The incoming waveˆ j l ( k ρ ) has momentum k such that k / E = k / (cid:126)µ applied to theground state of the two particle system ϕ ( ρ , ρ ). A. One dimensional phase shift To validate the proposed hybrid JM-ECS technique wefirst solve the one-dimensional scattering problem withan attractive Gaussian potential: (cid:18) − d dx + l ( l + 1)2 x − V e x /r − k (cid:19) ψ ( x ) = 0 , (48) where the parameters of the potential were chosen as V = 0 . r = 2. We determine the elastic scat-tering phase shift as a function of the momentum of theincoming particle. To estimate the accuracy of the hy-brid results we compare the phase shift to the one ob-tained with the Variable Phase Approach (VPA) [20],which yields the exact result with very high accuracy forthis simple problem. In Fig. 7 we display the results foroscillator length b = 0 . N ; this corresponds to an increasing size of the interac-tion region, and thus an increasing size of the truncatedoscillator basis. It is clear from this figure that the resultsstrongly depend on the choice of the matching point fora required accuracy. In Fig. 8 we display the results fordifferent values of the angular momentum l , all seen tobe of qualitatively comparable accuracy. kδ exactN=10N=50N=100 FIG. 7: s -wave phase shift for a Gaussian potential ( V = 0 . r = 2) as a function of wave number k for the problemof (48). Shown are the hybrid JM-ECS results for differ-ent matching points N , and the VPA (exact) result. Theoscillator length is b = 0 . To quantify the accuracy we display an error plot inFig. 9 where one notices an increasing accuracy with in-creasing size of the oscillator basis in the interaction re-gion. The error is seen to decrease in terms of N in boththe small- and large-energy region. For large energies theerror is mainly caused by the finite-differences approxi-mation to the second derivative. At small energies theerror mainly comes from an inaccurate approximate re-lation (27) which improves as the size of the oscillatorbasis increases.In table I we display the error value for different sizes N of the oscillator basis and for different partial waves.For each partial wave we consider the k -value where theerror is maximalWe observe that the error of the hybrid method hasrather intricate and oscillatory behavior in terms of thesize, but in general decreases with increasing oscillator1 kδ l = 0 l = 1 l = 2 FIG. 8: Phase shifts for a Gaussian potential ( V = 0 . r = 2) of equation (48) and different angular momentaobtained with JM-ECS (squares), compared to VPA results(solid lines). For all results b = 0 . N = 50. −3 k ∆ N=10N=50N=100N=300 FIG. 9: The absolute error in phase shift for different sizes ofthe oscillator basis (b=0.3, l=0) for model problem of fig. 7and 7 basis. Full convergence is however not yet obtained. Inany case it is seen that the approach provides relativelystable and accurate results for different values of the an-gular momenta, energy ranges, and sizes of the oscillatorbasis. B. s -wave benchmark with product two-bodypotentials As described above, the hybrid model introduced inthis paper is easily extended to systems with more de-grees of freedom, in contrast to the original JM method. N | ∆ | ( × − ) l = 0 l = 1 l = 2( k = 0 . 17) ( k = 0 . 62) ( k = 1 . b = 0 . We demonstrate this on a model breakup problem witha short-range potential taken from [21]: (cid:18) − ∂ ∂x − ∂ ∂y − V e − x − V e − y + W e − x − y − E (cid:1) Ψ = 0 (49)with V = 3 and W = 10. These interaction parametersyield one bound state for the attractive “one-particle”potential V with an energy E = − . E inc = 0 . E = 0 . . 471 to share between them.Since the two particles are indistinguishable the cross sec-tion is symmetric around 0 . / 2. In the figure we alsocompare with the results of the finite difference calcula-tions. We see a slight difference for equal energy sharing.To highlight the difference we have scaled the verticalaxis.The error between the finite difference calculations andthe hybrid method is given in table II where we look atthe error for a particular energy sharing. As the numberof oscillator states increases the results converge.In Table II we display the difference between the resultsof the hybrid method and those of a full finite-differencecalculation in terms of the the size of the oscillator basis N . C. More realistic s -wave benchmark with Gaussianpotential All potentials in (49) are exponentially decaying inboth coordinates, and such a model is not very realistic.In real problems the inter-particle potential depending onthe distance between particles should decay significantlyslower. We consider a more realistic three-dimensionalproblem, but similar to the one described above. Thiscould correspond to a scattering problem with two equal,2 C r o ss S e c t i on FIG. 10: SCDS results of s -wave scattering for model problem(49), for a total energy of 0.471 (solid line — finite differences,dashed line — JM-ECS). The parameters of the calculationare N =80 and b =0.3. The vertical axis has been rescaled toemphasize the difference. N | ∆ | 10 0.093720 0.083640 0.067860 0.059280 0.0184TABLE II: Absolute error ∆ of the hybrid method w.r.t. afull finite-difference approach for the SDCS results of problem(49). The total energy is 0.471, and there is equal energysharing between the particles. light, particles and a third heavy one: (cid:18) − 12 ∆ r − 12 ∆ r − V e −| r | − V e −| r | + W e −| r − r | − E (cid:17) Ψ = 0 (50)with r and r the coordinates of the two (light) par-ticles with respect to the third one. Here we considera Gaussian form of the potential for the interaction be-tween the light particles to simplify the subsequent calcu-lations, and to obtain a faster-decaying s -wave potential.The s -wave projection of this equation yields (cid:18) − ∂ ∂x − ∂ ∂y − V e − x − V e − y + W (cid:16) e − ( x − y ) − e − ( x + y ) (cid:17) xy − E Ψ = 0 (51)In Fig. 11 we show the SCDS for this problem, and com-pare with the results from a full finite difference calcu- lation. In these calculations we had to increase the sizeof the finite-difference grid to cover the interaction re-gion, which leads to an increase increasing of the oscil-lator parameter (see 29). We considered b = 1 . b . E C r o ss S e c t i on FIG. 11: SCDS results of s-wave scattering for model problem(51) with a Gaussian two body potential (solid line — finitedifferences, dashed line — JM-ECS). The parameters of thecalculation are N =70 and b =1.3. N | ∆ | 10 0.373020 0.022130 0.077650 0.068570 0.0153TABLE III: Absolute error ∆ of the hybrid method w.r.t. afull finite-difference approach for the SDCS results of problem(51). The total energy is 0.471, and there is equal energysharing between the particles. IV. APPLICATION TO NUCLEAR p -SHELLSCATTERING In this section we go beyond model problems and applyour method in the more realistic context of α -scatteringin p -shell nuclei [22], in particular Be [23].This is a commonly used benchmark problem in nu-clear cluster physics. We will perform an α − α scat-tering calculation including fully microscopic states forthe interaction region, and using a semi-realistic nucleon-nucleon interaction. We compare the results with thoseobtained in the JM approach.In the asymptotic region the system decays into twopoint particles, each corresponding to an α -particle. Atcloser distances, however, the internal structure of theclusters becomes of importance because of the micro-scopic interactions and the Pauli exchanges between nu-cleons.The cluster basis states have the formΨ nL = (cid:98) A { Φ ( α ) Φ ( α ) φ nL ( r ) } , (52)where (cid:98) A is the antisymmetrization operator, and Φ , ( α )is a translation invariant shell-model state built up of s -orbitals for the α -particles. The state φ nL ( r ) is a three-dimensional harmonic oscillator state for the relative mo-tion of the α -clusters φ nL ( r ) = φ n + n ,LM ( r ) (53)= N n + n ,L ρ L e − ρ / L L +1 / n + n (cid:0) ρ (cid:1) Y LM ( (cid:98) r ) ρ = | r | b , N nL = (cid:115) 2Γ ( n + 1)Γ ( n + L + 3 / n to denote the minimal value of the shell num-ber allowed by the Pauli exclusion principle n = (cid:40) ( N min − L ) / L < N min , L ≥ N min , (54)where N min is the number of oscillator quanta in thatshell. N min = (cid:40) A − , normal parity states: π = ( − A A − , abnormal parity states: π = ( − A +1 (55)In this form n = 0 , , ... numerates all the Pauli allowedstates for the cluster relative motion. For all further de-tails, we refer to [22] and [23].A common nucleon-nucleon interaction used for thissytem is a Volkov N1 (V1) potential, which can be writ-ten in the form V ij = (cid:88) k =1 V k (1 + mP rij ) exp {− ( r ij /a k ) } , (56)where m is the Majorana exchange parameter. Thestrength parameters V k determine the repulsive short-range core and the long-range attraction. We take the same parameters as in [23] and [22], but omit theCoulomb interaction between protons to simplify the ef-fective asymptotic α − α interaction.The calculations of matrix elements in the internal re-gion is identical as in the JM approach. The differencelies in the asymptotic region, where the explicit asymp-totic potential, expanded in the oscillator basis, is nowreplaced by the finite differences approach, as discussedin section II C, after proper matching of both regions;the asymptotic α − α interaction reduces to the point-like effective free particle kinetic energy between the twoclusters.In Fig. 12 we show the J π = 0 + , + and 4 + phase shiftsfor the Be system calculated with the hybrid JM-ECSand JM approaches. For the latter we take results asthey have been obtained in [23], where a modified JMapproach was considered, providing faster convergence[19]; we refer to such calculation as MJM for consistency.The overall results are very close to each other. Theonly noticeable discrepancies occur at very low energies.The difference between MJM and JM-ECS phase shiftsin this region is presented in Fig. 13. The oscillations inthis difference correlate with the value of the derivative ofthe wave function in the matching point between the tworepresentations. The source of the error is a “reflection”from this point, and an improvement of the method toreduce this discrepancy is currently under research. E, MeV δ , deg + + + FIG. 12: Scattering phase shifts for the α − α system for differ-ent values of total angular momentum calculated with MJM(solid lines) and hybrid JM-ECS method (closed squares). Allcalculations were made with 80 oscillator states V. DISCUSSION AND CONCLUSION This article reports on initial efforts to introduce ab-sorbing boundary conditions in quantum scattering prob-lems where the wave function is represented by a finitesum of oscillator states. We have introduced the exte-rior complex scaling (ECS) boundary conditions by con-4 E, MeV δ h y b − δ F D FIG. 13: Difference between MJM and JM-ECS scatteringphases for the 2 + state of the α − α system. structing a representation of the wave function that com-bines the oscillator states with a grid based representa-tion. The oscillator representation covers the inner regionof the problem, where the main interaction occurs, whilethe grid covers the near field. The two regions are nu-merically matched at an interface using an asymptoticexpression for the expansion of oscillator states.The far field amplitude, which gives the cross section,is extracted from the numerical wave function using asurface integral. Extra care is required in the calculationof this integral in this representation.We have numerically solved several benchmark prob-lems and compared our results with literature and plainECS calculations.Our results agree with existing methods and confirmthat the proposed method works. However, the accu- racy still needs to be improved. This can be done byincluding mixed potential matrix elements from differ-ent representations, by using a different finite differencegrid, or by using a higher order matching condition at theinterface between the oscillator representation and the fi-nite difference representation. Also a proper treatmentof the Coulomb interaction and other effective interac-tions with long-range behavior must be considered in thefuture development of the method.The building blocks presented in this papers are thestarting point for a fully coupled calculations similar to[13]. In such a calculation the wave function is a sum over l , m and l , m of 2D functions ψ l m ,l m combinedwith Y l m (Ω ) Y l m (Ω ). The linear system is then acoupled system where each diagonal block is a systemlike eq. (47). From such a representation of the wavefunction we can extract not only the single differentialcross section, but also the triple differential cross section,which gives the amplitude for certain breakup directions.In a similar way, the method can be applied to differentmulti-channel calculations of three-cluster systems anddifferent reaction processes with light nuclei. In all suchproblems asymptotic description can be greatly simpli-fied if we use absorbing boundary conditions.We conclude that it is possible to introduce a boundarylayer that absorbs all outgoing waves within a spectralbasis such as the oscillator representation. Acknowledgments This work is supported by FWO Flanders G.0.120.08. [1] B. D. Esry, C. D. Lin, and C. H. Greene, Physical ReviewA, , 394 (1996).[2] A. Cobis, D. V. Fedorov, and A. S. Jensen, PhysicalReview Letters, , 2411 (1997).[3] J. Broeckhove, F. Arickx, P. Hellinckx, V. Vasilevsky,and A. Nesterov, Journal of Physics G: Nuclear and Par-ticle Physics, , 1955 (2007).[4] B. Junker, Adv. Atom. Mol. Phys, , 207 (1982).[5] N. Moiseyev, Physics Reports, , 211 (1998).[6] A. Lane and R. Thomas, Reviews of Modern Physics, ,257 (1958).[7] E. Heller and H. Yamani, Physical Review A, , 1209(1974).[8] C. W. McCurdy, M. Baertschy, and T. N. Rescigno, J.Phys. B, , R137 (2004).[9] J. Berenger, Journal of computational physics, , 185(1994).[10] D. Givoli, Computational Acoustics of Noise Propaga-tion in Fluids-Finite and Boundary Element Methods,145 (2008).[11] B. Simon, Physics Letters A, , 211 (1979).[12] T. N. Rescigno, M. Baertschy, W. A. Isaacs, and C. W.McCurdy, Science, , 2474 (1999).[13] W. Vanroose, F. Martin, T. Rescigno, and C. McCurdy, Science, , 1787 (2005).[14] G. Shortley and R. Weller, Journal of Applied Physics, , 334 (1938).[15] M. Baertschy, T. N. Rescigno, W. A. Isaacs, X. Li, andC. W. McCurdy, Physical Review A, , 22712 (2001).[16] C. W. McCurdy, D. A. Horner, and T. N. Rescigno,Physical Review A, , 22711 (2001).[17] E. Heller and H. Yamani, Physical Review A, , 1201(1974).[18] A. Alhaidari, E. Heller, H. Yamani, and M. E. Abdel-monem, The J-Matrix Method (Springer, 2008).[19] W. Vanroose, J. Broeckhove, and F. Arickx, Phys. Rev.Lett., , 10404 (2001).[20] F. Calogero, Variable phase approach to potential scat-tering (Academic Press, New York, 1967).[21] T. N. Rescigno, C. W. McCurdy, W. A. Isaacs, andM. Baertschy, Phys. Rev. A, , 3740 (1999).[22] A. Sytcheva, J. Broeckhove, F. Arickx, and V. S.Vasilevsky, Journal of Physics G: Nuclear and ParticlePhysics, , 2137 (2006).[23] A. Sytcheva, F. Arickx, J. Broeckhove, and V. S.Vasilevsky, Physical Review C,71