Coupled-channel continuum eigenchannel basis

The goal of this paper is to calculate bound, resonant and scattering states in the coupled-channel formalism without relying on the boundary conditions at large distances. The coupled-channel solution is expanded in eigenchannel bases i.e. in eigenfunctions of diagonal Hamiltonians. Each eigenchannel basis may include discrete and discretized continuum (real or complex energy) single particle states. The coupled-channel solutions are computed through diagonalization in these bases. The method is applied to a few two-channels problems. The exact bound spectrum of the Poeschl-Teller potential is well described by using a basis of real energy continuum states. For deuteron described by Reid potential, the experimental energy and the S and D contents of the wave function are reproduced in the asymptotic limit of the cutoff energy. For the Noro-Taylor potential resonant state energy is well reproduced by using the complex energy Berggren basis. It is found that the expansion of the coupled-channel wave function in these eigenchannel bases require less computational efforts than the use of any other basis. The solutions are stable and converge as the cutoff energy increases.


Introduction
Considerable amount of efforts are devoted all around the world for studying the properties of unstable nuclei [1]. Because of this, new theoretical approaches, which takes into account the continuum explicitly, is called for revealing their properties. The coupled-channel method is a very powerful formalism for studying the structure of both strongly-bound nuclei [2,3] and loosely-bound nuclei [4] too. Here we propose a way to calculate the coupled-channel solutions in which all bound and continuum (resonant and non-resonant continuum) states are treated on the equal footing.
Complex eigenenergies, i.e. Gamow [5] or Siegert [6] states were calculated using the Green's function approach in momentum space in Refs. [7,8] for coupled channel problems. Gamow states for realistic deformed potentials were calculated first in Ref. [9] by solving the logarithmic derivative of the coupled equations with outgoing boundary condition. In Refs. [10] and [11] the coupled-channel Schrödinger equation with outgoing wave boundary condition were used to study the proton decay states in a rare-earth nucleus.
The complex scaling method has been successfully combined with the coupled-equation formalism to calculate resonances [12,13,14,15]. The extension of the Gamow Shell Model [16,17] to reaction problems in the framework of coupled-channel formalism was recently implemented in Ref. [18], where the low-lying states of 7 Li were calculated. The result of the direct integration of coupled equations was compared with that of the Berggren [19] expansion for the calculation of bound states of dipolar molecules in Ref. [20]. A full complex energy representation, was used in Ref. [21] for the calculation of the Isobaric Analog State by coupled Lane equations. The present paper extends the use of the continuum bases to the inelastic processes in coupledequation systems and to the calculation of scattering states.
The method presented in this paper allows the calculation of bound, resonances and scattering states in coupled systems on the same footing. All these states may be found by a single diagonalization simultaneously. Each channel wave function is expanded in an optimized basis set defined by the eigensolution of the corresponding uncoupled Schrödinger equation, that is the meaning of eigenchannel bases. Since this method prescinds from explicit boundary conditions, it might be useful for dealing with Coulomb breakup problems that appear, for instance in electron-impact ionization [22] or in breakup reactions important in astrophysics [23,24] or in studying the threebody Coulomb breakup reaction of 11 Li [25].
In section 2 we develop the method in which the coupled Schrödinger equations are expanded in the continuum bases of uncoupled channels. The first application of the method is done in section 3. It solves the problem of the exactly solvable two-channel Poeschl-Teller potential. This works as a test case. It shows the reliability of the method and it shows the relative importance of the continuum for the deep and for the loosely bound states. In Section 4, the method is applied to the bound and scattering states of the deuteron. The last application in section 5 is devoted to the simultaneous calculation of bound and resonant states. The outline for next applications and some remarks are given in the last section 6.

