8 Coupled Channels Methods

8.1 Coupled Reaction Channels

The numerical solution of Eqs. (16) is straightforward if the inter-channel couplings V_ij are local, as is the case for inelastic excitations of one or both nuclei. These are called `coupled channels' (CC) cases. Transfer couplings, however, couple different ${\bf R}_i$ and ${\bf R}_j$ values, giving what are called `coupled reaction channels' (CRC). The non-locality from a `finite range' treatment of recoil means that the coupled reaction equations must be solved either iteratively, or by a R-matrix treatment using square-integrable basis functions in an interior region. Local and iterative solution methods are presented in ref.[3], while the R-matrix methods presented in section 8.4, are common in atomic and molecular scattering research, but not so widely used in nuclear scattering problems. A variety of standard computer programs are available for evaluation of the couplings described above, and for solution of the coupled equations by the methods described below. The program PTOLEMY [17] can find coupled-channels solutions for local couplings or one-step non-local couplings from transfers, and ECIS [15] also solves coupled-channels equations. Both pay particular attention to the long-range couplings of Eq. (47) that arise from inelastic Coulomb excitations. The program FRESCO [3] includes these capabilities, as well as the iterative solutions of coupled equations with non-local couplings by all the methods to be now described. With non-local couplings from transfer channels, the Eqs. (16) may be solved iteratively, and the successive iterations amount to n-th DWBA solutions. As explained in section 3, the coupling matrix element $V_{ij} = \langle \phi_{pi}\phi_{ti} \vert H_m - E \vert\phi_{pj}\phi_{tj}\rangle$ has two different forms, depending on whether we use H_m=H_i (post form) or H_m=H_j (prior form). If we abbreviate $\Phi_i \equiv \phi_{pi}\phi_{ti}$, these give rise to the respective matrix elements

$$\left\langle\Phi_i \vert H_m - E \vert\Phi_j \right\rangle V_{ij}^{\rm post} + [T_i+U_i - E_i ] ~ K_{ij}$$ =

$$\mbox{ and } V_{ij}^{\rm prior} + K_{ij} ~ [T_j+U_j - E_j ]$$ = (72)

$$\begin{eqnarray*} \mbox{where } V_{ij}^{\rm post} &=& \langle \Phi_i \vert V_i ... ...\ \mbox{and } K_{ij} &= &\langle \Phi_i \vert \Phi_j\rangle . \end{eqnarray*}$$

The wave function overlap operator K_ij in equation (72) arises from the non-orthogonality between the transfer basis states defined around different centres in different mass partitions. We will see below that this term disappears in first-order DWBA, and can be made to disappear in second-order DWBA, if the first and second steps use the prior and post interactions respectively.

8.2 Multistep Born Approximations

