5 Numerical Solutions

This section discusses the methods used to solve the coupled reaction channels equations (30), when in general there are both local couplings $V_{\alpha :\alpha'}^\Gamma (R_\kappa)$ and non-local kernels $V_{\alpha :\alpha'} (R_\kappa , R_{\kappa'} )$. Now a group of m equations can be solved `exactly' (subject only to radial discretisation errors) by finding [53] a set of m linearly independent groups of solutions g<sub>i,k</sub> (R), and taking a linear combination of these which satisfies the required boundary conditions. This method is only practicable, however, if there are not too many equations (the numerical effort can rise as m<sup>3</sup>), and if there are only local couplings. For in that case the independent solutions can be found in a single outward `sweep' of m² radial functions. Non-local couplings mean, unfortunately, that the source terms at a given radius depend on the wave functions at other radii both larger and smaller, so that this `exact' method becomes impractical.

In many cases of interest in nuclear physics, however, the non-local couplings are not too strong, and can be treated as successive perturbations. They can then be applied iteratively until further applications have progressively smaller efffects, and the solutions have converged (to some preset criterion of accuracy). Some failures of convergence can remedied by the use of Padé approximants.

When both local and non-local couplings are present, and the local couplings are too strong to allow an iterative scheme to converge, a combination of the exact and iterative schemes is possible. In this approach, several channels can be `blocked' together, and treated as one unit during the iterations, while solving the couplings within the block by the exact method.

There are several other features of typical nuclear reaction calculations which simplify the numerical methods:

1. If the sum of the incoming projectile and target spins is greater than one, then solutions will often be required for the same set of CRC equations, only with different boundary conditions.
2. The diagonal potentials $U_\kappa (R_\kappa)$ usually have a significant imaginary component for small $R_\kappa$, and hence damp the solutions $f_\alpha (R_\kappa )$ in this region. This enables lower radial cutoffs to be used for $R_\kappa$ near zero, with little loss of accuracy.
3. The bound states $u_{\ell s j I} (r)$ used in transfer reactions decay exponentially outside the surface region of the nuclei. This means that the integrand in equation (113) for the transfer kernels will often decay exponentially both as | R - R' | increases, and as $u \equiv cos \theta \equiv \hat {{\bf R}} \cdot \hat {{\bf R}'}$ decreases from unity.

5.1 Integration of the Radial Equations

If the non-local interactions $V_{\alpha ,\alpha'} (R,R' )$ in the CRC equations (30) are present, then it will always be necessary to solve the coupled channels by iteration. With the local couplings $V_{\alpha ,\alpha'}^\Gamma (R),$ however, we have a choice whether to iterate, or to include them in the exact solutions of the close-coupling method. A simple option is to allow a specifiable number b of channels to be coupled exactly, with the remainder only being fed after one or more iterations. This would be useful, for example, if the channels for the low-lying states of a highly-deformed target were included in this block of b channels, and if the remaining channels (e.g. for transfers) were not fed by more than 2 or 3 steps beyond this initial block. Restricting these iterations to one is equivalent to solving a CCBA model.

Whether the coupled equations are of the simpler form of equation (30), or of the more complex form of section 2.3, a particular n'th iteration will require solving set of m equations of the form

$${d^2 \over dR^2} f_i (R) \sum _ {j=1} ^ b C_{ij} (R) f_j (R) + S_i (R) \mbox{ for } i=1 \cdots b ,$$ = (123)

$$\mbox{ and } {d^2 \over dR^2} f_i (R) C_{ii} (R) f_i (R) + S_i (R) \mbox{ for } i=b+1 \cdots m ,$$ = (124)

where S_i (R) is the source term constructed by means of the wave functions f_i^(n-1) (R) of previous iterations :

$$S_i (R) = \sum_{j = j_{\rm min}}^ m C_{ij} (R) f_j^{(n-1)} (R)$$ (125)

where $j_{\rm min} = b+1 \mbox{ if } i \leq b$ and $j_{\rm min} = 1$ if i> b.

These coupled differential equations can be solved, following the method of ref. [53] by forming the linearly independent solution sets g_i,k (R), where the k'th solution consists of a set of all channels ($i=1 \cdots m$) which is made independent of the other sets by having a distinctive starting value

$$g_{i,k} (R_{\rm min} - h) = 0, ~~~~~ g_{i,k} (R_{\rm min}) = {1 \over (2L_i + 1)!!} ~ (K_i R_{\rm min})^{L_i + 1} \delta_{i,k}$$ (126)

for the initial conditions in the radial integration of equations (123). For this integration, the code FRESCO uses the modified Numerov method, and other codes such as Tamura's JUPITOR [26] have used Euler's method to start with near the origin (R=0), and then Störmer's 6-point method to continue. A general discussion of numerical integration schemes is given in Melkanoff et al. [27], along with error analyses of the different methods.

The independent solutions g_i,k (R) are required for m+1 values of k. The solution vectors for $k=1 \cdots m$ are solved starting with equation (126) but with no source term in the equation (123): these will contribute to the complementary solution of the homogeneous equation. We also need a particular solution g_i,0 (R) of the inhomogeneous equation, solved with the source terms but with no non-zero values in equation (126). These partial solutions may be conveniently laid out as in figure 5.1. If, however, it is known that the wave functions of certain channels are not required (if, for example, they are only fed in the last iteration), then it is not necessary to store their components in the array, for their S-matrix elements can still be calculated.

The solutions f_i (R) are the linear combination of the g_i,k (R)

$$f_i (R) = \sum _{k=0} ^ m a_k g_{i,k} (R)$$ (127)

satisfying the boundary conditions of equation (32) at R = R_m and say R = R_m - 5h. The coefficient a₀ is always unity, to match the source terms correctly. The S-matrix elements are a by-product of the linear matching equations (32).

Figure 2: Independent Solution Vectors: Layout of the independent radial wave functions for solving a system of m equations, of which b are to coupled exactly. Each entry represents a vector of n radial points, and the entries in bold are those with a non-zero initial values for their outward radial integration.

Note that the independent solutions g_i,k (R) for $k \geq 1$ need only be calculated the first time this coupled channels set is used. If they are stored as in figure 5.1, subsequent iterations need only recalculate the first row (IT=1) as the source terms vary. Furthermore, if there are multiple incoming channels for fixed total spin J_T and parity $\omega$, then solutions after the first can also use the g_i,k (R) already stored. The first iteration for these subsequent incoming channels will in fact not require any radial integrations whatsoever, merely finding a new set of { a_k } from the new matching conditions, and recalculating the sum (127) if the wave functions are required.

Tolsma and Veltkamp [54] point out one difficulty with this method, which is that if the couplings C_i,j are strong for $i \neq j$, then the linear independence of the g_i,j (R) will be reduced as R increases through a classically forbidden region. This is because the components with negative local kinetic energy will generally consist of an exponentially growing part and an exponentially decreasing part. The former is responsible for the tendency to destroy the initially generated linear independence of the solution vectors. The longer the integration continues through a classically forbidden region, the stronger this tendency will be; for instance, it will occur in scattering problems of nuclear physics with energies near or below the Coulomb barrier. It will also occur if inelastic form factors are used which are not approximately derivatives of the diagonal potential, but which extend more than usual into the interior of the nucleus that is classically forbidden because of the centrifugal potentials.

Tolsma et al. [54] propose a stabilization procedure to monitor and if necessary re-orthogonalise the solution vectors. If this were not done, there would be large cancellations in the sum of equation (127), resulting if severe in complete loss of accuracy of the S-matrix elements and the solution wave functions.

A simpler approach is to increase the starting radius $R_{\rm min}$ at which the radial integrations begin. It is advisable in any case for reasons of stability at small radii to have a minimum radius proportional to some angular momentum $\overline L$ typical of the coupled channels set:

$$R_{\rm min} \geq c \overline L h$$ (128)

for some constant c in the region of 1 or 2, where h is the radial step size. This constant could be increased to avoid the loss of independence in the present problem, but this is not always satisfactory, as the absorptive effects of the optical potentials at intermediate radii might thereby be lost. An alternative remedy (adopted in ref.[34]) is to have a specifiable radial cutoff R⁽c)_min for the off-diagonal coupling terms only. This allows the absorption in the diagonal potentials to be effective at all radii outside $R_{\rm min}$ of equation (128), but does not allow any strong coupling terms to lead to loss of independence until some larger radius which can be adjusted to keep the loss of accuracy to an acceptable level. It thus should be a regular policy in a computer code to integrate the equations (123) to a precision of at least 12 to 16 significant figures, to monitor the degree of cancellation in equation (127), and to notify the user should this approach within 2 or 3 powers of ten of the precision limit of the computer. Note that it is not necessary for the coupling terms C_ij (R) (etc) to be accurate to full machine precision, only that they should be consistently precise when converted to that precision.

5.2 Convergence of the Iterative Method

The iterative method of solving the CRC equations (5, 30) will converge if the couplings are sufficiently small. The procedure will however diverge if the the couplings are too large, or if the system is too near a resonance. 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.

5.2.1 Improving the Convergence Rate

There are several ways of dealing with this problem:

1. Solving some of the local couplings exactly by the methods of section 5.1, and iterating only on the non-local couplings and the remaining local couplings.
2. Solving all the channels simultaneously via a very lar ge system of linear equations, with each radial point in each channel as a separate unknown [28].
3. Find a separable expansion for the non-local kernels, so that they can be included exactly in the coupled-channels solution [49].
4. Expand the wave functions with a range of basis states of Coulomb and (say) Gaussian [50] or Airy [5] functions, and take the coefficients in this basis as the unknowns in a system of linear equations.
5. Use Padé approximants to accelerate the convergence of the sequence $S^{(n)} _\alpha$ of S-matrix elements [51,52].
6. Iterating the equations sequentially as in [51] and [52], rather than all equations as a block.
7. The inwards-outwards method of refs. [29], [30] and [37].

For the range of heavy and light-ion reactions that we are considering here, the methods (1) and (5) above have been adopted.

The method (2) is not used because of the size of the matrix that results. Initially, the matrix would be sparse, with selected elements away from the diagonal being non-zero because of the coupling potentials. The kinetic energy operators occupy a band of width three along the diagonal. Although a Gaussian elimination procedure would allow potentials of arbitrary coupling strength to be included, it will fill in large regions of the matrix as the solution proceeds, and this makes the storage requirements prohibitive.

The separable expansion method (3), while useful for light-ion reactions, is unsatisfactory for heavy-ion transfers. This is because if the masses of the initial and final nuclei become large relative to the mass of the transferred particle, the form factor for the transfers becomes more nearly local. As we approach the no-recoil limit (which makes the form factors exactly local) a separable expansion of a nearly-local kernel will require a large number of terms. In the limit of a local form factor, the separable expansion will require the same number of terms as there are radial points.

The method (4) of expanding the wave functions in Gaussians could have been used, provided the characteristic widths in R-space of the basis states were chosen in accordance with the wave number $K _\alpha$ in the relevant channel. This requirement is less severe with light-ion reactions, where the wave numbers are typically $\leq 1$ fm^-1. For heavy-ion reactions, however, the oscillation rates are much larger, and a more sensible method is to expand in terms of Airy functions that are depend explicitly on the local wave number over some radial region.

It is very useful to be able to iterate the coupled equations in a conventional manner, as then 1, 2 and 3 step DWBA results (etc.) can be recovered by stopping the iterations short of full convergence. This recovery of DWBA results is more difficult with sequential iteration (6), but both that method and the method of (7) would be definitely advantageous when, say, exciting a long rotational band by successive application of a quadrupole coupling. Using Padé acceleration has the advantages that it need only be employed if ordinary iterations are seen to diverge, and that it transforms the previously-divergent results with little new computational effort.

5.2.2 Padé Approximants for Sequence Extrapolation

A given sequence $S_0 , S_1 , \cdots$ of S-matrix elements that result from iterating the CRC equations can be regarded as the successive partial sums of the polynomial

$$f(\lambda) = S_0 + (S_1 - S_0) \lambda + (S_2 - S_1) \lambda ^ 2 + \cdots$$ (129)

evaluated at $\lambda$=1. This polynomial will clearly convergence 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. This is accomplished by finding a rational approximation

$$P_{[n,m]} (\lambda) = {p_0 + p_1 \lambda + p_2 \lambda^2 + \cdots + p_n \lambda^n \over 1 + q_1 \lambda + q_2 \lambda^2 + \cdots + q_m \lambda^m }$$ (130)

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+m}$.

There are many different ways [48] of evaluating the coefficients {p_m , q_n}, but for the present problem we can use Wynn's $\epsilon$-algorithm [31], 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$.

5.2.3 Wynn's epsilon Algorithm

Initialising $\epsilon_0^{(j)} = S_j$ and $\epsilon^{(j)}_{-1} = 0$, we form an array using the relation $\epsilon^(j)_{k+1} = \epsilon^{(j+1)}_{k-1} + ( \epsilon_k^{(j+1)} - \epsilon_k^{(j)} ) ^ {-1} .$ Thus we can construct the table given the second column from the initial sequence S_j. The table then gives the transposed upper right half of the Padé table, including the diagonal:

$$\epsilon^{(j)}_{2k} = P_{[k,k+j]} (1) .$$ (131)

Experience has shown that for typical sequences the most accurate Padé approximants are those near the diagonal of the Padé table, and these are just the right-most $\epsilon^{(0)}_{2k}$ in the $\epsilon$ table.

When accelerating a vector S-matrix elements ${\bf S}_j$, with a component for each coupled channel, then it is important to accelerate the vector as a whole. Wynn [32] pointed out that this can be done using the Samuelson inverse

$${\bf x}^{-1} = ({\bf x} \cdot {\bf x}^* ) ^ {-1} {\bf x}^*$$ (132)

where ${\bf x}^*$ is the complex conjugate of ${\bf x}$. Otherwise there will be problems when iterating (say) a two-channel system with alternating backwards and forwards coupling, because of zero divisors in the $\epsilon$ algorithm.

5.3 Transfer Form Factors

5.3.1 The Cancellation Problem

As discussed in section 4.4.1, the summations over T in equation (114) involve large cancellations, and as the degree of cancellation gets worse for large $\ell$ and $\ell'$, this places a limit on the maximum value $\ell + \ell'$ of the transferred angular momentum.

Typically, however, the transfer form factors are only needed to be accurate to around 0.1 to 1%, so if computer arithmetic is used which is accurate to 14 or 16 decimal digits, then cancellations up to 12 or 13 orders of magnitude should in principle not result in catastrophic errors in the transfer rates. With careful programming, this accuracy can be achieved. What is necessary is to be careful that all quantities in the equations (114, 113) above which depend on the Legendre order T are calculated to the full computer precision. It is not necessary, for example, for the channel wave functions $f _\alpha (R)$, the bound state wave functions $u_{\ell s j I} (r)$ or the quadrature of the integral (113) to be accurate to full precision (which in any case would be impossible). It is only necessary that all these quantities have exactly the same computer precision when the coefficients over T (the ${\bf q}_{\ell , \ell'}^T (R,R') )$ are evaluated, and when the sums over T (in equation 114) are performed. This will require principally that the `radial framework' that gives ${\bf r}$ and ${\bf r}'$ in terms of ${\bf R}$ and ${\bf R}'$ be accurate to full machine precision, as also the Racah algebra coefficients in equation (114). In fact, the channel wave functions $f _\alpha (R)$ and the bound state wave functions $u_{\ell s j I} (r)$ may be calculated with reduced precisions using shorter computer words and faster arithmetic should these be available. It is also not necessary for the coefficients and sums over T be consistent to full accuracy for different R and R' values, as the large cancellations only occur between different T values for each separate R and R' combination.

Since the accuracy of the quadrature in the equation (113) is not critical to the overall accuracy of the transfers, calculations may be speeded up if we economise on the range of the u variable and on the number of intermediate steps required. Even in light ion reactions it is not necessary to integrate u to -1 ($\theta$ to 180$^\circ$) as was done in the code LOLA [72] for example. An efficient procedure to adopt is that used in the DWBA code DAISY[55], where, for each successive R value, the code monitors the rate of decay of the integrand as $\theta$ increases. For a given accuracy criterion, an estimate can then be made of an adequate upper limit for the $\theta$ integration at the next R value. Typically, the upper limits of $\theta$ decrease monotonically as R increases from 0 to the upper limit R_m. Because the integrand is largest for $\theta$=0, the accuracy of the angular integration for small $\theta$ is improved by a change of variable from u to x as in ref.[55]:

$$\theta = \frac{1}{4} (3 x^2 + 1) x \theta_{\rm max}$$ (133)

for $0 \leq x \leq 1.$ The quadrature over u of equation (113) then becomes

$${\bf q}^T_ {\ell , \ell'} (R,R') = \frac{1}{2}\int_ 0 ^ 1 {\bf V}... ... sj'} (r')\over {r'}^{\ell' +1}} ~ P_T (u) \sin(\theta) {d \theta \over dx} dx.$$ (134)