Formalism
Let us denote by H the Hamiltonian which describes a collision between two nuclei being in bound states (a, A). We split H into two parts: (1) the Hamiltonian H ′ α that is left when the two initial fragments are far away from each other and (2) V = i∈a,j∈A V ij which includes the projectile(a)target(A) interaction. Changing in H ′ α to relative coordinates in each fragments and then changing to the relative coordinates between the fragments [26], we end up with H ′ α = H α + T (we have set to zero the centroid kinetic energy), where T = − 2 2µ ∇ 2 r is the relative kinetic energy, µ is the projectile-target reduced mass, and H α = H a + H A , where H a and H A are the intrinsic Hamiltonians of the projectile and target, respectively. Then, the total Hamiltonian reads H = H α + T + V . The residual interaction V = V d + V od is split into a diagonal part V d and an off-diagonal one V od [27]. The eigenfunction ψ J π M of H is expanded into different channels using the channel basis functions Φ J π M α defined as Then, Substituting the channel expansion (2) into the Schrödinger equation H ψ J π M (r, a, A) = E ψ J π M (r, a, A) and projecting into a certain channel Φ α we get (omitting the index J π M) where we have separated the diagonal matrix elements V αα and we have defined the single particle channel Hamiltonians, with V αα ′ = Φ α |V |Φ α ′ raA , where the suffix indexes mean integration over the angular coordinater of the relative motion and the internal coordinates of the projectile a and target A nuclei, respectively. Notice that the structure of the Eqs. (3) and (4) are the same to that of Eq. (25) of Ref. [27]. Although, in principle any complete set of states will allow the computation of the interaction matrix elements, in practice, a judicious choice of the basis states will minimize the number of matrix elements to be calculated and reduce the computation time needed. Here we use the diagonal part V d of the residual interaction V = V d +V od to generate the basis. Notice that the basis does not correspond to the one generated without residual interaction V = 0.
In the next step, we expand the wave functions u α (r) in each channel in the basis generated by its own channel Hamiltonian h α h α u (0) α,n (r) = ε (0) α,n u (0) α,n (r) (5) where the summation includes integration over the continuum part of the spectrum of h α . Replacing the expansion of u α (r) (Eq. (6)) in Eq. (3) and projecting over u where N denotes the number of channels and M α is the number of single particle basis states for the channel α.
The matrix elements between different channels contain only the interaction V ii ′ given by Using the basis generated by the diagonal part of the channel interaction one can save the calculation of interaction matrix elements. The number of these matrix elements increases rapidly as the number of open channels N and the dimension of the basis M α increases.
There are two advantages of using a basis expansion method instead of using the asymptotic boundary conditions. The matrix diagonalization does not diverge even if the coupling terms are large. This might happen in the direct numerical integration [3] of the coupled equations. The matrix diagonalization does not face any instability of the numerical integration of the coupled equations. The disadvantage of using a basis expansion is that one has to deal with the completeness problem of the basis. A difficulty of using basis expansion is that one needs an efficient and accurate method to solve the single particle Schrödinger equation, that is, to find real and complex poles as well as the real and complex energy scattering states. The real and complex energy scattering states were calculated by using a piecewise perturbation method [28]. The code implements the so called Ixaru's method [29]. The real and complex energy poles were also calculated by using a modified version of the program [28]. This version has a higher precision than the GAMOW code [30], which however is more flexible.

Application to the Poeschl-Teller potential: bound state calculation using bases composed of bound states and real energy continuum
In this section we compare the exact solution of the two-channel Poeschl-Teller potential with the numerical solution using the same eigenchannel bases for both channels. The bases are composed of bound and real energy scattering states.
Let us consider the Schrödinger equations with two-channels and with = 2µ = 1, l 1 = l 2 = 0, ε 1 = ε 2 = 0 and V αα ′ (r) given by Ixaru [31] For this special interaction the coupled equations (3) can be collected into two uncoupled equations, . Then, one may choose V d and V od such that V ± have exact solutions. In this way we find the eigenvalues E + and E − of h + and h − which will be also eigenvalues of the original coupled-equation, i.e.

