2 Solution of coupled hyperradial equations

For coupled equations arising for A independently-moving bodies, at total energy E we have to solve a set of N coupled equations in a radial variable r for channels labelled $i = \{K,\gamma\}$:

$$\Bigl( - \frac{\hbar^2}{2m} \Bigl[ \frac{d^2}{dr^2} - \frac{{\cal... ...psilon_i - E \Bigr) \psi_{ij}(r) = - \sum_{i' \not= i} V_{ii'}(r) \psi_{i'j}(r)$$ (1)

where ${\cal L}_i = K+\Delta$ with $\Delta = (3A-6)/2$, and $\epsilon_i$ is the sum of the internal energies of the A bodies in channel i. There are sometimes orthogonality conditions $\langle u_m \vert \psi \rangle = 0$ where u_m is (for example) from the wave function of occupied states in some binary subsystem.

For two-body problems, r is the separation of the bodies and m is the reduced mass, while for three- or more-body problems this equation arises when r is the hyperradius $\rho$, and m is a unit scaling mass.

The $\gamma$ refers to a set of additional quantum numbers. For three-body problems with inert particles and a spin-zero core, $\gamma = \{ l_x, l_y, L, S \}$, while in general the quantum numbers of the A interacting bodies will be included in the set $\gamma$.

We use the channel momenta $k_i^2 = 2m/\hbar^2~ (E-\epsilon_i)$. The radial wave functions $\psi _{i }(r )$ for bound states then have the standard boundary conditions of

$$\psi_{i}(r \rightarrow 0) \sim r^{{\cal L}_i+1} \mbox{ ; } \... ... \psi_{i }(r \rightarrow \infty ) \sim \exp(-\vert k_i\vert r)$$ (2)

and for continuum wave functions the conditions as $r \rightarrow \infty$ of

$$\psi_{ij}(r ) \sim H^{-}_{{\cal L}_i}(k_i r)\delta_{ij} - H^{+}_{{\cal L}_i}(k_i r )~ S_{ij}\ .$$ (3)

Here $H^{-}_{\cal L}$and $H^{+}_{\cal L}$ are the Coulomb functions of index ${\cal L}$, with asymptotics $\sim \exp(\mp ik_i r )$ describing the in- and out-going spherical waves. The appropriate value of the Coulomb parameter $\eta$ is given in section 2.5. The S_ij is the S-matrix for the outgoing amplitude in channel i from an incoming plane wave in channel j. We only deal with Borromean or democratic types of many-body problems, so we do not treat explicitly any sub-system bound states in the asymptotic regions.

Note that for zero potentials S=1, and the scattering wave functions become

$$\psi_{ij}( r ) \sim \frac{2}{i} F_{{\cal L}_i} (k_i r) \ ,$$ (4)

since $H^\pm = G\pm iF$. These are the plane wave parts of all solutions $\psi_{ij}( r )$.

2.1 Previous methods

To solve the coupled equations (1), we previously [1,4,5,3] integrated N linearly independent solutions from r=0 to r=r_m, for some radius r_m beyond which the couplings V_ii'(r) are assumed to be negligible. A linear combination of these solutions was then found, in order to satisfy the boundary conditions (3). This was straightforward for limited numbers of channels (N up to 10 or 12), but for larger sets the solutions become linearly dependent. This is because the large range of centrifugal barriers ${\cal L}({\cal L} +1)/r^2$ means that some channels are, for a large period, being integrated in their classically forbidden domain, and (as discussed in [6]) the linear independence of the exponentially decreasing solutions is lost.

2.2 Diagonalisation methods

To avoid numerical instabilities with large sets, we use an expansion on Sturmian or Sturmian-like radial basis functions. For bound-state wave functions, it is convenient to use Sturmian states, which are eigenstates at a fixed energy for varying multiples of some diagonal potential. This has been developed for bound states in deformed nuclei [7,8,9], and we have also found them to be useful for hyperradial bound states. For scattering states, however, the Sturmian states are not suitable as they have few oscillations at large distances up to r_m. We therefore use a basis set of `energy eigenstates' of the diagonal terms of eqns. (1):

$$\Bigl( - \frac{\hbar^2}{2m} \Bigl[ \frac{d^2}{dr^2} - \frac{{\cal L}({\cal L} +1)}{r^2} \Bigr] + V_{ii}(r) - \alpha_q \Bigr) f^q_{i}(r) = 0$$ (5)

for eigenenergies $\alpha_q$, with the basis functions all having fixed logarithmic derivatives $\beta = d\ln f^q_{i}(r)/dr$ at r_m (rmax). The constancy of the logarithmic derivatives $\beta$ means that (for each i channel separately) the f^q_i form an orthogonal basis set over the interval [0,r_m], and over this range they can be normalised to unity. Then, as many q=1...Q of the orthonormal f^q_i per i channel can be used as desired for accuracy. For $\beta < 0$ the basis states with low $\alpha_q<0$ will be similar to the Sturmian states with low potential multipliers, and hence suitable for bound state expansions. Higher $\alpha_q > 0$ basis states will oscillate out to r_m, and can be used in expansions of continuum states.

The wave functions of the coupled problem (1) can now be solved completely over the interior range [0,r_m], by using the orthonormal basis set of the {f^q_i(r)} with coefficients to be determined. The coefficients are found in two stages: first by finding all the eigensolutions $\phi^p_{i}(r)$ of the equations (1) using the above orthonormal basis, and then expanding the scattering wave functions in terms of these $\phi^p_{i}(r)$.

The first stage, the diagonalisation of interior Schrödinger equation (1) yields P=QN eigenenergies e_p with corresponding multichannel eigenstates

$$\phi^p_{i}(r) = \sum_{q=1}^Q ~ c^{pq}_{i} f^q_{i}(r)$$ (6)

Eigenstates here with e_p<0 are the bound states, while solutions with e_p > 0 contribute to the scattering solutions. Certain of the e_p > 0 solutions may correspond to low-lying resonances if those are present, but the majority of the positive eigenenergies have no simple physical interpretation. These $\phi^p_{i}(r)$ form of course another orthonormal basis in the interior region.

For scattering states at arbitrary energy E, the coupled solutions are then expanded in terms of the multichannel eigenstates as

$$\psi _{ij} = \sum_p A^p_{ij} \phi^p_{i} \ .$$ (7)

We define an R-matrix at energy E by

$$\psi _{ij}(r) = \sum_{i'} r R_{ii'}(E) \Bigl [ \psi{'} _{i'}(r) - \beta \psi _{i'j}(r) \Bigr] \ .$$ (8)

In the limit of $r\rightarrow r_m$ from above, the R-matrix $\bf R$ can be calculated directly from the eigenstates by standard methods [10]

$$R_{ii'}(E)= \frac{\hbar^2}{2mr_m} \sum_{p=1}^P \frac{ \phi^p_{i}(r_m) \phi^p_{i'}(r_m) } { e_p - E} \ .$$ (9)

Using eqns. (3) and (8), and writing the Coulomb functions ${\bf H^{\pm}}$ as diagonal matrices, the scattering S-matrix is given in terms of $\bf R$ by

$${\bf S} = [{\bf H^+} - r_m{\bf R}({\bf H'^+} - \beta {\bf H^... ...{-1} [ {\bf H^-} - r_m{\bf R}({\bf H'^-} - \beta {\bf H^-}) ]$$ (10)

and the expansion coefficients for the wave functions are

$$A_{ij}^p = - \frac{\hbar^2}{2m} \frac{1}{e_p-E} \sum_{i'} \phi^p_{i'}(r_m) \Bigl [ \delta_{ij} (H'^-_{\cal L}(k_ir_m)-\beta H^-_{\cal L}(k_ir_m))$$

$$S_{ij} (H'^+_{\cal L}(k_ir_m)-\beta H^+_{\cal L}(k_ir_m)) \Bigr ] \ .$$ - (11)

The coefficients c^pq_i and energies e_p in eqn. (6) satisfy matrix equations

$$\alpha_q c^{pq}_{i} + \sum_{q'i'} \langle f^q_{i} \mid V_{ii'} \mid f^{q'}_{i'} \rangle c^{pq'}_{i'} = e_p c^{pq}_{i}$$ (12)

for each eigenstate p, which are of the matrix form

$${\bf A}c = e c \ .$$ (13)

2.3 Buttle Correction

The R-matrix calculated by eqn. (9) is only exact when the sum over p extends to all energies e_p. To improve the accuracy of calculations with finite Q (and hence finite P), the Buttle correction [19] is added to the right hand side of eqn. (9). This modifies the diagonal terms R_ii(E) to reproduce for each uncoupled problem the exact scattering solution $\chi_{i}(r)$ after this has been integrated separately. Since we use the energy eigenstates f^q_i(r) rather than Sturmian basis states, the R-matrix sum from (9) for each uncoupled channel is simply

$$R^N_{i}(E)= \frac{\hbar^2}{2mr} \sum_{q=1}^Q \frac{ f^q_{i}(r_m)^2 } { \alpha_q - E}$$ (14)

and the exact one-channel R-matrix is

$$R^0_{i}(E) = r^{-1} ~ \chi_{i}(r_m)~(\chi'_{i}(r_m) - \beta \chi_{i}(r_m))^{-1}\ .$$ (15)

The Buttle-corrected full R-matrix to be used in eqn. (10) is then

$$R^c_{ii'}(E) = R_{ii'}(E) + \delta_{ii'} \Bigl [~ R^0_{i}(E_1) - R^N_{i}(E_1)\Bigr ] \ .$$ (16)

The energy E₁ can be equal to E, or chosen just near to it if necessary to avoid the poles in eqn. (14), since the Buttle correction varies smoothly with energy.

2.4 Orthogonality Constraints

This set of equations is easily modified to take into account the requirement from the Pauli Principle of orthogonality to a set of occupied nucleon states $\vert u_m\rangle$ in the core nucleus. Sometimes [3] an l_x-dependent potential is used to approximately produce this orthogonality, by using repulsive potentials in the partial wave channels with occupied states. Now, however, we can also use a more accurate orthogonalisation to specific occupied wave functions. The occupied states $u_m({\bf x})$ are first multiplied by one of a complete set of spline functions $s_n({\bf y})$ (as in [2]) to make a many-body wave function $U_{mn}({\bf x},{\bf y}) = u_m({\bf x}) s_n({\bf y})$. This wave function (for total spin state $\vert JM_J\rangle$) is expanded in the hyperspherical basis f^q_i(r) as

$$U_{mn}({\bf x},{\bf y}) = \sum_{i q} w^q_{mni} f^q_{i}(r) \Upsilon^{l_xl_y}_{JKLSM_J}(\Omega_5) \ ,$$ (17)

where $\gamma = \{ l_x, l_y, L, S \}$ as before. The orthogonality requirement of $\langle u_m({\bf x})\vert \psi_{ij}\rangle = 0$ is now satisfied by requiring that each interior multichannel eigenstate p has no forbidden component:

$$\sum_{qi} w^q_{mni} c^{pq}_{i} = 0$$ (18)

for each forbidden state m, each spline function n. In matrix form, this is $\langle w_{mn} \vert c\rangle = 0$ for each eigenstate p. We therefore construct a projection operator

$$P = \sum_{mn} \vert \tilde{w}_{mn} \rangle \langle \tilde{w}_{mn} \vert\ ,$$ (19)

where the set $\{ \tilde{w}_{mn}\}$ is an orthonormalised basis set constructed by the Gramm-Schmidt process from the { w_mn}.

We now diagonalise in the allowed subspace of P c = 0, and ensure this by replacing matrix equations (13) by

$$(1-P){\bf A}(1-P)c + E_0 P c = e c$$ (20)

where E₀ is a large positive energy to which the forbidden eigenstates are moved. All the solutions of (20) with $e \ne E_0$ satisfy P c = 0. In this way, both the bound state and continuum wave functions can be made orthogonal to the required set of occupied core states.

2.5 Asymptotic Propagations

The radial wavefunction describing the scattering with in the outer region of configuration space, where exchange effects and orthogonality conditions may be neglected, satisfies the set of coupled differential equations Eq. (1), which we now write as

$$-\frac{\hbar^2}{2m} \left (\frac {d\sp 2}{dr\sp 2} - \fra... ...=1}\sp N V_{ii'}(r)\psi_{i'j}(r) = 0, \;\;\; i,j=1, \ldots ,N$$ (21)

We here separate a Coulomb term proportional to Z_i. For two-body scattering of particles with charges Z_a and Z_b, this $Z_i = - Z_aZ_be^2 m/\hbar^2$, while for many-body problems Z_i is given by a more complicated integral.

Beyond a certain radius, the coupling interactions V_ii' in addition to the Coulomb potential are assumed to be extrapolated by means of a multipole expansion

$$V_{ii'}(r) = \frac{\hbar^2}{2m} \sum_{\lambda=1}^\Lambda F\sp {\lambda }_{ii'} r\sp {-\lambda -1}.$$ (22)

2.5.1 R-matrix Propagation

R-matrix propagator techniques have been widely used to integrate the equations describing atomic and molecular scattering problems. In essence the radial interval over which the equations are to be integrated is divided into sectors. The known R-matrix on one boundary of a sector is then used to obtain the R-matrix on the other. The process is continued so that the R-matrix is moved from one extreme of the interval of interest to the other. The principal advantage of this approach over alternative methods is very high numerical stability. A second important feature is that propagator methods are very economic in cases where large numbers of scattering energies must be considered. Much of the calculation which must be carried out within a given sector may be saved for energies subsequent to the first by saving appropriate sector information.

Equations (21) may be rewritten for a given solution, j, as

$$\frac{d^2\psi_{ij}(r)}{dr^2} = \sum_{i'=1}^{N} \tilde{V}_{ii'}(r)\psi_{i'j}(r)$$ (23)

by redefining the interaction potential as $\tilde{V}$. Assume that the radial interval between the R-matrix internal region, r=a and some asymptotic matching point, r=b, is divided into sectors. The R-matrix at the left-hand boundary of the p-th sector is defined by the equation

$$\psi_{ij}^p(r_L^p) = \sum_i' {\cal R}_{ii'}^p r \left . \frac{d\psi_{i'j}^p(r)}{dr} \right \vert _{r=r_L^p}$$ (24)

with an analogous result for the right-hand boundary. If the Green's function corresponding to equation (23) defined by equation

$$\frac{d^2{\cal G}_{ij}(r,r^\prime)}{dr^2}=\delta_{ij}\delta(r-r^\prime)+ \sum_{k} \tilde{V}_{ik}(r){\cal G}_{kj}(r,r^\prime)$$ (25)

is known within the sector then the R-matrices on the right- and left-hand boundaries are related by the equation

$$-{\cal R}_{ij}^p = {\cal G}_{ij}^p(r_R^p,r_R^p)+ \sum_{kl}{\... ...p(r_L^p,r_L^p) \right ]_{kl}^{-1} {\cal G}_{lj}^p(r_L^p,r_L^p)$$ (26)

It is clear from this equation that given suitable sector Green's function matrices evaluated at the sector boundaries the R-matrix propagation involves simple linear algebra which may be computed very rapidly using optimized blas routines.

Various propagator methods exist which differ on how the the sector Green's functions are determined. As the present package is designed to treat scattering in a region outside the range where the interaction potential is slowly varying, a potential following approach should be in general the most appropriate. The primary propagator method used is therefore one proposed by Light and Walker [12], [13], [14] which assumes that the interaction potential is constant within each sector. With this choice the sector Green's function matrices are rapid to compute and accuracy is maintained by varying the sector size. In practice it is found that sector sizes may become extremely large so this choice provides enormous computational advantages over solution-following methods.

2.5.2 Light-Walker Propagator (LWP)

In the Light-Walker (LW) approach the interaction potential is diagonalized at the centre of each sector. This defines an orthogonal transformation which may be used to decouple the Green's function equations within the sector. Therefore for the p-th sector centred at r=r_C^p if

$$\left ( {\bf T}^p \right )^T {\bf V}(r_C^p){\bf T}^p = ({\bf\lambda}^p)^2,$$ (27)

each of the uncoupled Green's functions, $\bar{\cal G}_i^p$, then satisfy the equations

$$\bar{\cal G}_i^p(r,r^\prime) = \delta(r-r^\prime)+(\lambda_i^p)^2 \bar{\cal G}_i^p(r,r^\prime)$$ (28)

throughout the sector. The uncoupled Green's function matrices required in equation (26) then may be written in the form

$$\bar{\cal G}_i^p(r,r^\prime) = \left \{ \begin{array}{lr} \f... ...R^p)\phi_i^p(r^\prime-r_R^p),&r>r^\prime, \end{array} \right .$$ (29)

where

$$\phi_i^p(r-r^\prime)=\left \{ \begin{array}{lr} \cosh\lambd... ...e),&(\lambda_i^p)^2 = -(\gamma_i^p)^2 <0. \end{array} \right .$$ (30)

and W represents the Wronskian of the functions $\phi_i^p$.

The sector size is determined by estimating the error resulting from variations in the potential over a sector. For example, the size of the (p+1)-th sector, L^p+1, is given by

$$\beta\left [ \left \{ 2\sum_{i=1}^N\frac{(\lambda_i^p)^2-(\l... ...p-1})^2} {N(L^p+L^{p-1})} \right \}^2 \right ]^{-\frac{1}{6}}.$$ (31)

where $\beta$ is a constant.

As discussed in Stechel et al [13] channels which do not contribute to the scattering can be dropped during propagation by simply truncating the R-matrix between one sector and the next. Such channels are deeply bound with slight coupling with open and less deeply bound channels. An estimate of the coupling of a channel j with channels (1,2,...,j-1) is given by the quantity q^p_j where

$$q^p_j = \frac {1}{L^p(2j-2)} \left [\sum_{k=1}^{j-1} \left \{(Q^p_{jk})^2 + (Q^p_{kj})^2 \right \} \right ]^{1/2}$$ (32)

where

$${\bf Q}^p = ({\bf T}^{p-1})^T {\bf T}^p$$ (33)

In propagating across a given sector we retain (1) all open channels, (2) a minimum number of closed channels (minnop in the program), (3) those less deeply bound than some criterion (deplim in the program) and (4) those deeply bound which have a value of q_j^p greater than some criterion (cuplim in the program). In addition the number of retained channels may decrease by only one from one sector to the next.

2.5.3 Asymptotic Expansions

The radial scattering equations must be integrated outwards from the boundary of the R-matrix internal region to a point where they can be matched to the asymptotic form of the solution with the appropriate boundary conditions. This integration region can be greatly reduced (and sometimes eliminated) by obtaining and matching to the solution of the scattering equations obtained analytically in the form of an asymptotic series which can be summed by the application of series acceleration techniques. In the present program an asymptotic series suggested by Gailitis[15] is summed using a simplified version of the Padé algorithm used by Noble and Nesbet [16]. The expansion selected has a number of advantages. First, a large part of any residual Coulomb interaction is removed from the asymptotic expansion and the recursion relation used to calculate the expansion coefficients contains terms which depend on the factor Z/k rather than on (Z/k)². Second, the expansion may be determined directly from a single set of recursion relations rather than indirectly by expanding and dividing out the Coulomb factors. Finally, the expansion is conveniently expressed in a form which is highly vectorizable provided the number of scattering channels is large.

Assume that a solution of each component of equations (21) may be written in the form

$$\psi_{ij}(r) = R(r) \sum_{\mu = 0} a_{ij}^{\mu} r^{-\mu} + \dot{R}(r) \sum_{\mu = 0} b_{ij}^{\mu} r^{-\mu}$$ (34)

where R(r) denotes a solution of the Coulomb equation

$$\left ( \frac{d^2}{dr^2} - \frac{L(L+1)}{r^2} + \frac{2Z}{r} + k^2 \right ) R(r) = 0$$ (35)

and

$$\dot{R}(r) = \frac{1}{k} \frac{dR(r)}{dr}.$$ (36)

Substituting the expansion given by equation (34) into equation (21) yields recursion relations for the expansion coefficients $a_{ij}^{\omega}$ and $b_{ij}^{\omega}$ of the form (assuming Z = Z_i)

$$(k_i^2 - k^2)a_{ij}^{\omega} + \left [ (\omega-2)(\omega-1)+L(L+1... ...{\cal L}_i+1)\right ] a_{ij}^{\omega - 2} + 2k(\omega - 1)b_{ij}^{\omega - 1} -$$

$$\sum_{j=1}^N \sum_{\lambda = 1}^{\Lambda} F_{ii'}^{\lambda} a_{i'... ...ga - 3)b_{ij}^{\omega - 2} - \frac{2}{k}L(L+1)(\omega -2) b_{ij}^{\omega-3} = 0$$ (37)

and

$$(k_i^2 - k^2)b_{ij}^{\omega} + \left [ (\omega - 2)(\omega - 1)+L(L+1)-{\cal L}_i({\cal L}_i+1)\right ] b_{ij}^{\omega - 2} -$$

$$\sum_{{i'}=1}^N \sum_{\lambda = 1}^\Lambda F_{i{i'}}^{\lambda} b_{i'j}^{\omega - \lambda - 1} - 2k(\omega - 1)a_{ij}^{\omega - 1}= 0$$ (38)

The required 2N open-channel solutions are found by choosing k² = k_i² for $i=1,\ldots,N$ and taking the regular Coulomb solution $R(r) \equiv F_{\cal L}(kr)$ for solutions corresponding to the boundary condition

$$\psi_{ij}(r) \mbox{\raisebox{-2.mm}{$\stackrel{\textstyle\si... ...tstyle r \rightarrow \infty }$}} \delta_{ij} \sin \theta_i(r)$$ (39)

and the irregular solution $R(r) \equiv G_{\cal L}(kr)$ for the boundary condition

$$\psi_{ij}(r)\mbox{\raisebox{-2.mm}{$\stackrel{\textstyle\sim... ...tstyle r \rightarrow \infty }$}} \delta_{ij} \cos \theta_i(r)$$ (40)

where

$$\theta_i(r) = k_ir - \eta_i \ln 2k_ir - \frac{1}{2}{\cal L}_i\pi - \sigma_{\eta_i}({\cal L}_i).$$ (41)

The Coulomb parameter, $\eta_i$, is given by

$$\eta_i = - \frac{Z_i}{k_i}$$ (42)

and $\sigma_{\eta}({\cal L}) \equiv \arg\Gamma({\cal L}+i\eta+1)$ is the Coulomb phase shift.

The coefficients a_ij⁰, b_ij⁰ are arbitrary if k_i² = k² or zero if $k_i^2 \neq k^2$. We therefore assume that $a_{ij}^0 = \delta_{ij}$ and b_ij⁰ = 0 so that

$$\psi_{ij}(r)\mbox{\raisebox{-2.mm}{$\stackrel{\textstyle\sim}{\scriptstyle r \rightarrow \infty }$}} \delta_{ij} R(r)$$ (43)

and take $L = {\cal L}_i$ with Z = Z_i.

In the case of closed channels, k_i² < 0, the solution R is taken equal to R(r) = H⁺(k_ir). Also $k_i \equiv i\vert k_i\vert$ and $\eta_i = i\vert\eta_i\vert$. With this choice the closed channel solutions asymptotically assume the usual form of Whittaker functions since

$$H^+_L(k_ir) \mbox{\raisebox{-2.mm}{$\stackrel{\textstyle\sim... ...right ]} W_{\vert\eta_i\vert,L+\frac{1}{2}}(2\vert k_i\vert r)$$ (44)

If we replace $b_m^\omega$ by $ib_m^\omega$ in the previous equations, the recursion relations for both open and closed channels may be expressed as

$$(k_i^2 \mp \vert k_j\vert^2)a_{ij}^\omega = -2\vert k_j\vert (\omega-1)b_{ij}^{\omega-1} \pm 2\vert\eta_j\vert(2\omega-3)b_{ij}^{\omega-2}$$

$$- \left [(\omega-2)(\omega-1)+{\cal L}_j({\cal L}_j+1)-{\cal L}_i({\cal L}_i+1) \right ] a_{ij}^{\omega-2}$$

$$\mp \frac{2}{\vert k_j\vert}{\cal L}_j({\cal L}_j+1)(\omega-2)b_{... ...um_{\lambda=1}^{\Lambda}F_{ii^\prime}^\lambda a_{i^\prime j}^{\omega-\lambda-1}$$ (45)

and

$$(k_i^2 \mp \vert k_j\vert^2)b_{ij}^\omega = - \left [(\omega-2)(\... ...1)+{\cal L}_j({\cal L}_j+1)-{\cal L}_i({\cal L}_i+1) \right ] b_{ij}^{\omega-2}$$

$$\mp 2\vert k_j\vert (\omega-1)a_{ij}^{\omega-1} + \sum_{i^\prime=... ...um_{\lambda=1}^{\Lambda}F_{ii^\prime}^\lambda b_{i^\prime j}^{\omega-\lambda-1}$$ (46)

for non-degenerate channels where $k_i^2 \neq k_j^2$. For degenerate channels, k_i²=k_j²,

$$2\vert k_j\vert \omega a_{ij}^{\omega} = \pm \left [\omega(\omega... ...\sum_{\lambda=1}^{\Lambda}F_{ii^\prime}^\lambda b_{i^\prime j}^{\omega-\lambda}$$ (47)

and

$$2\vert k_j\vert \omega b_{ij}^{\omega} = - \left [\omega(\omega-1)+{\cal L}_j({\cal L}_j+1)-{\cal L}_i({\cal L}_i+1) \right ] a_{ij}^{\omega-1}$$

$$\pm \frac{2}{\vert k_j\vert}{\cal L}_j({\cal L}_j+1)(\omega-1)b_{... ...\sum_{\lambda=1}^{\Lambda}F_{ii^\prime}^\lambda a_{i^\prime j}^{\omega-\lambda}$$ (48)

Given the expansion coefficients obtained from eqs (45-48) the solutions are defined by Eq. (34). The traditional methods for evaluating asymptotic series of this form are able to provide accurate solutions only at large radial distances r. Noble and Nesbet [16] have shown that series acceleration techniques may be used to evaluate the wavefunctions and their derivatives at much smaller radial distances. In the present package a simplified version of the algorithm used by Noble and Nesbet is employed.

An important feature of the approach we have adopted is that, given the coefficients of the asymptotic expansion, the wavefunctions and their derivative may be determined rapidly at a series of radial points in order to determine the minimum distance necessary to achieve the required level of accuracy. The accuracy of the results obtained is determined from the computed value of the multichannel Wronskian relations [17].

2.5.4 Scattering matrices

If all N scattering channels are open the asymptotic form of the radial wavefunction describing the scattering may be written in the form

$$\Psi(r) \mbox{\raisebox{-2.mm}{$\stackrel{\textstyle\sim}{\scriptstyle r \rightarrow \infty }$}} {\bf F} + {\bf G} {\bf K}$$ (49)

This asymptotic expression defines the scattering K-matrix, ${\bf K}$, appropriate for standing wave boundary conditions.

It is necessary to match the R-matrix obtained by outwards propagation from the internal region boundary to the asymptotic radial solutions obtained using the techniques described in the previous subsection. Given the asymptotic form of the radial wavefunctions defined by Eq. (49) it follows that the real symmetric K-matrix, K, is given by

$${\bf K} = - [ {\bf G} - r{\bf R}({\bf G}'-\beta{\bf G}) ]^{-1} [ {\bf F} - r{\bf R}({\bf F}'-\beta{\bf F}) ]$$ (50)

where ${\bf F}$, ${\bf G}$ are the regular and irregular full matrix solutions of Eq. (34) corresponding to Eqs. (39) and (40) respectively, evaluated at some matching point r. The matrix ${\bf R}$ represents the R-matrix evaluated at the same point. The K-matrix is related to the scattering S-matrix by the equation

$${\bf S} = \frac{{\bf 1} + i{\bf K}}{{\bf 1} - i{\bf K}} \ .$$ (51)

This equation is equivalent to Eq. (10).