The parameter $\theta_{\rm max}$ is adjusted for each successive value of R. according to the rate at which the integrand is observed to decay as $\theta$ increases, as described earlier.

5.3.2 Radial Grids

The methods used to calculate, store and use the non-local form factors $qu^T_{\ell , \ell'} (R,R')$ (equation 113) and $V_{\alpha ,\alpha'} (R,R' )$ (equation 108) have to be efficient in a wide variety of reactions, from light-ion reactions such as ³He(³H,⁴He)²H or ¹⁶O(²⁰Ne,²⁴Mg)¹²C to heavy-ion reactions, such as nickel on tin one-nucleon transfers. In the former cases, the radial form factors $V_{\alpha ,\alpha'} (R,R' )$ will be non-zero over large regions of the R-R' space, so (following ref [56]) interpolation procedures should prove effective.

However, when small masses are transferred between two larger nuclei the form factor is nearly local, and only large around $R \sim R' .$ If the whole (R,R') array had to be calculated and stored in these cases, modelling heavy-ion transfers would become inefficient, even with interpolation methods. The form factor now varies rapidly as a function of $\delta R \equiv R - R'$ (especially for heavy ion reactions, as the Jacobian b³ in equation (114) becomes large), and varies only slowly with R (if $\delta R$ is constant), as this variation follows the radial dependence of the bound state wave functions. The best procedure is thus [56] to first change to the coordinate pair $\delta R$ and R, and then to use different interpolatory intervals $h_\delta$ and h_R in the two directions respectively. Then, when nuclear masses become large compared with the mass of the transferred particle, $h_\delta$ can become smaller, perhaps even smaller than h, the basic radial step size.