If the coupling interactions V_i in Eq. (72) are weak, or if the back coupling effects of these interactions are already included in the optical potentials of the prior channel, then it becomes reasonable to use a n-step distorted wave Born approximation (DWBA). This approximation always feeds flux `forwards' in the sequence $1 \rightarrow 2 \rightarrow \cdots \rightarrow n+1$ neglecting the back couplings. In the elastic channel the wave function is governed by the optical potential defined there, and the wave function in the i'th channel is governed by the equation

$$\left [ E_i - T_i-U_i \right ]\psi_i ( {\bf R}_i ) = \sum _ {j=1}... ...ft\langle\Phi_i \vert {\cal H}-E \vert\Phi_j \right\rangle \psi_j ( {\bf R}_j )$$ (73)

Initial channel:

$$\left [ E_1 - T_1-U_1 \right ]\psi_1 ( {\bf R}_1 ) = 0$$

Second channel:

$$\left [ E_2 - T_2-U_2 \right ]\psi_2 ( {\bf R}_2 ) = \left\langle\Phi_2 \vert H_m-E \vert\Phi_1 \right\rangle \psi_1 ( {\bf R}_1 )$$ (74)

If the prior interaction is used, the right hand side can be simplified to

$$\begin{eqnarray*} &=&\langle\Phi_2 \vert {\cal V}_1 \vert\Phi_1\rangle\psi_1 ... ... \psi_1 \mbox{ is on-shell.}\\ &=& V_{21}^{\rm prior} \psi_1 \end{eqnarray*}$$

Final channel: c=n+1

$$\left [ E_c - T_c-U_c \right ]\psi_c ( {\bf R}_c ) = \sum _ {j=1}... ...} \left\langle\Phi_c \vert H_m-E \vert\Phi_j \right\rangle \psi_j ( {\bf R}_j )$$ (75)

If the post interaction had been used for all the couplings to this last channel, then there is again a simplification:

$$\left [ E_c - T_c-U_c \right ]\psi_c ( {\bf R}_c )$$

$$\sum_ {j=1} ^ {j=c-1} \langle\Phi_c \vert {\cal V}_c \vert\Phi_j\... ... [T_c+U_c - E_c ]~ \sum_{j=1} ^ {j=c-1} \langle\Phi_c \vert\Phi_j\rangle \psi_j$$ =

so

$$\left [ E_c - T_c-U_c \right ] \chi_c ( {\bf R}_c ) = \sum_{j=1} ^ {j=c-1} V_{nj}^{\rm post} \psi_j$$ (76)

where

$$\chi_c ( {\bf R}_c ) =\psi_c + \sum_{j=1} ^ {j=c-1} \langle\Phi_c \vert\Phi_j\rangle\psi_j =\langle\Phi_c \vert\Psi\rangle \nonumber$$

Note that, as all the $\Phi_i \equiv \phi_{pi}\phi_{ti}$ are square-integrable and hence decay faster than r -1 at large radii, the $\psi_c$ and $\chi_c$ are the same asymptotically. They differ only by an `off-shell transformation', and hence yield the same (on-shell) scattering amplitudes. The equation for $\chi_c$ has no non-orthogonality terms once the post interaction is used in the final channel: this is what is meant by saying that the final channel is `effectively on-shell'. These results imply that in n-step DWBA, some non-orthogonality terms can be made to disappear if `prior' interactions are used for the first step, and/or if `post' interactions are used for the final step. This means that the non-orthogonality term never appears in the first-order DWBA, irrespective of the choice of prior or post forms. In second-order DWBA, the prior-post combination must be chosen[14] to avoid the non-orthogonality terms. It is clear that non-orthogonality terms will have to be evaluated if the DWBA is continued beyond second order.

8.3 Iterative Solutions

The iterative method of solving the CRC equations (16) proceeds by analogy with the n-step DWBA iterations until the series converges. Convergence is readily obtained if the couplings are sufficiently small, and different iterative strategies may be employed [17]. The procedure will however diverge if the the couplings are too large, or if the system is too near a resonance or a bound state pole. On divergence, the successive wave functions $\psi_i^{(n)}$ will become larger and larger as n increases, and not converge to any fixed limit. Unitarity will of course be violated as the S-matrix elements will become much larger than unity. In this case we may use Padé approximants to accelerate the convergence of the sequence $S^{(n)} _\alpha$ of S-matrix elements[15,17]. A given sequence $S^{(0)} , S^{(1)} , \cdots$ of S-matrix elements that result from iterating the coupled equations can be regarded as the successive partial sums of a `vector valued' polynomial

$$f(\lambda) = S^{(0)} + (S^{(1)} - S^{(0)}) \lambda + (S^{(2)} - S^{(1)}) \lambda ^ 2 + \cdots$$ (77)

evaluated at $\lambda$=1. This polynomial will clearly converge for $\lambda$ sufficiently small, but will necessarily diverge if the analytic continuation of the $f (\lambda)$ function has any pole or singularities inside the circle $\vert \lambda \vert<1$ in the complex $\lambda$-plane. The problem that Padé approximants solve is that of finding a computable approximation to the analytic continuation of the $f (\lambda)$ function to $\lambda$=1. This is accomplished by finding a rational approximation

$$P_{[N,M]} (\lambda) = {{p_0 + p_1 \lambda + p_2 \lambda^2 + \cdot... ... \lambda^N } \over {1 + q_1 \lambda + q_2 \lambda^2 + \cdots + q_M \lambda^M }}$$ (78)

which agrees with the $f (\lambda)$ function in the region where the latter does converge, as tested by matching the coefficients in the polynomial expansion of $P_{[N,M]} (\lambda)$ up to and including the coefficient of $\lambda^{n}$ for n=N+M. There are many different ways[18] of evaluating the coefficients p_i , q_j, but for the present problem we can use Wynn's $\epsilon$-algorithm[19,20], which is a method of evaluating the upper right half of the Padé table at $\lambda$=1 directly in terms of the original sequence $S^{(0)} , S^{(1)} , \cdots$. Experience has shown that for typical sequences the most accurate Padé approximants are those near the diagonal of the Padé table. We use $\overline{S^{(n)}} \equiv P_{[N,M]} (1)$ for N=[(n+1)/2] and M=[n/2] in calculating the Padé-resummed cross sections.

8.4 R-matrix Solutions

The radial stepping methods of solving the coupled equations only allow local couplings to be treated properly, and non-local couplings from transfers have to be included iteratively. The R-matrix method[21] is an equivalent way of solving the coupled equations, and has the advantages of being more stable numerically, and also allowing non-local components of the Hamiltonian in an interior region to be included to all orders. It has recently been revived in nuclear physics applications[22,23] for these reasons. Both transfer and non-orthogonality non-localities may be included non-perturbatively, and resonances and bound states may be described without difficulty. This method uses a basis set of `energy eigenstates' of the diagonal parts of Eqs. (34):