Basis expansion
The channel wave functions u 1 (r) and u 2 (r) in Eq. (3) are expanded in the same basis, since (11), is shown if Fig. 1. The basis is formed by the five bound states u 5 = −1.25622, while the continuum is discretized using Gauss-Legendre partition ε i ∈ (0, ε max ) with weights ω i .
The matrix elements V αn,α ′ n ′ (with α, α ′ = 1, 2 and n, n ′ = 1, . . . , 5 + N c ) were calculated using the potential V od (r) Eq. (12) shown in fig. 1. The integration was performed by using Gauss-Legendre quadrature with 40 mesh points between r = (0, 10). The convergence of the solutions was studied as a function of the cutoff energy ε max and the number of mesh points N c . Table 1 shows the convergence of the energies for ε max = 70 as the function of the number of continuum states N C . We can see a fast convergence for all states except the state being closest to the threshold. It is worthwhile to mention that all ten perturbed bound states were found in a single diagonalization by using bases with only five (unperturbed) bound states.
4. Application to the proton-neutron system: bound state and eigenphase shifts calculations using real energy continuum bases In this section we solve the coupled-equation for the deuteron, and calculate the eigenphases δ S and δ D using the soft core potential of Ref. [32]. Only a single adjustable parameter, the cutoff energy is used here.
For this system the quantities which appear in the coupled equations (3) are: l 1 = 0, l 2 = 2, ε 1 = ε 2 = 0, 1 µ = 1 mp + 1 mn (with m n and m p the neutron and proton mass, respectively), and 2µ/ 2 = 0.0241138 (MeV.f m 2 ) −1 . The potentials (see fig. 2) are given by the following expressions, Table 1: Two-channels Poeschl-Teller energies obtained in diagonalization as function of the number N c of continuum states in the bases. The cutoff energy is ε max = 70. E exact refers to the exact energies of Eqs. (9) and (10) The radial coordinate r is given in fm while the interactions are given in MeV units and x = (0.7f m −1 ) × r is dimensionless. For ground state of the deuteron, the channel wave functions u 1 (r) ≡ u S (r) and u 2 (r) ≡ u D (r) correspond to the 3 S 1 and 3 D 1 components of the wave function. While for the n − p scattering states u 1 (r) and u 2 (r) are standing-waves [33] which asymptotically behave like u 1,α (r) ∼ sin(kr + δ α ) and u 2,α (r) ∼ sin(kr − π + δ α ) with α = S, D and δ S and δ D the eigenphase shifts [34,35].

Basis expansion
Since none of the diagonal potentials hold any bound state (the potential V 11 has an anti-bound state at the energy −5.671 MeV), the two bases are formed from discretized continuum states Here ω i are the weights of the Gauss-Legendre mesh points at the energies ε i ∈ (0, ε max ). We took the same continuum basis for both channels, i.e. both channels were expanded using the same number of mesh points up to the same cutoff energy ε max .
The interaction matrix elements were calculated using Gauss-Legendre quadrature with r ∈ (0, r max ). It was checked that the bound state solution was stable when we varied the cutoff radius between 16 fm to 24 fm. For the calculation we took 20 fm for the cutoff radius and 100 for the number of mesh points.
The ground state energy E d and wave function of the deuteron were calculated as a function of the cutoff energy ε max . The real energy basis states were defined by the following vertices (in MeV): (0, 10), (10, 50), (50, 100), (100, 250), (250, 500), (500, 750), (750, 1000), (1000, 2000),. . . (9000, 10000), where . . . means that an interval of 1000 MeV have been used. Six mesh points for each interval up to 1000 MeV and ten mesh points from there on, have been taken. While the same energy partition was used for both channels the scattering functions were not the same since h 1 and h 2 were different. Figure 3 shows that the deuteron ground state energy converges very slowly to the experimental value E exp = −2.224 MeV as the cutoff energy increases. The final value E d = −2.210 MeV gives 93.7 % and 6.3 % for the 3 S 1 and 3 D 1 partial wave amplitudes, respectively. Both, the energy and the wave function content fit well to the experimental values.  The ground state wave function was built from the eigenvector of E d , Eq. (6). Fig. 4 shows the convergence of the channel wave functions as the cutoff energy increases. The D component of the deuteron wave function shows oscillations due to the oscillations in the high energy basis states. The magnitude of the oscillations decreases by increasing the mesh points energies or by increasing the cutoff energy. In short, the oscillations is due to the incomplete representation used in the expansion of Eq. (6). The oscillations in the S component are much smaller and they are not visible at the scale used in the figure. The difference in the magnitudes of the oscillations between the two components of the wave function could be attributed to the large differences between the values of the diagonal potentials in the tail region (see fig. 2). From the diagonalization we obtained, beside the ground state, the scattering states. As for the bound ground state, we get the corresponding eigenvectors for each positive eigenvalue. Then, Eq. (6) gives the scattering states expanded in the continuum basis. These scattering solutions can be used to calculate the eigenphase shifts [34,35] δ S and δ D . From each scattering state we built the sum of the two channels wave functions and fit it to the function A[cos(δ) * F 0 (kr) + sin(δ) * G 0 (kr)] + B[cos(δ) * F 2 (kr) + sin(δ) * G 2 (kr)], for r > 9 fm. The functions F l (kr) and G l (kr) are the regular and irregular Coulomb functions, respectively. They were calculated using the program [36]. The eigenphases were calculated using the Levenberg-Marquad code from Numerical Recipes [37]. Using the same basis that for the deuteron we calculated the eigenphase shifts δ S and δ D and compared them with the results obtained by doubling the mesh points for energy larger than 1000 MeV; it was found that the values of δ S changed around 1% while the values of δ D changed in the third figure. Then, we doubled and tripled the mesh for the energies below 1000 MeV. There, it was found that the eigenphase δ D has a smooth behavior, while the eigenphase δ S shows oscillations around and above 100 MeV. In order to have a more uniform distribution, we made intervals of 10 MeV from zero up to 350 MeV and took the same number of point n c in each interval. We calculated the eigenphases for increasing n c . For n c = 1, 2, 3, we found the same qualitative behavior than for the deuteron basis, while for n c = 4 almost no δ S was found for energies larger that 120 MeV within the required error. Instead δ D was fitted all the range except between 290 MeV and 310 MeV. Using these last basis the fit was done with the restriction that for energies larger than 100 MeV only the amplitude of the most important channel together with the eigenphase shift was fitted. Figure 5 shows the results of these two calculations. We can notice that δ D values resulted by both calculations are very similar, while the two parameters fit calculation smoothly connects to the calculation of δ S using the three parameters fit. 5. Application to the Noro-Taylor potential: resonant state calculations using complex energy basis states The Noro-Taylor potential [39] is a two-channel model with a strong repulsion in the second channel and with a strong coupling. This system has a narrow resonance at the energy E r = 4.7682 (Γ = 0.001420). The parameters (in atomic units) which appear in the coupled equations (3) are: l 1 = l 2 = 0, ε 1 = 0, ε 2 = 0.1, µ = = 1. The potentials are given by the following expressions,

