Improved convergence of scattering calculations in the oscillator representation
IImproved convergence of scattering calculations in the oscillatorrepresentation.
Y. Bidasyuk a,b , W. Vanroose a a Departement Wiskunde-Informatica, Universiteit Antwerpen, Antwerpen, Belgium b Bogolyubov Institute for Theoretical Physics, Kyiv, Ukraine
Abstract
The Schr¨odinger equation for two and tree-body problems is solved for scattering states in a hybrid repre-sentation where solutions are expanded in the eigenstates of the harmonic oscillator in the interaction regionand on a finite difference grid in the near– and far–field. The two representations are coupled through ahigh–order asymptotic formula that takes into account the function values and the third derivative in theclassical turning points. For various examples the convergence is analyzed for various physics problems thatuse an expansion in a large number of oscillator states. The results show significant improvement over theJM-ECS method [Bidasyuk et al, Phys. Rev. C , 064603 (2010)]. Keywords: quantum scattering, oscillator representation, Schr¨odinger equation, absorbing boundaryconditions, asymptotic analysis
1. Introduction
The accurate prediction of the breakup of a many-particle system into multiple fragments is one of themost challenging problems in quantum mechanics. Not only the relative motion of the particles needs to bemodeled, but also the internal structure of the target and the products need to be described accurately. Thisleads in many cases, to a high–dimensional Schr¨odinger equation posed on a huge domain. For example, inbreakup or collisions of nuclear clusters the cross section depends on a delicate interplay of the forces thathold the clusters together and the forces between the clusters.The quantum state that describes the internal structure of a bound many particle system is often representedas linear combination of eigenstates of the quantum mechanical oscillator. These states form an L -basiswhich reduces the problem to finding the correct expansion coefficients. Examples are the correlated manyelectron state of a quantum dot [1, 2, 3] and the many nucleon state in nuclear physics [4, 5]. The variousapplications of the oscillator representation are discussed in [6]. Such a representation is very efficient fortightly bound ground state or lowest excited states as any smooth potential well is close to parabolic shapenear it’s bottom. As a result, if the oscillator parameter is optimized to match this local parabolic potential,the lowest states of a system can be efficiently represented with only a few oscillator functions [7].However, the representation is inefficient to describe scattering or break-up processes. Scattering states arenot square integrable and many oscillator states are required to represent the interaction and asymptoticregion. Furthermore, many potential matrix elements need to be calculated and this results in a dense linearsystem that has a complexity of N to arrive at the solution, where N is the number of oscillator states usedin the representation, omitting the cost of calculating the potential matrix elements. Email addresses:
[email protected] (Y. Bidasyuk),
[email protected] (W. Vanroose)
Preprint submitted to Journal of Computational Physics October 9, 2018 a r X i v : . [ phy s i c s . c o m p - ph ] J un he J -matrix method offers a way of calculating cross sections and other scattering observables in theoscillator representation. It was proposed by Heller and Yamani in the seventies [8, 9] and mainly appliedto atomic problems. The method exploits the tridiagonal structure of the kinetic energy operator. The J -matrix has been under constant development since its inception and a review of the recent developmentscan be found in [10]. For problems with Coulomb interactions the Coulomb-Sturmian basis is preferredover the oscillator representation. Recently the Coulomb-Sturmians have been used to describe multiphotonsingle and double-ionization [11] and electron impact ionization [12].At the same time, and mostly independent, the Algebraic Resonating Group Method was developed fornuclear scattering problems [4, 13, 14, 15, 16]. This method exploits the same principles as the J -matrixmethod for description of nuclear cluster systems, where the oscillator representation is efficiently used todescribe the internal structure of clusters. If the same representation is used for intercluster degrees offreedom, then the nucleon symmetrization rules become straightforward.One shortcoming of the methods based on an L basis is that the asymptotic solutions need to be knownexplicitly before a system for the wave function in the inner region and the scattering observables canbe written down. This is a serious limitation since it is hard to find the asymptotic wave function forbreakup reactions with multiple fragments. The Complex Scaling method, which scales the full domain intothe complex plane, can be easily implemented in the oscillator representation by taking a complex valuedoscillator strength. This has been used to calculate energy spectrum or extract the resonances [17], but thecalculation of cross sections and other scattering observables may be quite difficult.In recent decades significant progress has been made in the numerical solution of scattering processes de-scribed by the Helmholtz equation. In contrast to many particle systems, where the potential V is oftennon-local, the wave number k ( x ) in the Helmholtz equation depends only on the local material parameterssuch as the speed of sound in acoustics or the electric permittivity and magnetic permeability for electro-magnetic scattering. For these problems grid based representations such as finite difference, finite element[18, 19], or Discontinuous Galerkin [20, 21] are preferred since they lead to sparse matrices that can be easilysolved by preconditioned Krylov subspace methods [22, 23].Another technique that has found widespread application is the use of absorbing boundary conditions. Theseboundary conditions allow a scattering calculation without prior knowledge about the asymptotic wave form.Exterior Complex Scaling (ECS) and Complex Scaling (CS) are widely used in atomic and molecular physics[24, 25, 26, 27]. Perfectly matched layers (PML) are used for electromagnetic and acoustic scattering [28],which can also be interpreted as a complex stretching transformation [29]. There are many other excellentabsorbing boundary conditions [30, 31, 32].In [33] the JM-ECS method was introduced that combines the J -Matrix method with a grid based ECS. Themethod describes the scattering solution in the interior region with an oscillator representation and in theexterior region with finite differences. The two representations are matched through a low order asymptoticformula with an error that scales as O ( N − / ), where N is the size of the oscillator basis describing theinner region. Once the grid and oscillator representation are matched, it is easy to introduce an absorbingboundary layer since the grid representation can be easily extended with an ECS absorbing layer or anyother absorbing boundary condition.The resulting method was illustrated for one– and two–dimensional model problems representative for realscattering problems with local interactions. Furthermore, the representation was used for nuclear p -shellscattering. However, the accuracy of the calculations was unsatisfactory due the low–order matching condi-tion. While the grid representation on the exterior has an accuracy of O ( N − ), the matching is only accurateto order O ( N − / ).The main contribution of the current paper is to increase the accuracy of the asymptotic formula that allowsa better matching of the grid and oscillator representations. The better asymptotic formula takes function2alues and its third derivative into account. It can bring the matching error down to the level of the accuracyof the grid representation.A higher–order asymptotic approximation of the oscillator representation was already discussed by S. Igashovin [34]. However, the formula was not used to increase the accuracy of scattering calculations.This paper is outlined as followed. In Section 2 we shortly describe the process of scattering calculationsin the oscillator representation. In Section 3 we derive a higher–order asymptotic formula that takes intoaccount the behavior of the function in the classical turning points in coordinate and Fourier space. InSection 4 and 5 we use this asymptotic formula to solve scattering problems.
2. Review of scattering calculations in the oscillator representation
In this section we discuss the most important properties of the oscillator representation and how they can beused to perform scattering calculations with the J -matrix method. We also recall the working of the hybrid J -matrix and ECS method proposed in [33]. The aim is to solve the Schr¨odinger equation in atomic units ( m = 1, (cid:126) = 1) that describes a scatteringprocess of two particles for an energy E ∈ R . This equation written in relative coordinates is (cid:20) −
12 ∆ + V ( r ) − E (cid:21) Ψ( r ) = F ( r ) , ∀ r ∈ R (1)where F ( r ) is the function describing the initial state or the source term and V ( r ) is the potential. Thecoordinate r can be written in spherical coordinates ( ρ, θ, φ ). In case of spherically symmetric potential( V ( r ) = V ( ρ )), Eq. (1) can be reduced to a one-dimensional radial equation using partial wave decompositionΨ( r ) = Ψ( ρ, θ, φ ) = (cid:88) l,m ψ l ( ρ ) ρ Y l,m ( θ, φ ) , F ( r ) = F ( ρ, θ, φ ) = (cid:88) l,m f l ( ρ ) ρ Y l,m ( θ, φ ) , where Y l,m ( θ, φ ) is a spherical harmonic, l = 0 , , , . . . is the orbital angular momentum of the relativemotion, m is the projection of this angular momentum. The resulting reduced radial equation becomes (cid:20) − d dρ + l ( l + 1)2 ρ + V ( ρ ) − E (cid:21) ψ l ( ρ ) = f l ( ρ ) . (2)For the problem we are interested in there is a range a > ρ > a both V ( ρ ) and χ ( ρ ) arezero.We solve the Eq. (2) for E > ψ l scattering observables such as the crosssections or phase shift. In order to solve the Eq. (2) we represent the solution as ψ l ( ρ ) = ∞ (cid:88) n =0 c n,l ϕ n,l ( ρ ) , (3)where ϕ n,l ( ρ ) are orthogonal L functions, in particular we will use reduced oscillator functions, whoseproperties will be explained in the next section. After projection of (3) on ϕ n,l , Eq. (2) results in an infinitelinear system ∞ (cid:88) n =0 (cid:16) T ( l ) kn + V ( l ) kn − E (cid:17) c n,l = b k,l , (4)3here T ( l ) kn denotes the elements of the kinetic energy matrix and V kn the elements of the potential energymatrix. They are the integrals T ( l ) kn = (cid:90) ∞ ϕ k,l ( ρ ) (cid:18) − d dρ + l ( l + 1)2 ρ (cid:19) ϕ n,l ( ρ ) dρ,V ( l ) kn = (cid:90) ∞ ϕ k,l ( ρ ) V ( ρ ) ϕ n,l ( ρ ) dρ and b k,l = (cid:82) ∞ ϕ k,l ( ρ ) χ l ( ρ ) dρ . As the considered problem is effectively one-dimensional we will use x insteadof ρ to denote radial relative coordinate in all following sections devoted to two-body problem. Before we explain the strategy to solve Eq. (4), we repeat the main properties of the oscillator representationand the function ϕ n,l that will be used in the expansion.The reduced radial equation analogous to (2) for the quantum harmonic oscillator with an oscillator strength ω is (cid:20) − d dx + 12 l ( l + 1) x + 12 ω x (cid:21) ϕ n,l ( x ) = E n,l ϕ n,l ( x ) (5)with x ∈ [0 , ∞ [ a radial coordinate. The boundary conditions are ϕ n,l (0) = 0 and lim x →∞ ϕ n,l ( x ) = 0. Theeigenvalues are E n,l = (cid:18) n + l + 32 (cid:19) ω, (6)and eigenstates ϕ n,l ( x ) = ( − n N n,l b − / (cid:16) xb (cid:17) l +1 exp (cid:18) − x b (cid:19) L l +1 / n (cid:18) x b (cid:19) , (7)where L l +1 / n are Laguerre polynomials. The normalization is N n,l = (cid:112) n ! / Γ( n + l + 3 / n ∈{ , , . . . , } and oscillator length is defined as b = (cid:112) /ω .The classical turning point associated with each state is R n,l = b √ n + 2 l + 3, defined as the point wherethe potential energy equals the total energy of the system.The functions (7) form a complete orthonormal basis and any wave function ψ l ( x ) that behaves as x l in x = 0 can be represented using the infinite sum: ψ l ( x ) = ∞ (cid:88) n =0 c n,l ϕ n,l ( x ) , where c n,l = (cid:90) ∞ ϕ n,l ( x ) ψ l ( x ) dx. (8)In the following section we will use that the radial oscillator equation, (5), can be rewritten in terms of b ,the oscillator length, and R n,l , the classical turning point as (cid:34) − d dx + l ( l + 1) x + x b − R n,l b (cid:35) ϕ n,l ( x ) = 0 . (9)An important property that forms the basis of the results in this paper is that the n -th oscillator state has n oscillations between the origin and the classical turning point R n,l = √ n + 2 l + 3. Beyond this turningpoint the function is exponentially decaying without additional oscillations. This means that as n increases,the frequency of the oscillation between the origin and R n,l grows proportional to √ n . This property willbe used to derive the asymptotic formula in Section 3.The Bessel transform ˜ F l ( k ) of a function F l ( x ) is defined as˜ F l ( k ) := (cid:114) π (cid:90) ∞ F l ( x )ˆ j l ( kx ) dx. (10)4 x x φ n , l ( x ) R n,l x x φ n , l ( x ) R n,l Figure 1: The reduced oscillator state ϕ n,l ( x ) with a classical turning point R n,l for l = 0 and two values of n : n = 10 (left) and n = 30 (right). The function oscillates n times on the interval [0 , R n,l ]. Since R n,l only grows with O ( √ n ) the state becomesrapidly oscillating as n grows. where ˆ j l ( kx ) is a Riccati–Bessel function, the regular solution of the radial Schr¨odinger equation withoutpotential, (2). It is connected with the ordinary Bessel function in a following way ˆ j l ( x ) = (cid:112) πx/ J l +1 / ( x ).The Bessel transform of an oscillator state is again an oscillator state˜ ϕ n,l ( k ) = ( − n b ϕ n,l ( kb ) , (11)with a classical turning point K n,l = R n, /b = √ n + 2 l + 3 /b . This relation that can easily be derived using e.g. formula (7.421.4) from [35]. The oscillator states therefore form an orthonormal basis of L ([0 , ∞ [) thatdiagonalizes the Bessel transform. Similar properties hold for the Cartesian oscillator state based on Hermitepolynomials, see, for example, [3].An other important property that will be used in the remainder of the text is the following: Proposition 2.1.
Let ψ l a function that behaves as x l in x = 0 . The projection on the oscillator state ϕ n,l can be calculated with either ψ l or its Bessel transform ˜ ψ l as: c n,l = (cid:90) ∞ ϕ n,l ( x ) ψ l ( x ) dx = ( − n b (cid:90) ∞ ϕ n,l ( x ) ˜ ψ l ( x ) dx. (12) Proof.
Since Parseval’s theorem holds we can calculate c n either with ϕ n,l or its Bessel transform c n,l = (cid:90) ∞ ϕ n,l ( x ) ψ l ( x ) dx = (cid:90) ∞ ˜ ϕ n,l ( k ) ˜ ψ l ( k ) dk. (13)But since ˜ ϕ n,l ( k ) = ( − n bϕ n,l ( kb ), the oscillator state in the latter integral is the same as in the firstintegral, after substitution of variables and up to a constant.This result will be used in the next section to derive the asymptotic formula for c n,l with large n . Thissymmetry between ψ l ( x ) and its Bessel transform ˜ ψ l ( k ) should also hold in the asymptotic formula.The kinetic energy operator T li,j is a tri-diagonal matrix. For the oscillator basis the non-zero elements are T li,j = (cid:90) ∞ ϕ i,l ( x ) (cid:18) − d dx + l ( l + 1)2 x (cid:19) ϕ j,l ( x ) dx = (cid:0) i + l + (cid:1) ω for j = i, − (cid:113) i (cid:0) i + l + (cid:1) ω for j = i − , − (cid:113) ( i + 1) (cid:0) i + l + (cid:1) ω for j = i + 1 . (14)5or the remainder of the text we will consider only the case of zero angular momentum and drop l from thenotation. All results, however, are valid for arbitrary l . In [36] an asymptotic formula was derived for the projection of a smooth radial function ψ on the state ϕ n .The derivation uses stationary phase arguments to exploit the increase of oscillations as n → ∞ . We repeathere the derivation of the results.There are two main contributions to the value of the integral c n = (cid:90) ∞ ϕ n ( x ) ψ ( x ) dx ≈ I + I R n . (15)There is a first contribution, denoted I , of the left integration boundary near zero. There is a secondcontribution, denoted I R n , from the point of stationary phase, which is the classical turning point wherethe oscillations stop (see Figure 1). Here we provide only general steps of calculating this integral in theasymptotic region. More details can be found in [36].The contribution from the classical turning point can be calculated by approximating the oscillator state ϕ n near the classical turning point by an Airy function as ϕ n ( x ) ≈ b (cid:18) b R n (cid:19) / Ai (cid:34)(cid:18) R n b (cid:19) / ( x − R n ) (cid:35) if | x − R n | (cid:28) . (16)The integral representation of this Airy function and the stationary phase approximation leads to the con-tribution of the turning point to the integral I R n ≈ b (cid:114) R n ψ ( R n ) . (17)The contribution of the left integration boundary, I , can be derived by approximating the oscillator statenear the origin by a Riccati–Bessel function ϕ n ( x ) ≈ ( − n √ b (cid:114) πK n ˆ j ( K n x ) , (18)where K n is the classical turning point of the oscillator state in the momentum space. Then the contributionfrom the origin becomes I ≈ ( − n b (cid:114) K n ˜ ψ ( K n ) . (19)And the resulting asymptotic approximation of the oscillator coefficients becomes c n ≈ ( − n b (cid:114) K n ˜ ψ ( K n ) + b (cid:114) R n ψ ( R n ) if n (cid:29) . (20)This relation has a contribution from the turning points in the coordinate space, R n and the Fourier space, K n . Note that it satisfies the symmetry observed in 2.1 above.Note that Eq. (20) as well as the similar equation in [36] does not provide the order of the approximation.This is one of the shortcomings that will be addressed in the current paper. We now discuss how a finite linear system can be obtained that solves for the wave function in the interactionregion and the phase shift, describing the solution in the asymptotic region. The presentation here is based6n the asymptotic formula and differs from how the method was derived historically. For more detail aboutthe formulation of the original J -matrix method we refer to [9] and [14].Let ψ ( x ) be a smooth two-body radial scattering state with a bounded energy, depending on one spatialcoordinate (can be always reduced using center-of-mass relative coordinates and spherical coordinates). Sinceit is a scattering state, the function does not go to zero as x → ∞ . However, its Bessel transform ˜ ψ ( k ) goesto zero as k → ∞ . This means that, as n grows, the only contribution to the expansion coefficient c n comesfrom the classical turning point in coordinate space. So, for a scattering state ψ with total energy E = k / j l ( kr ) + tan( δ l )ˆ n l ( kr ) where δ l is the phase shift [37] holds that c n ≈ b (cid:114) R n ψ ( R n ) = b (cid:114) R n (cid:16) ˆ j l ( kR n ) + tan( δ l )ˆ n l ( kR n ) (cid:17) , (21)where ˆ j l and ˆ n l are the regular and irregular solutions of the free-particle equation. The c n becomes therepresentation of ψ on the grid of classical turning points R n .The aim of a one-dimensional radial scattering calculation is to find a numerical approximation to tan( δ l )for a given potential V . This can be achieved by writing the solution in the oscillator representation as c n,l = c n,l + j n,l + tan( δ l ) n n,l n < Nj n,l + tan( δ l ) n n,l n ≥ N, (22)where j n,l is the regular solution of the homogeneous three-term recurrence relation T n,n − j n − ,l + ( T n,n − E ) j n,l + T n,n +1 j n +1 ,l = 0 , (23)where as n (cid:29) j n,l = b (cid:112) /R n ˆ j l ( kR n ). The n n,l is the irregular solution of the recurrencerelation that goes as n n,l = b (cid:112) /R n ˆ n l ( kR n ) when n becomes large.With the form (22) we reduce the infinite linear system to a finite dimensional problem with a set of N + 1unknowns { c i,l , tan δ l } . Once solved we simultaneously obtain the wave function of the system, { c i,l } , andthe scattering information, tan( δ l ).A detailed description of the construction of this linear system, known as the J -matrix, can be found in [14]. J -matrix and ECS method In [33] the asymptotic formula (20) was used to introduce a hybrid oscillator and grid representation for theSchr¨odinger equation that is useful for scattering calculations where an asymptotic form such as (21) is notexplicitly known. Such a situation appears in breakup reactions of three or more particles.Let ψ ( x ) be again a smooth one-dimensional radial scattering state with a bounded energy such that c n ≈ b (cid:114) R n ψ ( R n ) . (24)holds. The c n becomes the representation of ψ on the grid of classical turning points R n . The grid distancebetween these points becomes smaller when n increases since h n = R n − R n − ≈ b R n . (25)However, the asymptotic ψ does not necessarily need to be represented on the grid of turning points R n .Another option is, for example, a regular grid of equally spaced points.The hybrid JM-ECS method represents the one-dimensional radial wave function as a vector Ψ in C N + K ,where Ψ = ( c , c , . . . , c N − , ψ ( R N ) , ψ ( R N + h ) , . . . , ψ ( R N +( K − h )) . (26)7he first N elements represent the wave function in the oscillator representation. While the remaining K elements represent ψ on an equidistant grid that starts at R N , the N -th classical turning point, and runsup to R N + ( K − h with a grid distance h equal to the difference of the last two turning points of theoscillator representation h = R N − R N − . It is assumed that the matching point that connects the oscillatorto finite–difference representation corresponds to a large index N such that the asymptotic formula, (24),applies.Again, the kinetic energy operator in this hybrid representation is tridiagonal since in both finite–differenceand oscillator representations it is tridiagonal. One should only be careful in matching both representations.To obtain the kinetic energy in the last point of the oscillator representation, the tridiagonal kinetic energyformula (14) is used. It involves a recurrence relation connecting the three terms c N − , c N − and c N .The latter, the coefficient c N , is unknown. Only ψ ( R N ) is available on the grid. Using the asymptoticrelation (24), however, it is possible to calculate the required matrix element as follows:( T c ) N − = T N − ,N − c N − + T N − ,N − c N − + T N − ,N b (cid:112) /R N ψ l ( R N ) . To calculate the kinetic energy in the first point of the finite difference grid, the second derivative of the wavefunction has to be known. To approximate the latter with a finite difference formula, one needs the wavefunction in the grid points R N − , R N and R N + h . Again it is possible to apply (24) to obtain ψ ( R N − ) interms 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 the matching point is sketched in Figure 2, together withthe terms involved to determine the correct matching.Oscillator Finite Differences
T c 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) Figure 2: Illustration on the calculation of the kinetic energy matrix elements that are calculated in the last point R N − of theoscillator representation and in the first point R N of the finite difference representation. To calculate T applied to a solutionvector we need to translate the oscillator representation to the grid. Vice versa for the application of the finite difference stencilapproximation to the second derivative. However, in the next section we will show that the asymptotic matching condition that is used in [33] to coupleboth representations, is only accurate to O ( R − N ), while the outer grid is accurate to order O ( h ) = O ( R − N ).Therefore the largest error is made in the matching condition that couples both representations.Note that introducing an absorbing layer is easy once the oscillator representation is coupled to a gridrepresentation. For example ECS is implemented by extending the grid with a complex scaled part [24].
3. Higher–order asymptotic formula
In this section we derive a higher–order asymptotic formula. It not only takes into account the functionvalues in the turning point R n but also the third derivative in this point. A similar asymptotic formula wasderived by S. Igashov in [34], however, it did not include the contributions from the Fourier space.8 roposition 3.1. Let ψ ( x ) be a regular scattering state that behaves as x l in x = 0 , that is infinitelydifferentiable and ϕ n ( x ) is a reduced oscillator state, solution of (9) . Then the projection of the scatteringstate on the oscillator state can be approximated as c n = (cid:90) ∞ ϕ n ( x ) ψ ( x ) dx = b (cid:114) R n (cid:18) ψ ( R n ) + b R n ψ (cid:48)(cid:48)(cid:48) ( R n ) (cid:19) + O ( R − / n ) . (27) This expresses the expansion coefficient in terms of the wave function and its derivatives in the classicalturning point R n .Proof. The oscillator Eq. (9) can be rewritten as (cid:20) − d dx + 2 R n ( x − R n ) b + ( x − R n ) b (cid:21) ϕ n ( x ) = 0 . (28)For R n (cid:29) | x − R n | (which means either n (cid:29) | x − R n | ≈ ϕ n ( x ) ≈ ϕ (0) n ( x ) = 2 b (cid:18) b R n (cid:19) / Ai (cid:34)(cid:18) R n b (cid:19) / ( x − R n ) (cid:35) , (29)where the normalization is chosen such that it coincides with the oscillator state near the classical turningpoint R n .This can be further improved by writing ϕ n ( x ) = (1 + u n ( x )) ϕ (0) n and inserting this in the oscillator equation.The equation for u n ( x ) becomes − u (cid:48)(cid:48) n ϕ (0) n − u (cid:48) n ϕ (0) n (cid:48) − u n ϕ (0) n (cid:48)(cid:48) + (cid:18) R n ( x − R n ) + 12 ( x − R n ) (cid:19) u n ϕ (0) n = −
12 ( x − R n ) ϕ (0) n . (30)Using that ϕ (0) n (cid:48) = (cid:0) R n /b (cid:1) / Ai (cid:48) (cid:104)(cid:0) R n /b (cid:1) / ( x − R n ) (cid:105) and ϕ (0) n (cid:48)(cid:48) = (cid:0) R n /b (cid:1) / Ai (cid:48)(cid:48) (cid:104)(cid:0) R n /b (cid:1) / ( x − R n ) (cid:105) the equation (30) becomes − u (cid:48)(cid:48) n ϕ (0) n − u (cid:48) (cid:18) R n b (cid:19) / Ai (cid:48) (cid:34)(cid:18) R n b (cid:19) / ( x − R n ) (cid:35) − u n (cid:18) R n b (cid:19) / Ai (cid:48)(cid:48) (cid:34)(cid:18) R n b (cid:19) / ( x − R n ) (cid:35) + (cid:18) R n ( x − R n ) + 12 ( x − R n ) (cid:19) u n ϕ (0) n = −
12 ( x − R n ) ϕ (0) n . (31)In the limit n → ∞ , the term with R n ( x − R n ) dominates compared to (1 / x − R n ) and the otherterms since both the first and second derivative of Ai( x ) around zero are bounded. We find that u n ( x ) ≈− ( x − R n ) / (2 R n ) is a solution of the remaining equation.At this time we will not consider the contributions to the integral for the boundary at zero. This contributionwill be equal to the contribution from the Fourier space as in Eq. (20) and for scattering states and large n this contribution is negligible. So we can lower the integration boundary and write c n ≈ √ a n b (cid:90) ∞−∞ ψ ( x + R n )Ai( x/a n ) (cid:18) − x R n (cid:19) dx, (32)where a n := ( b / R n ) / and we have substituted the integration variable x → x + R n .9he value of this integral will be determined by the behavior around R n , the point of stationary phase. Sincethe function ψ is infinitely differentiable we can Taylor expand c n ≈ √ a n b (cid:90) ∞−∞ ∞ (cid:88) m =0 x m m ! ψ ( m ) ( R n ) (cid:18) − x R n (cid:19) Ai( x/a n ) dx (33)All integrals in these series can be calculated explicitly using specific properties of Airy function (see page52 of [38]). (cid:90) ∞−∞ Ai( x ) x k dx = (3 k )!3 k k ! and (cid:90) ∞−∞ Ai( x ) x k +1 dx = (cid:90) ∞−∞ Ai( x ) x k +2 dx = 0 . (34)We finally get c n ≈ a / n b ∞ (cid:88) k =0 a kn k k ! (cid:18) ψ (3 k ) ( R n ) − a n R n ψ (3 k +2) ( R n ) (cid:19) = b (cid:114) R n ∞ (cid:88) k =0 k ! (cid:18) b R n (cid:19) k (cid:18) ψ (3 k ) ( R n ) − b R n ψ (3 k +2) ( R n ) (cid:19) . (35)In lowest order (in terms of 1 /R n ) this gives us exactly the initial relation (24). Including the next ordercorrection we get a relation (27). For testing purposes we also include terms of the next order 1 /R n c n = b (cid:114) R n (cid:18) ψ ( R n ) + b R n ψ (cid:48)(cid:48)(cid:48) ( R n ) + 1 R n (cid:18) b ψ (6) ( R n ) − b ψ (cid:48)(cid:48) ( R n ) (cid:19) + O ( 1 R n ) (cid:19) (36)Note that a similar result is readily obtained for Cartesian harmonic oscillator states that are based onthe Hermite polynomials. There, however, there is a contribution from both turning points, one from R n = √ n + 1 and R n = −√ n + 1. Corollary 3.2.
Let ψ ( x ) be a function that behaves as x l in x = 0 that is infinitely differentiable, then theasymptotic expansion coefficient is c n = (cid:90) ∞ ϕ n ( x ) ψ ( x ) dx = b (cid:114) R n (cid:18) ψ ( R n ) + b R n ψ (cid:48)(cid:48)(cid:48) ( R n ) (cid:19) + ( − n b (cid:114) K n (cid:18) ˜ ψ ( K n ) + 16 b K n ˜ ψ (cid:48)(cid:48)(cid:48) ( K n ) (cid:19) + O ( n − / ) . (37) Proof.
The asymptotic formula should have the same result when we interchange ψ ( x ) with its Bessel trans-form ˜ ψ ( k ) as a result of Proposition 2.1. Indeed, the integral c n = ( − n b (cid:90) ∞ ϕ n ( k ) ˜ ψ ( k ) dk will have a contribution from integral boundary near 0 and a contribution from the turning point K n inFourier space. The latter is, with the help of the previous results, ( − n b (cid:112) /K n ( ˜ ψ ( K n ) + ˜ ψ (cid:48)(cid:48)(cid:48) ( K n ) / b K n ).While the contribution from boundary become the contribution from the turning point R n in coordinatespace. Example 3.3.
We illustrate the convergence of the asymptotic formula (37) with an application to thefunction ψ ( x ) = exp( − ax ) sin( x ) . We give results for two choices of a . As the scale of the representationis defined by the oscillator length b , the result will change with its choice. If the product ab is large, the −8 −6 −4 −2 n E rr o r n −5/4 −10 −8 −6 −4 −2 n E rr o r n −5/4 Figure 3: Convergence of the general asymptotic formula (37) with only coordinate terms (crosses), only Fourier terms (circles),both coordinate and Fourier terms (solid line) for test function ψ ( x ) = exp( − ax ) sin( x ) with different values of a : left figure — a = 0 .
2, right figure — a = 1 . b = 1 in both cases) Fourier components will dominate in the expansion coefficient. When ab is small the function ψ resemblesa scattering state, on a scale defined by b , and the coordinate approximation will dominate. However, asFigure 3 illustrates the combined formula always gives the correct result. We see an overall convergence with n − / . This is smaller than O ( h n ) which goes as O (1 /n ) . In general, it is possible to construct a function,for which both coordinate and Fourier terms fail to represent the exact oscillator coefficient, but the combinedformula remains valid even in this case (see Figure 4). Example 3.4.
We have verified the asymptotic relations derived in the previous section with a numericalexperiment using the wave function ψ ( x ) = sin( kx ) . We calculate oscillator coefficients directly and comparethem against the values obtained with asymptotic formulas of different order derived previously (see Fig. 5).3.1. Inverse relation The asymptotic expansion coefficient is an approximation with the help of the functions values of ψ evaluatedat certain grid points. It is also useful to derive in inverse relation that allows us to construct the functionvalue in certain grid points given the expansion coefficients.To approximate c n with Eq. (27) the function value and its third derivative need to be calculated at theturning point R n . When only the function values of ψ are available at the turning points R n , the thirdderivative can be approximated by finite differences.The coefficient is then calculated as c n = b (cid:114) R n (cid:32) ψ ( R n ) + b R n (cid:88) k D (3) k ψ ( R n + k ) (cid:33) + O ( R − / n ) + R − / n (cid:15) n , (38)11 −8 −6 −4 −2 n E rr o r n −5/4 Figure 4: Convergence of the general asymptotic formula (37) for the test function ψ ( x ) = exp( − ax ) cos( x ) with a = 0 .
2, forwhich neither coordinate (crosses) nor Fourier terms (circles) give accurate result, while the combined formula (solid line) does. −8 −6 −4 −2 n E rr o r n −5/4 n −3/4 n −7/4 −5 −4 −3 −2 n E rr o r n −1 n −1/2 Figure 5: Convergence of the asymptotic formula (36) (left) and the inverse asymptotic formula (41) (right) with differentnumber of terms included (crosses — one term, solid line — two terms, circles — three terms). The test function has the form ψ ( x ) = sin( kx ) with k = 1 . b = 0 . D (3) is the matrix with coefficients that approximates the third derivative and k indicates the stencilpoints and (cid:15) n is the error term of this approximation. In case of 5-point finite difference approximation k ∈ {− , − , , , } . Note that the grid of turning points R n is an irregular grid and the coefficients willdepend on the local distances between the neighboring grid points. The error of the approximation (cid:15) n willalso depend on the local grid distances as O ( h p ) ∼ O ( R − pn ), where p is the order of the approximation.The resulting transformation can be represented by a banded sparse matrix Uc n = (cid:88) k U nk ψ ( R k ) + O ( R − / n ) . (39)Note that the matrix elements of U will only depend on the values of R n . Indeed, The differential operator D (3) only depends on the grid distances and so do the coefficients of the asymptotic expression.The relation (38) translates ψ on the grid of turning points to corresponding c n . It is also useful to derive theinverse relation that gets ψ ( R n ) from known values of c n . We can not use a direct inversion of Eq. (27) sinceit involves the third derivative, a dense operator in the oscillator representation. However, an approximateinverse relation can be obtained by rearranging terms in (38) as ψ ( R n ) = 1 b (cid:114) R n c n − b R n (cid:88) k D (3) k ψ ( R n + k ) + O ( R − n ) + R − n (cid:15) n (40)and replacing values of the wave function in the right–hand side by oscillator coefficients using only the firstterm of (24), i.e. ψ ( R n ) = (1 /b ) (cid:112) R n / c n + O ( R − n ). This only introduces an error of the order of O ( R − n )but combined with the 1 /R n this gives the approximate inverse relation ψ ( R n ) = 1 b (cid:114) R n c n − b R n (cid:88) k D (3) k (cid:114) R n + k c n + k + O ( R − n ) (41)that is accurate to O ( R − n ). An example of the resulting numerical accuracy of this inverse asymptoticrelation is shown on the right panel of figure 5. Also this transformation can be presented as a sparse bandedmatrix multiplication ψ ( R n ) = (cid:88) n W nk c k + O ( R − n ) . (42)Again the matrix elements of W only depend on the values of the turning points R n .Combining the two relations leads to an approximate partition of unity: U W = I + O ( R − n ) . (43)It is important to note that since these transformation matrices only depend on values of R n they can, inprinciple, both be defined for an arbitrary grid. With the help of these two transformation matrices U and W we can now build an approximate oscillatorrepresentation of an operator Q using its finite difference representation. Let Q ( fd ) be the finite differencerepresentation of Q on the grid R n . Then the application of the Q on ψ can be written as( Qψ ) ( osc ) n = (cid:88) m Q ( osc ) nm c m = (cid:88) m ˇ Q ( osc ) nm c m + O ( n − ) , where ˇ Q ( osc ) nm := [ U Q ( fd ) W ] nm . We illustrate the accuracy of this relation on a second derivative operator D ( −6 −4 −2 n E rr o r n −1 Figure 6: Convergence of the approximate second derivative operator ˇ D = UD ( fd ) W (crosses) and the approximate identityoperator ˇ I = UW (circles) with a test function of the form ψ ( x ) = sin( kx ) ( k and b are the same as in Figure 5) ˇ D (2)( osc ) = U D (2)( fd ) W, (44)where D (2)( fd ) is a finite difference matrix of the second derivative. To analyze the accuracy of this approx-imation Figure 6 shows the result of the error operator defined as ( D ( osc ) − ˇ D ( osc ) ) acting on the vector ofoscillator coefficients of the test function ψ ( x ) = sin( kx ). For comparison we show the error of the approx-imate identity operator, that can be defined as ( I − ˇ I ) = ( I − U W ). We see that both considered erroroperators decay as O ( n − ) in high- n region. This means that our approximate operators are asymptoticallyequal to the exact ones and we can expect the same order of convergence in a solution of the scatteringproblem.It is important to note that all matrices in (44) can be built for any arbitrary spatial grid, though it willnot be an approximate oscillator representation. But we can modify this grid only in the asymptotic region( e.g. to implement the absorbing boundary with ECS). Then the coefficients c n have no physical meaningas oscillator coefficients but the coordinate wave function can sill be reconstructed with (42).
4. High-order hybrid representation for scattering problems
The aim is now to build a hybrid representation for scattering calculations where the oscillator basis is used inthe internal region and finite differences in the outer region. The matching of the two should use the higher–order asymptotic formula (27). However, the strategy displayed on Fig. 2 cannot be easily applied since thethird derivative of the wave function needs to be estimated from several neighboring points symmetricallydistributed on both sides of the matching point. But in the matching point we only have the requiredfunction values on one side of the matching point. This arguments holds both for the oscillator and for thefinite difference representation. We therefore present a matching strategy that differs from the proposal of[33].Consider a arbitrary grid ¯ R n in [0 , ∞ ] with n = 1 , , . . . that differs from the grid of turning points R n . Onthis grid we can construct the differential operator D (3) and construct the operators ¯ U and ¯ W by replacing14very appearance of R n in equations (38) and (41) by ¯ R n . Similarly, for a function ψ sampled in the gridpoints ¯ R n we can calculate ¯ c n = (cid:88) k ¯ U nk ψ ( ¯ R k ) . (45)These coefficients ¯ c n are not the expansion coefficients of the function ψ in the oscillator function basis. Onlywhen ¯ R n equals R n , the oscillator turning points, then the ¯ c n are approximations to the oscillator expansioncoefficients c n .To build the hybrid representation, we choose the grid ¯ R n such that first N + 1 points correspond to theturning points R n . The remaining points of ¯ R n , for n > N + 1, are chosen to correspond to points on aequidistant finite difference grid. To solve a scattering problem the equation (2) needs to be discretized withfinite differences on the grid ¯ R n and then transformed with the help of ¯ U and ¯ W into (cid:88) j (cid:0) ¯ U H fd ¯ W − E ¯ U ¯ W (cid:1) ij ¯ c j = ( ¯ U f ) i (46)to arrive at an equation for ¯ c i .The first N coefficients of ¯ c n correspond now to approximations to oscillator expansion coefficients c n .However, because n is low, they are only a poor approximation. The idea is now to replace the first N coefficients with the exact coefficients. At the same time we replace the first N rows of the matrix with theexact operator in the oscillator representation. The linear system then becomes H ( osc )00 − E . . . H ( osc )0 N H ( osc )0 N +1 . . .H ( osc )10 . . . H ( osc )1 N H ( osc )1 N +1 . . . ... ... ... H ( osc ) N . . . H ( osc ) NN − E H ( osc ) NN +1 . . . [ ¯ UH ( fd ) ¯ W ] N +1 , . . . [ ¯ UH ( fd ) ¯ W ] N +1 ,N [ ¯ U ( H ( fd ) − E ) ¯ W ] N +1 ,N +1 . . . ... ... ... c c ... c N ¯ c N +1 ... = f ( osc )0 f ( osc )1 ... f ( osc ) N (cid:80) j ¯ U N +1 j f ( x j )... , (47) where H ( osc ) is the representation of the Hamiltonian in the oscillator representation and H ( fd ) in the finitedifference representation.We emphasize the difference with [33]. Here we do not match two regions by using the asymptotic formula inone point. Now the representation in the asymptotic region is a fairly good approximation of the oscillatorrepresentation. Therefore the structure of the discretized wave function is simpler than (26). Now we have¯ Ψ = ( c , . . . c N , ¯ c N +1 , . . . ¯ c K ) , where in the initial version of hybrid method, the latter elements of the vector are the function values of ψ inthe grid points. Now the internal region is covered by an exact oscillator representation, and the asymptoticregion is covered by approximate representation which is based on finite differences and includes the ECStransformation. We illustrate the method for a one-dimensional radial example. After the solution of (47), we first reconstructthe coordinate wave function using (42). Outside the range of the potential V the solution can be written asa linear combination ψ sc = A ˆ h + l ( kx )+ B ˆ h − l ( kx ), where ˆ h ± l are the in- and outgoing Riccati-Hankel functions.The coefficient A is then extracted with A = W (cid:16) ψ sc ( x ) , ˆ h − l ( kx ) (cid:17) /W (cid:16) ˆ h + l ( kx ) , ˆ h − l ( kx ) (cid:17) , (48)where x is outside the range of the potential but still on the real part of the ECS domain. The Wronskianis calculated as W ( u, v ) = u (cid:48) v − v (cid:48) u , where the derivatives can be implemented with finite differences. Fromthe Wronskians for A and B we can extract the phase shift of the solution.15 E δ , deg E E rr o r , deg Figure 7: Scattering phase shift of the model problem of Section 4 with a Gaussian potential (left) and the absolute error ofthe scattering phase shift calculated with JM-ECS (dashed line) and our new approach (solid line) depending on the energy ofthe system (right). Both calculations were made with 100 functions in the oscillator basis with the oscillator length b = 0 . −6 −5 −4 −3 N osc E rr o r N −1/2 N −1 −3 −2 −1 N osc E rr o r N −1/2 N −1 Figure 8: Accuracy of the scattering phase shift calculated with JM-ECS (crosses) and the new method (circles) depending onthe size of the oscillator basis. Calculations were made for two values of energy: E = 0 .
2, where δ = 43 . ◦ (left) and E = 3,where δ = 19 . ◦ (right) It is important to note that the accuracy of the method does not directly depend on the size of the asymptoticregion as long as this region is located outside the range of the potential ( V = 0 in the considered part ofˇ H ( osc ) ). In this region we use a grid of 150 equidistant points and then the ECS layer that spans 10dimensionless length units and applies a complex rotation of 45 ◦ to the coordinate axis.We first consider a model problem with a attractive potential in Gaussian form, V ( x ) = − exp( − x ), and welimit our consideration to the case of zero total angular momentum l = 0. Nevertheless, all the conclusionsare applicable to higher angular momenta as well. The right panel of Figure 7 shows the error in thescattering phase shift of the considered problem calculated with the original JM-ECS and the new approach.To obtain the reference phase shift we use the highly accurate variable phase approach (VPA). We see thatfor all energies the new method gives much more accurate and less oscillatory results. The convergence as afunction of the number of oscillator states in the inner region of both methods is shown on Figure 8. We seethat we arrived at the desired the convergence rate to N − in both low-energy and high-energy regions.16 −3 −2 −1 N osc E rr o r , deg N −1 −2 −1 N osc E rr o r , deg N −1 N −1/2 −4 −3 −2 N osc E rr o r , deg N −1/2 N −1 −4 −3 −2 N osc E rr o r , deg N −1 N −1/2 Figure 9: Same as Figure 8 but for the Morse potential (two upper panels) and Yukawa potential (two lower panels)
The Gaussian shape of the interaction potential was chosen as the main model problem as it is particularlywell adapted to the harmonic oscillator basis and J -matrix method for this potential converges for anyoscillator length b . However we can also test our approach against other short-range potentials. But ifthe asymptotic behavior of the potential is not Gaussian then the results will converge only for specificvalues of the oscillator length, that match the range of the potential. Also, if the potential has a specialpoint in the origin (like Yukawa potential) then the convergence will be much slower due to the importanceof the Fourier contribution in the potential matrix elements. For additional convergence tests we havechosen two potentials with non Gaussian asymptotics and different behavior in the origin: Morse potential V ( x ) = exp( − | x/b | ) − exp( −| x/b | ) and Yukawa potential V ( x ) = − exp( −| x/b | ) /x . In both potentials thepotential range is chosen to match the oscillator length of the basis. Figure 9 shows the convergence of thephase shift for these two potentials. We see that in high-energy region the convergence pattern is mostlysimilar for all potentials as the kinetic energy is much higher than the potential energy in this case. Thebehavior of the error is much less clear for low energies. The convergence rate appears to be the same herefor both methods with Morse potential, for Yukawa potential the convergence rate is blurred due to highoscillations which clearly indicates the importance of Fourier terms. Though we can generally conclude thatthe new approach always gives better and generally less oscillatory result.17 . Extension to systems with more particles The results of the previous sections can be generalized to systems of three and more particles. Let r and r be two relative coordinates describing a three–body problem. The 6D wave function can be expanded inpartial wavesΨ( r , r ) = Ψ( ρ , ρ , θ , θ , φ , φ ) = (cid:88) l ,m ,l ,m ψ l ,l ,m ,m ( ρ , ρ ) ρ ρ Y l o ,m ( θ , φ ) Y l ,m ( θ , φ ) . (49)Such an expansion is used for example to describe double ionization processes in atomic physics [39]Let ψ ( x, y ) now be an infinitely differentiable two–dimensional radial scattering wave function that is ex-panded in the bi-oscillator basis as ψ ( x, y ) = ∞ (cid:88) n =0 ∞ (cid:88) m =0 c nm ϕ n ( x ) ϕ m ( y ) , (50)where x and y should be interpreted as two radial coordinates. The expansion coefficient is then calculatedas a double integral that can be approximated by applying the asymptotic relation twice. First in the x -direction and then in the y -direction c nm = (cid:90) ∞ (cid:90) ∞ ϕ n ( x ) ϕ m ( y ) ψ ( x, y ) dxdy = (cid:90) ∞ ϕ n ( x ) (cid:104) √ R − / m ψ ( x, R m ) + ( b / R − / m ∂ yyy ψ ( x, R m ) + O ( R − / m ) (cid:105) dx =2 ( R m R n ) − / ψ ( R n , R m ) + b √ R − / n R − / m ∂ xxx ψ ( R n , R m ) + b √ R − / n R − / m ψ yyy ( R n , R m )+ b R − / n R − / m ∂ xxx,yyy ψ ( R n , R m ) + O ( R − / n R − / m ) + O ( R − / m R − / n ) . (51)It is important to note that the accuracy depends on both indices n and m . For example, when the low orderapproximation is used for both integrals, there will be error terms of the order R − / m R / n and R − / n R / m .Along the diagonal, i.e. when R m = R n , it is accurate up to order R − / n R − / m = R − n = O ( n − ), which is ofthe order h n , with h n the distance between the classical turning points. For the higher–order approximationthe error term along the diagonal, where R n = R m , will be R − n = O ( n − / ). In the left panel of Figure 10we compare the exact expansion coefficients with the first order and second order approximation along thediagonal for n = m . The results show an improvement with the higher–order expression over the first orderapproximation.However, when one of the indices remains small, for example m , then the error becomes asymptotically O ( R − / n ), as expected from (51). As n becomes large, all error terms with a higher–order accuracy decayuntil the term O ( R − / n R − / m ) dominates for the higher–order formula, or O ( R − / n R − / m ) for the low orderaccuracy. This is shown in the right panel of Figure 10. In a similar way as in the previous sections we can use the asymptotic formulate to build a hybrid repre-sentation that can solve the radial equation of two radial variables that arises after an expansion of a fullthree–body scattering problem equation in spherical harmonics (49). This typically leads to a set of coupledtwo–dimensional radial partial waves. A diagonal block of this coupled set reads (cid:18) − d dx − d dy + l ( l + 1)2 x + l ( l + 1)2 y + V ( x, y ) − E (cid:19) ψ ( x, y ) = χ ( x, y ) (52)18 −7 −6 −5 −4 −3 −2 n,m E rr o r n −1 n −3/2 −5 −4 −3 n E rr o r n −1/2 Figure 10: Convergence of the asymptotic formula for two-dimensional radial functions. The error is calculated by comparingthe exact expansion coefficient with the asymptotic approximations. The function f ( x, y ) = sin( kx ) sin( ky ) along the diagonal(n=m, left) and along the line with fixed m=20 (right). Crosses show the first order and solid lines correspond to the secondorder asymptotic expansion coefficient ( k = 1, b = 1). with boundary conditions ψ ( x,
0) = ψ (0 , y ) = 0.The solution of this equation can be represented in a bi-oscillator basis, (50), using two indices n and m .Then similar to the two–body case we use approximate oscillator coefficients in the asymptotic region. Butthere are now three asymptotic regions: first, the region where n is large, second, the region where m islarge and finally, the region where both n and m are large. For each of these regions we define approximateoscillator coefficients that take the form: c nm ≈ ` c nm := (cid:88) i U ni (cid:90) ∞ ϕ m ( y ) ψ ( R i , y ) dy + O ( R − / n ) where n (cid:29) c nm ≈ ´ c nm := (cid:88) i U mi (cid:90) ∞ ϕ n ( x ) ψ ( x, R i ) dx + O ( R − / m ) where m (cid:29) c nm ≈ ˇ c nm := (cid:88) i,j U ni U mj ψ ( R i , R j ) + O ( R − / n R − / m ) + O ( R − / m R − / n ) where n, m (cid:29) , (55)where the error terms in the last expression are explained in (51).Then the coefficient matrix c nm of hybrid representation is divided into four blocks corresponding to differentregions c nm = c . . . c M − ´ c M . . . ´ c K (cid:48) c . . . c M − ´ c M . . . ´ c K (cid:48) ... ... ... ... c N − . . . c N − M − ´ c N − M . . . ´ c N − K (cid:48) ` c N . . . ` c N M − ˇ c N M . . . ˇ c N K (cid:48) ` c N +1 0 . . . ` c N +1 M − ˇ c N +1 M . . . ˇ c N +1 K (cid:48) ... ... ... ...` c K . . . ` c K M − ˇ c K M . . . ˇ c K K (cid:48) , (56)where we define N and M as sizes of the oscillator bases in each direction and K and K (cid:48) are the totalnumber of variables in each direction. 19o build a linear system corresponding to (52), we reshape the matrix c nm to a vector. Then the Hamiltonianmatrix of a 2D problem is constructed as a Kronecker sum of the two-body Hamiltonian’s constructed asin (47) and a two-body potential matrix. In the case of our approximate oscillator representation this takesthe form ˇ H ( osc )2 D = ˇ H ( osc )1 ⊗ ( U y W y ) + ( U x W x ) ⊗ ˇ H ( osc )2 + ( U x ⊗ U y ) V ( fd )12 ( W x ⊗ W y ) (57)so the Kronecker sum contains the approximated unity operator instead of the real one. Where U x and W x are the transformation matrices for the x coordinate and similarly for the y coordinate. The total size of theHamiltonian matrix is K × K (cid:48) . We first evaluate the method with the help of a model Helmholtz problem with a constant wave number orenergy E = k / l = 0 and l = 0. The equation is then (cid:18) −
12 ∆ − E (cid:19) ψ ( x, y ) = ϕ ( x ) ϕ ( y ) , (58)with boundary conditions ψ ( x,
0) = 0 and ψ (0 , y ) = 0. The form of the right–hand side here was chosenfor simplicity as it results in a vector (1 , , , . . .
0) in bi-oscillator representation. This problem is exactlysolvable with the help of the Greens function G ( x, y ; x (cid:48) , y (cid:48) ) = i (cid:16) H (1)0 ( k (cid:112) ( x − x (cid:48) ) + ( y − y (cid:48) ) ) − H (1)0 ( k (cid:112) ( x + x (cid:48) ) + ( y − y (cid:48) ) ) (59) − H (1)0 ( k (cid:112) ( x − x (cid:48) ) + ( y + y (cid:48) ) ) + H (1)0 ( k (cid:112) ( x + x (cid:48) ) + ( y + y (cid:48) ) ) (cid:17) , (60)where H (1)0 is a 0-th order Hankel function of the first kind. The scattering solution is then ψ ( x, y ) = (cid:90) ∞ (cid:90) ∞ ϕ ( x (cid:48) ) ϕ ( y (cid:48) ) G ( x, y ; x (cid:48) , y (cid:48) ) x (cid:48) y (cid:48) dx (cid:48) dy (cid:48) (61)For simplicity we are not going to compare any scattering information extracted from the wave function asit involves additional operations like surface integration, which can lead to additional loss of accuracy andneed to be studied separately. We will compare the values of the wave function in a fixed spatial point. Asthe spatial grid is different for every size of the oscillator basis we can not ensure that any fixed spatial pointlies exactly on the grid point at every calculation, so we need to interpolate. We have used a cubic splineinterpolation, and from comparison to the other interpolation methods we can expect that the additionalerror introduced by this operation is negligible compared to the errors of the method (of course, that is onlyif the function is smooth enough, that means not very high values of E ).As we do not have any potential in this model problem, the convergence of the proposed hybrid oscillatorrepresentation will depend mainly on the accuracy of the asymptotic relations used. This means that wecan choose the size of the spatial domain on which we solve the two–dimensional radial problem as small aspossible but not smaller than the region spanned by the exact oscillator basis. As the maximal size of theoscillator basis we use in this calculations is 350 and the oscillator length was chosen as b = 0 . R N ≈ ,
28) when we increasethe basis size simultaneously and at (28 ,
12) when we keep the basis size fixed in the y -dimension.In Figure 11 we compare the numerical solution of the linear system with the one obtained from (61) whichwe consider as exact. We see that, similar to Figure 10, if we expand the basis only in one direction theconvergence is of the order N − / . 20 −4 −3 −2 N,M E rr o r N −1 −4 −3 −2 N E rr o r N −1/2 Figure 11: Convergence of the two-dimensional scattering wave function for the problem in Eq. (58) at the fixed spatial point.Left figure represents the case where the basis was increased in two dimensions simultaneously ( N = M ), right figure shows theconvergence with a fixed M = 50. Both calculations are made with b = 0 . E = 2. Next we test the method on a model potential scattering problem described by (52) with zero partial angularmomenta and the Gaussian interaction potential in the form V ( x, y ) = − − x − − y + 10e − x − y (62)The one-body potential, here V ( x ) = − − x ), can support one bound state. This means that usingthe potential (62) we can model a three-body breakup problem using the right–hand side in the form χ ( x, y ) = − ( V ( x, y ) − V ( x )) f ( x )ˆ j ( y ) , (63)where f ( x ) is the wave function of the bound state, ˆ j ( y ) is the Riccati-Bessel function and together theyrepresent the initial state of the system. For the considered problem there are no analytical results tocompare with, so on the Figure 12 we present the solution of the linear system for c nm and the spatial wavefunction, reconstructed from c nm by applying the transformation matrix W ⊗ W . On both figures we canclearly recognize the patterns of elastic scattering (the plane waves along the axes) and the breakup (theradial waves with lower frequency as part of the energy was spent to break the bound state). The in-depthanalysis of the three–body scattering results is the subject of future work.
6. Conclusions and outlook
This paper focuses on scattering calculations in the oscillator representation, where the solution is expandedin the eigenstates of the harmonic oscillator. The oscillator representation is not the most natural represen-tation to describe scattering processes since it involves a basis set that is designed to describe bound states.This leads to a rather low convergence and may result in a linear system with very large dense matrix. It isoften more natural to describe a scattering process with the help of a grid representation. These grid–basedcalculations have been very successful in describing scattering and breakup processes in atomic and molecularphysics. The Helmholtz equation is also efficiently solved on these grids.However, internal structure of the nuclear clusters and other many–particle systems are efficiently describedby such an oscillator basis, since the bottom of the potential can be mimicked by the oscillator potentialleading to an efficient description of the internal structure.21 igure 12: The solution of three–body scattering problem in oscillator representation (left figure) and the reconstructed spatialwave function (right figure). The calculation is made with b = 1 and E = 1. The problem has Dirichlet boundaries on all sides. In this paper we combine the advantages of grid–based calculations with those of the oscillator representation.The method was originally proposed in [33], as a further development of the J -matrix or algebraic method forscattering. There the method combined the grid and oscillator representation with a low–order asymptoticformula. In this paper we have improved this matching with a higher–order approximation and this bringsthe overall error down to the level of the discretization error of the finite–difference grid.Although a similar asymptotic formula appeared earlier in the work of S. Igashov [34], we believe that theasymptotic formula presented in the current paper is more generally applicable. Furthermore, we have usedthe asymptotic formulas to improve the accuracy and convergence of the hybrid simulation method. Theconvergence of the method is illustrated with various examples from two–body and three–body scattering.In preparation of this papers, our initial efforts involved application of the strategy showed in Figure 2 withthe higher–order formula. However, this required the use of a forward and backward stencils to estimate thethird derivative in the first grid point of the finite difference representation. This strategy, however, gives riseto eigen modes with a negative energy localized near the interface. This destroys the positive definitenessof the Laplace operator. These modes are avoided in the proposed method that uses a symmetric stencilaround the matching point. This gives satisfactory convergence.We have applied the asymptotic technique to radial oscillator state that are based on Laguerre polynomials.Similar results can be derived for the 1D oscillator states that are based on the Hermite polynomials.These Hermite polynomials are closely related to Gaussian Type Orbitals (GTO) that are frequently used incomputational quantum chemistry [40] . It is worth to explore if the asymptotic formulas make it possibleto combine GTO’s with grid based calculations to describe molecular scattering processes.In the future further research is necessary on asymptotic expressions of products of two functions. This willinvolve convolution integrals. Better expressions for products will also help to take into account asymptoticbehavior of the potentials in the grid representation.
7. Acknowledgments