The method adopted in FRESCO is to let the user specify $h_\delta$ and h_R as multiples or submultiples of h. The value of h_R is very often always 3 to 5 times larger than h, as this reflects the typical rate at which bound state wave functions vary. If the bound state wave functions have many internal nodes, then the interpolation interval h_R cannot be so large (this is often the case with $\alpha$-particle bound states).

The $h_\delta ,$ on the other hand, will be larger than h for light-ion reactions (as described in [47]), but will be comparable with or smaller than h for few-nucleon transfers between heavy ions. The user also specifies the maximum and minimum values of the range of $\delta R ,$ which again will be large ( $\sim$ nuclear radii) for light ions, and small ( $\sim$ 1 or 2 fm.) for heavy ion reactions. The accuracy of these choices is checked retrospectively by collecting statistics on the distributions of the functions ${\bf q}^T_{\ell , \ell'} (R, \delta R),$ averaging over R. and all partial waves T, $\ell$, and $\ell' .$

When $h_\delta$ or h_R are multiples h, then (say) cubic splines in two dimensions can be used to expand the form factors for the integrals of equation (106). If, however, $h_\delta$ is a submultiple h, as is the case in many heavy-ion reactions, then a more efficient procedure is possible.

Suppose, say, we wish to evaluate the numerical integral

$${\cal I} = \sum _ j V(x_j) \overline f (x_j) ,$$ (135)

where the $\overline f (x_j)$ are the interpolated values of the function f(x) between its stored values f_i at x = (i-1)h . Let the interpolation method be linear:

$$\overline f (x) = \sum _ m a_m (x) f_m$$ (136)

for some x-dependent coefficients a_m (x) from (say) fitting cubic splines over some range (most of the a_m will be zero except for $m\sim i\pm 2$). Then $\cal I$ can be evaluated directly in terms of the f_m :

$${\cal I} \sum _ j V(x_j) \sum _ m a_m (x_j) f_m$$ = (137)

$$\sum _ m \overline V_m f_m$$ = (138)

where

$$\overline V_m \equiv \sum _ j V(x_j) a_m (x_j)$$ (139)

is a new effective form factor This means that when $h_\delta$ is a submultiple of h, we do not need to store a form factor at intervals of $h_\delta ,$ only at intervals of h, if we use the `preemptive interpolation' of equation (139). This has the further advantage that as the no-recoil limit is approached (as the mass of the transferred particle becomes a smaller fraction of the interacting nuclei), then the form factors $\overline{\bf q}^T_{\ell , \ell'} (R, \delta R)$ and $\overline V_{\alpha ,\alpha'} (R, \delta R)$ need fewer grid points in the $\delta R$ direction. Less arithmetic is needed to evaluate the source functions of equation (106), which change from

$$S _\alpha (R) = \int_ {\delta R_{\rm min}} ^ {\delta R_{\rm max}}... ... V_{\alpha ,\alpha' } (R, \delta R) ~ f_{\alpha'} (R - \delta R) d( \delta R) .$$ (140)

to

$$S _\alpha (i h_R) = h \sum _ j \overline V_{\alpha ,\alpha' } (ih_R , j h) f_{\alpha'} ((in_R - j)h) \mbox{ where } n_R \equiv h_R / h$$ (141)

even when the original kernel functions vary rapidly as $\delta R$ changes in steps of h (with R constant).

Simultaneous Two-Nucleon Transfers: A similar `preemptive' summation is possible when calculating the form factors for the simultaneous transfer of two nucleons between states of the form of equation (50) in the projectile and in the target. As mentioned in section 3.4, two-nucleon transfer can be viewed as the transfer of a `structured particle' with internal coordinates $( \ell , (s_1 s_2) S ) j$ and $\rho$, the distance between the two nucleons. A transfer is only possible if the initial and final states have identical values for these `internal coordinates'. The angular momentum quantum numbers can be matched exactly, but source terms can either be constructed for each $\rho$ value and summed in equation (106), or the separate $\rho$ products can be summed as early as equation (113). Because the separate $\rho$ values are only used in a summation, it is most economical to use Gaussian quadrature, as for a given accuracy this reduces by a half the number of rho_i values at which the wave functions of equation (50) need to be calculated and stored. If the rho_i are chosen to be the Gaussian quadrature points over some chosen range, and if w_i are the corresponding weights, the equation (113) becomes

$${\bf q}^T_ {L , L'} (R,R') = \frac{1}{2}\int_ {-1}^{+1} r^{-L-1} ... ...[ \sum _ i w_i {\bf V} u_{12} (r, \rho_i) ~ u'_{12} (r' , \rho_i) ] P_T (u) du.$$ (142)

Equation (114) remains unchanged, and this means that two nucleon transfers can be calculated efficiently with little more computational work than that required for single-particle transfers.

Acknowledgements

I would like to thank Dr. Nagarajan, John Lilley, Ray Mackintosh, Victoria Andres and a referee for valuable questions and comments at various stages in this project. Manchester University, Daresbury Laboratory, the Niels Bohr Institute, and Bristol University are thanked for their hospitality when these ideas were being developed, and the Department of Engineering Mathematics at Bristol is thanked for their assistance in the preparation of this paper.