Basis expansion
The diagonal potential V 11 has two bound states at energies ε The matrix elements V 1n,2n ′ were calculated using Gauss-Legendre quadrature with 100 mesh points for the radial coordinate from 0 to 35.
Since the potential V 22 is repulsive the channel wave function u 2 is expanded only by continuum basis states. We considered two different bases: (1) only real positive energy states are included, which we call real energy representation, and (2) complex energy states are also included, which we call complex energy (Berggren) representation. Using the real representation only bound states can be found while using a properly chosen complex energy representation we can get also resonant states.
The first channel wave function u 1 is expanded in terms of the wave functions of the two bound states ε (0) 1,1 and ε (0) 1,2 plus a set of discretized continuum states along the real axis up to a cutoff energy ε 1,max . The second channel wave function u 2 can either be expanded by using a set of discretized real energy scattering states up to an cutoff energy ε 2,max (real energy representation) or alternatively by the resonant state ε (0) 2,1 plus a set of discretized complex energy scattering states taken along a contour in the complex energy plane.
The resonance in the second channel affects the selection of the mesh points since its presence requires a denser mesh of the scattering states in the vicinity of the resonant energy. The second column of table 2 shows the contour used in the real energy representation (RR), while the fourth column shows the contour for the complex energy representation (CR). Notice that the real representation requires a larger ε 2,max .  (50,0) 50 Table 3 shows the perturbed energies calculated using either the real or the complex energy representations. For comparison we give the results obtained in Ref. [40] using the Jost function method combined with complex energy rotation, this corresponds to the last column (under the title E exact ). It is found that only the complex energy representation (Berggren basis [19]) is able to reproduce simultaneously bound and unbound perturbed states in this coupled channels system. Table 3: First five poles of the Noro-Taylor potential calculated from the diagonalization in real (E RR ) and complex (E CR ) representations compared to E exact given in Ref. [40].

Conclusions
Loosely and deeply bound states, resonant states, and eigenphase shifts have been calculated in an optimized continuum basis.
Each channel defines its own basis (eigenchannel basis) through the diagonal parts of the channel potentials, in this way the number of matrix elements to be calculated is reduced considerably.
The matrix diagonalization gives the poles and scattering solutions simultaneously, i.e. the basis is energy-independent.
Since the solutions of the coupled-channel equations do not rely on boundary conditions, the method could be convenient for studying systems where the boundary conditions cannot be treated easily.
In summary, we presented a method for describing nuclear reactions involving weakly bound or unbound system. The next step is to apply the continuum eigenchannel basis expansion for studying deuteron elastic breakup process.
The author thanks Prof. T. Vertse for valuable discussions. This work has been partially supported by the National Council of Research PIP-77 (CONICET, Argentina).