$$\Bigl( T_{iL}(R_i) + U_i(R_i)+ e_{pi}+e_{ti}- \varepsilon_q \Bigr) w^q_{i}(R_i) = 0$$ (79)

for eigenenergies $\varepsilon_q$, with the basis functions all having fixed logarithmic derivatives $\beta = d\ln w^q_{i}(R_i)/dR_i$ at R_m. The constancy of the logarithmic derivatives $\beta$ means that (for each i channel separately) the w^q_i form an orthogonal basis set over the interval [0,R_m], and over this range they can be normalised to unity. Without this constancy, a Bloch operator must be added to the kinetic energy to make it Hermitian. The wave functions of the coupled problem (34) can now be solved completely over the interior range [0,R_m], by using the orthonormal basis set of the {w^q_i(R_i)} with coefficients to be determined. The coefficients are found in two stages: first by finding all the eigensolutions g^p_i(R_i) of the equations (34) using the above orthonormal basis, and then expanding the scattering wave functions in terms of these g^p_i(R_i). In the traditional R-matrix method, the diagonalisation of the N-channel Hamiltonian in equation (34) yields P=QN eigenenergies e_p with corresponding multichannel eigenstates

$$g^p_{i}(R_i) = \sum_{q=1}^Q ~ c^{pq}_{i} w^q_{i}(R_i)$$ (80)

Eigenstates here with e_p<0 are close to 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 g^p_i(R_i) 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 _{ii_0} = \sum_p A^p_{ii_0} g^p_{i}$. If we define an R-matrix at energy E by

$$\psi _{i}(R_i) = \sum_{i'} R_{ii'}(E) \Bigl [ \psi{'} _{i'}(R_i) - \beta \psi _{i'}(R_i) \Bigr]$$ (81)

in the limit of $R_i\rightarrow R_m$ from above, then the R-matrix $\bf R$ can be calculated from the eigenstates by standard methods[21,24]

$$R_{ii'}(E)= \frac{\hbar^2}{2\mu_i'} \sum_{p=1}^P \frac{ g^p_{i}(R_m) g^p_{i'}(R_m) } { e_p - E} \ .$$ (82)

The coefficients c^pq_i and energies e_p in Eq. (80) satisfy matrix equations

$$\varepsilon_q c^{pq}_{i} + \sum_{q'i'} \langle w^q_{i} \... ... i'} \mid w^{q'}_{i'} \rangle c^{pq'}_{i'} = e_p c^{pq}_{i}$$ (83)

for each eigenstate p, where V_ii' refers to all the off-diagonal couplings. These equations are of the matrix form

$${\bf H} {\bf c} = e {\bf c} \ .$$ (84)

There is an alternative method[25,23] for finding the R_ii', which does not diagonalise the matrix on the left side of Eq. (84), but solves a set of linear equations. We need the solution of $({\bf H}-E) {\bf x} = {\bf w}(R_m)$ for the right hand side consisting of the values of the basis functions at the R-matrix boundary. Then we can solve directly

$${\bf R} = \frac{\hbar^2}{2\mu_i} {\bf w}^T(R_m)({\bf H}-E)^{-1} {\bf w}(R_m) \ .$$ (85)

This has the advantage of naturally continuing the R-matrix method to complex potentials, avoiding the diagonalisation of non-Hermitian matrices. Using Eqs. (36) and (81), 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^+}\! -\! {\bf R}({\bf H'^+}\! -\! \beta ... ...[ {\bf H^-}\! -\! {\bf R}({\bf H'^-}\! -\! \beta {\bf H^-}) ]$$ (86)

and the expansion coefficients for the wave functions are

$$A_{ii_0}^p = - \frac{\hbar^2}{2\mu_i'} \frac{1}{e_p-E}$$

$$\sum_{i'} g^p_{i'}(R_m) \Bigl [ \delta_{ii_0} (H'^-_L(k_iR_m)-\beta H^-_L(k_iR_m))$$

$$- S_{ii_0} (H'^+_L(k_iR_m)-\beta H^+_L(k_iR_m)) \Bigr ] \ .$$ (87)

8.4.0.1 Buttle Correction

The R-matrix calculated by Eq. (82) 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[26] is added to the right hand side of Eqs. (82,85). This modifies the diagonal terms R_ii(E) to reproduce for each uncoupled problem the exact scattering solution $\chi_{i}(R_i)$ after this has been integrated separately. From the definition of the energy eigenstates w^q_i(R_i), the R-matrix sum from (82) for each uncoupled channel is

$$R^u_{i}(E)= \frac{\hbar^2}{2\mu_i} \sum_{q=1}^Q \frac{ w^q_{i}(R_m)^2 } { \varepsilon_q - E}$$ (88)

and the exact one-channel R-matrix is $R^0_{i}(E) = \chi_{i}(R_m)/(\chi'_{i}(R_m) - \beta \chi_{i}(R_m))$. The Buttle-corrected full R-matrix to be used in Eq. (86) is then

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

The energy $\tilde E$ can be equal to E, or chosen just near to it if necessary to avoid the poles in Eq. (88), since the Buttle correction varies smoothly with energy.

8.4.0.2 CRC matrix elements

The solution of the CRC equations (16) with all the non-orthogonality terms in Eq. (72) requires in Eq. (83) the matrix element integrals of the form

$$\langle w^q_{i}\vert V_{ii'}\vert w^{q'}_{i'} \rangle = \la... ...vert H_m - E \vert\Phi_{i'}\rangle \mid w^{q'}_{i'} \rangle$$ (90)

for m=i (post) or m=i' (prior). In the post form, H_m contains T_i+U_i, and since w^q_i is just the eigenfunction of this operator with eigenvalue $\varepsilon_q$, we can operate to the left to obtain

$$\langle w^q_{i}\vert V_{ii'}\vert w^{q'}_{i'} \rangle_{\rm ... ..._i) \langle w^q_{i} \Phi_i \vert\Phi_{i'} w^{q'}_{i'} \rangle$$ (91)

with the similar prior form

$$\langle w^q_{i}\vert V_{ii'}\vert w^{q'}_{i'} \rangle_{\rm ... ...'}) \langle w^q_{i} \Phi_i \vert\Phi_{i'} w^{q'}_{i'} \rangle$$ (92)

The wave function overlaps in the second term $\langle \Phi_i \vert\Phi_{i'} \rangle$ go to zero asymptotically, and may be assumed small when R_i, R_i'>R_m. The standard R-matrix theory therefore still applies in the asymptotic region.