Coupled Reaction Channels Calculations in Nuclear Physics

Ian J. Thompson Department of Engineering Mathematics, Bristol University, Bristol, BS8 1TR, UK

Address after 1st August 1988: Department of Physics, Surrey University, Guildford, GU2 5XH, UK

Computer Physics Reports 7 (1988) 167-212 Latex version of Jan 2000. With corrections to original article. Reprinted May 7, 2004

Abstract

This paper describes the components and methods of a comprehensive code for coupled reaction channels calculations in nuclear physics. Procedures are described which are common to the modelling of reactions induced by light and medium-mass ions, and which are sufficient to calculate the effects of successive processes to any order.

1 Introduction

When two nuclei approach each other they may interact in several ways. In the first approximation they may be regarded as clusters of nucleons, and their primary interaction results from the inter-nucleon twobody force, which has an average effect found by folding the force over the internal configurations of the two clusters. However, these configurations are not static, and one or more rearrangement processes may occur during the time in which nuclei are together during a collision. As well as elastic scattering, with the projectile and the target remaining in their ground states, various kinds of non-elastic interactions may have time to operate.

Inelastic excitations may occur, for example when one or both of the nuclei are deformed or deformable, with the result that higher-energy states of the nuclei may become populated. Single-particle excitations are another kind of inelastic process, when a particle in one of the nuclei is excited during the reaction from its initial bound state to another state which may be bound or unbound. Nucleons may also transfer from nucleus to the other, either singly, or as the simultaneous transfer of two nucleons as a particle cluster.

In this paper I will consider some mathematical models sufficient to describe these processes, and the principal interest will be in calculating the effects of their occurring successively as multi-step processes. One-step processes have been traditional described with the Distorted Wave Born Approximation (DWBA), and although second-order DWBA expressions can be written down and computed, I shall be mainly concerned with coupled-channels formalisms, in order to predict the effects of multi-step processes to any or all orders.

The present work is in the framework of Direct Reaction theory, which attempts to solve the Schrödinger equation for a specific model of the components thought to be important in the reaction, and of their interaction potentials. In direct reaction theories, the phases describing the coherence of all components of the wave function are coherently maintained, and the potentials typically include imaginary components to model how flux is lost from the channels of the model to other channels. By contrast, a theory of compound-nucleus processes would make approximations as to the statistical distribution of the inelastic excitations. Direct reaction theory would describe these effects with an imaginary potential, with the argument that because the compound nucleus channels are incoherent with respect to each other, their effects back on the direct-reaction channels are also incoherent, and may hence be represented as a statistical loss of flux that occurs when the nuclei overlap each other to any significant extent.

Comprehensive accounts of the physical assumptions, methods and results of direct reaction theory is given in the papers by Tamura et al. [2], and in the books by Austern [3] and Satchler [33]. The aim of the present paper is to show how a large subset of the direct reaction mechanisms can be modelled in a general purpose computer program. For definiteness, I am following the methods used in the recent code FRESCO, while also mentioning, where appropriate, additional features that could well be included within its framework. Brief descriptions will also be given of alternative methods, and the relative merits of the different procedures will be discussed.

The code of ref.[34] has not been developed to include any special treatment of the long-range Coulomb mechanisms that are significant when heavy ions are incident on strongly-deformed nuclei. For methods of dealing with these processes efficiently, the reader is referred to refs. [35], [36], [4], and [5], [37].

The organisation of this paper is as follows. Section 2 will give a derivation of the coupled reaction channels (CRC) equations within the framework of the Feshbach formalism for direct reactions, and show how one-step and two-step DWBA (etc.) are special cases of the CRC equations. In both the CRC and DWBA formalisms, particular attention is paid to the treatment of the so-called ‘non-orthogonality terms’ which arise with couplings between different mass partitions.

The wave functions needed to specify scattering states and nuclear eigenstates are given in section 3. Single-nucleon wave functions are defined for both bound and unbound energies, and a method for solving the coupled-channels bound state eigen-problem is presented. Two-nucleon wave functions are described in both the centre-of-mass and independent-radii coordinates. Finally, the formulae are given for calculating the observable cross sections and polarisations in terms of the S-matrix elements of the scattering wave functions.

Section 4 specifies some of the different kinds of potentials that exist between nuclei, or can couple together the excited states of a single nucleus because of the interaction with its reaction partner. Details of the rotational model, single-particle excitations and particle transfers are given.

The methods used to solve the CRC equations are described in section 5, along with the procedures for calculating the transfer form factors in terms of a two-dimensional kernel function. Appendix A defines some of the notation and phase conventions used, while Appendix B summarises the more widely-known coupled-channels codes which have been written to solve problems in nuclear physics.

2 Coupled Reaction Channels Formalism

The coupled reaction channels (CRC) model of direct reactions in nuclear physics proceeds by constructing a model of the system wave function, and solving Schrödinger’s equation as accurately as possible within that model space. The model used here projects the complete wave function \overline{Ψ} onto a product ϕ_{i} ≡ ϕ_{ip} ∗ ϕ_{it} of projectile and target states with a wave function ψ_{i}(R_{i}) describing their relative motion:

P \overline{Ψ} ≡ Ψ = ∑_{i}^{N} ϕ_{i}ψ_{i}(R_{i})  (1)

The basis states ϕ_{ip} and ϕ_{it} can be bound states of their respective nuclei, or they may be discrete representations of continuum levels. In the former case we have a ‘bound state approximation’, and in the second case we have a ‘coupled discrete continuum channels’ [39, 40] (CDCC) approximation. The states ϕ_{i} can be in different mass partitions, or they can be different excited states of the projectile and/or the target in any one of the partitions. What is essential to the CRC framework is that there be a finite set (N say) of square-integrable basis states, as this leads to a finite set of equations coupling the channel wave functions ψ_{i}(R_{i}) as unknowns.

For a complete Hamitonian \overline{H} and total energy E, Schrödinger’s equation [\overline{H} − E]\overline{Ψ} = 0 becomes [H − E]Ψ = 0 in the model space with [6]

H = P \overline{H}P − P \overline{H}Q\frac{1}{Q\overline{H}Q − E − iϵ} Q\overline{H}P,  (2)

where Q ≡ 1 − P and ϵ is a positive infinitesimal quantity whose presence ensures that the excluded channels have a time-retarded propagator, and hence only remove flux from the model space. The second term as a whole describes the effects of the excluded channels on the model subspace P \overline{Ψ}. These effects could be, for example, from compound nucleus formation, which we have excluded from explicit consideration within direct reaction theory. In the absence of detailled knowledge of these effects, we construct our model Hamiltonian H using effective potentials which we believe approximate (in some average manner) the processes described by equation (2). The effective potentials will often be optical potentials with real and imaginary components fitted to some simpler kinds of reactions, and the effects of compound nucleus formation on these potentials is to contribute to their imaginary component.

The model Hamiltonian H for the CRC system can now be projected onto the individual basis states ϕ_{i}. If E_{i} is the asymptotic kinetic energy in the i’th channel, then the channel-projected Hamiltonian H_{i} satisfies

H_{i} − E_{i} = ⟨ϕ_{i}|H − E|ϕ_{i}⟩  (3)

and will be composed of a kinetic energy term and a diagonal optical potential. The ‘interaction potential’ V_{i} is then defined to be everything in H not included in H_{i}, so

H_{i} − E_{i} + V_{i} = H − E.  (4)

This construction gives V_{i} which have vanishing diagonal matrix elements ⟨ϕ_{i}|V_{i}|ϕ_{i}⟩ = 0.

2.1 Coupled Equation Set for N bound state pairs

If we take the model Schrödinger’s equation [H − E]Ψ = 0, and project separately onto the different basis states ϕ_{i}, we derive the set of equations

[E_{i} − H_{i}] ψ_{i}(R_{i}) = ∑_{j≠i} ⟨ϕ_{i}|H − E|ϕ_{j}⟩ ψ_{j}(R_{j}).  (5)

which couple together the unknown wave functions ψ_{i}(R_{i}). The matrix element ⟨ϕ_{i}|H − E|ϕ_{k}⟩ has two different forms, depending on whether we expand

H − E = H_{i} − E_{i} + V_{i} (the ‘post’ form)
= H_{j} − E_{j} + V_{j} (the ‘prior’ form).

Thus

⟨ϕ_{i}|H − E|ϕ_{j}⟩ = V_{ij}^{post} + [H_{i} − E_{i}]K_{ij}(post)  (6)
or = V_{ij}^{prior} + K_{ij}[H_{j} − E_{j}](prior)

where

V_{ij}^{post} ≡ ⟨ϕ_{i}|V_{i}|ϕ_{j}⟩, V_{ij}^{prior} ≡ ⟨ϕ_{i}|V_{j}|ϕ_{j}⟩, K_{ij} ≡ ⟨ϕ_{i}|ϕ_{j}⟩.  (7)

The overlap function K_{ij} = ⟨ϕ_{i}|ϕ_{j}⟩ in equation (6) arises from the well-known non-orthogonality between the basis states ϕ_{i} and ϕ_{j} if these are 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.

2.2 N-step DWBA

If the coupling interactions V_{i} in equation (6) 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 distorted wave Born approximation (DWBA). This approximation always feeds flux ‘forwards’ in the sequence 1 → 2 → ⋯ → 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

[E_{i} − H_{i}] ψ_{i}(R_{i}) = ∑_{j=1}^{j=i−1} ⟨ϕ_{i}|H − E|ϕ_{j}⟩ ψ_{j}(R_{j})  (8)

Initial channel:

[E_{1} − H_{1}] ψ_{1}(R_{1}) = 0

Second channel:

[E_{2} − H_{2}] ψ_{2}(R_{2}) = ⟨ϕ_{2}|H − E|ϕ_{1}⟩ ψ_{1}(R_{1})  (9)

If the prior interaction is used, the right hand side becomes

= ⟨ϕ_{2}|V_{1}|ϕ_{1}⟩ψ_{1} + ⟨ϕ_{2}|ϕ_{1}⟩[H_{1} − E_{1}]ψ_{1}
= ⟨ϕ_{2}|V_{1}|ϕ_{1}⟩ψ_{1} as ψ_{1} is on-shell.  (10)
= V_{21}^{prior}ψ_{1}  (11)

Final channel: (c = N + 1)

[E_{c} − H_{c}] ψ_{c}(R_{c}) = ∑_{j=1}^{j=c−1} ⟨ϕ_{c}|H − E|ϕ_{j}⟩ ψ_{j}(R_{j})  (12)

If the post interaction had been used for all the couplings to this last channel, then

[E_{c} − H_{c}] ψ_{c}(R_{c}) = ∑_{j=1}^{j=c−1}⟨ϕ_{c}|V_{c}|ϕ_{j}⟩ψ_{j} + [H_{c} − E_{c}] ∑_{j=1}^{j=c−1}⟨ϕ_{c}|ϕ_{j}⟩ψ_{j}  (13)

so

[E_{c} − H_{c}] χ_{c}(R_{c}) = ∑_{j=1}^{j=c−1} V_{cj}^{post}ψ_{j}  (14)

where

χ_{c}(R_{c}) = ψ_{c} + ∑_{j=1}^{j=c−1}⟨ϕ_{c}|ϕ_{j}⟩ψ_{j}
= ⟨ϕ_{c}|Ψ⟩

Note that, as all the ϕ_{j} are square-integrable and hence decay faster than r^{−1} at large radii, the ψ_{c} and χ_{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 χ_{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 [7] to avoid the non-orthogonality terms. It should be also clear that non-orthogonality terms will have to be evaluated if the DWBA is continued beyond second order.

2.3 Full CRC solution by iteration

There are a number of different ways of solving the CRC equations with the non-orthogonality terms present: for discussions of different approaches see refs. [8], [41] and the survey of ref.[ch. 3][33].

There are schemes available which can iterate all channels with an arbitrary choice of post or prior interactions for all the couplings. Define

θ_{ij} = 0 or 1 : presence of post on the j → i coupling,  (15)
so 1 − θ_{ij} = 1 or 0 : presence of prior.  (16)

The following iterative scheme [42] (n=1,2,..) on convergence then solves the CRC equations (5):

For n =0, start with

ψ_{i}^{(0)} = δ(i, i_{0})ψ_{elastic}  (17)
δS_{i}^{(0)} = δψ_{i}^{(0)} = 0  (18)

For n = 1 → N + 1 (for N-step DWBA) solve

[H_{i} − E_{i}]χ_{i}^{(n)} + S_{i}^{(n−1)} = 0  (19)

with

S_{i}^{(n−1)} = ∑_{j}[θ_{ij}V_{ij}^{post} + (1 − θ_{ij})V_{ij}^{prior}] ψ_{j}^{(n−1)} − δS_{i}^{(n−1)}  (20)

then calculate for subsequent iterations

δψ_{i}^{(n)} = ∑_{j} θ_{ij}⟨ϕ_{i}|ϕ_{j}⟩ψ_{j}^{(n−1)}  (21)
δS_{i}^{(n)} = ∑_{j}(1 − θ_{ij})⟨ϕ_{i}|ϕ_{j}⟩[S_{j}^{(n−1)} + [H_{j} − E_{j}]δψ_{j}^{(n)}]  (22)
ψ_{i}^{(n)} = χ_{i}^{(n)} − δψ_{i}^{(n)}  (23)

This scheme avoids numerical differentations except in an higher-order correction to δS_{i} that arises in some circumstances.

When the non-orthogonality terms are included properly, it becomes merely a matter of convenience whether post or prior couplings are used, for one, two, and multistep calculations. The equivalence of the two coupling forms can be confirmed in practice (see, for example, refs.[42], and [9]), and used as one check on the accuracy of the numerical methods employed.

3 Wave Functions for Scattering and Bound States

In order to describe details of the nuclear transitions realistically, it is necessary to specify in sufficient detail the initial and final states of the nuclei involved. To start with, the excitation energies, spins and parities of all the states in each mass partition need to be specified, along with the nuclear masses, charges and relative Q-values of the partitions. Each pair projectile and target excited states is then a distinct channel with its own scattering wave function and boundary conditions. The initial projectile and target states will select one such channel as the ‘incoming channel’, with its boundary conditions specifying an incoming plane wave. All channels (including the incoming channel) will have outgoing spherical waves. Particular attention must be given to the consistent placement of i^{L} factors in these definitions.

The individual nuclear states are then specified in sufficient detail for the particular reaction mechanisms involved. It is not necessary to specify the full quantum mechanical states of all the nucleons in the nucleus, but rather, only the states of those changed in the reactions being considered. In particular, one and two-nucleon wave functions will have to be described, if those nucleons are to be transferred to other nuclei. If a nuclear state consists of a particle of spin s bound outside a nucleus with possible core states ϕ_{I}, then the bound state radial wave functions u_{ℓsjI}(r) will have to be found by solving a coupled-channels set of equations for negative energy eigen-solutions. If the particle is not bound, in the other hand, then its continuum range of energies must be discretised into a finite collection of ‘bin’ states which can be scaled to unit normalisation. If the nuclear state consists of two particles of intrinsic spins s_{1} and s_{2} outside a core, then it is usually specified by a shell-model or by a Sturmian-basis calculation in terms of the independent coordinates r_{1} and r_{2}. To calculate transfer rates, however, the two-particle wave functions need to be given in terms of the collective coordinates (usually r = \frac{1}{2}(r_{1} + r_{2}) and ρ = r_{1} − r_{2}). In order to use the states in a reaction calculation, therefore, equations are given for the transformation from the independent coordinates.

When we have calculated the scattering wave functions, or at least their asymptotic parts in terms of their S-matrix elements, we can find the cross sections for each outgoing pair of projectile and target states in each partition. Furthermore, if the initial projectile has non-zero spin J_{p}, then the effect on these cross sections of polarisation of the projectile is specified by the tensor analysing powers T_{kq} (for 1 ≤ k ≤ 2J_{p} and 0 ≤ q ≤ k). Integrated cross sections and fusion polarisations can also be found using the S-matrix elements.

3.1 Total wave function

In each partition κ of the system into a projectile of mass A_{κp} and a target of mass A_{κt}, the coupling order is

L + J_{p} = J; J + J_{t} = J_{T} ,  (24)

which may be defined by writing

Ψ_{κJ_{T}}^{M_{T}} = | (L J_{p})J, J_{t}; J_{T} ⟩  (25)

where J_{p} = projectile spin, J_{t} = target spin, L = orbital partial wave, and J_{T} = total system angular momentum.

The set {κpt, (L J_{p})J, J_{t}; J_{T} } will be abbreviated by the single variable α. Thus, in each partition,

Ψ_{κJ_{T}}^{M_{T}} (R_{κ}, ξ_{p}, ξ_{t}) = ∑ ; \begin{matrix}LJ_{p}JJ_{t} \\ Mμ_{p}M_{J}μ_{t} \end{matrix} ϕ_{J_{p}}(ξ_{p}) ϕ_{J_{t}}(ξ_{t}) i^{L}Y_{L}^{M}(R̂_{κ})\frac{1}{R_{κ}} f_{(LJ_{p})J,J_{t}}^{κJ_{T}}(R_{κ})

⟨LMJ_{p}μ_{p}|JM_{J}⟩ ⟨JM_{J}J_{t}μ_{t}|J_{T} M_{T} ⟩  (26)

here R_{κ} = radial coordinate from the target to the projectile in partition κ, ξ_{p} = internal coordinates of projectile, ξ_{t} = internal coordinates of target, and

f_{(LJ_{p})J,J_{t}}^{κJ_{T}}(R) ≡ f_{α}(R)  (27)

are the radial wave functions. The i^{L} factors arise from the spherical Bessel expansion of the incoming plane wave. Some formalisms include extra powers of i in the equation (26), in order to make the coupling interactions V_{α:α′} real. Inelastic coupling interactions can be made real (for integer-valued spins, at least), by including a factor i^{J_{p}+J_{t}} in the definition of Ψ, and transfer couplings can be made real by including a factor i^{ℓ} for orbital angular momentum ℓ of the bound particle state in this partition. In a general purpose code [34], however, there may be clashes between these different conventions. The ground state of ^{7}Li, for example, would have a factor of i^{3/2} on the rotational model convention, but a factor of i^{1} if the state were regarded as a ℓ = 1 bound state of a α core and a triton cluster. It seems simplest, therefore, to omit these addition phase factors completely. The coupling interactions can very often be made real, nevertheless, if the i^{L} factors are included explicitly in the CRC equations, as in the next section.

The wave function Ψ could also have been defined using the ‘channel spin’ representation (as in [43]) Ψ = |L, (J_{p}J_{t})S; J_{T} ⟩, which is symmetric upon projectile → target interchange except for a phase factor (−1)^{S−J_{p}−J_{t}}. This would simplify the subsequent description of the coupling elements in section 4, as the formulae for projectile mechanisms and target mechanisms would differ only by this phase factor. However, the channel spin representation has the disadvantage that the projectile spin-orbit force is not diagonal in this basis. This would not matter if coupled-channels solutions were always sought, but one of the advantages of sometimes solving the CRC equations iteratively is that the DWBA solutions of first and second order (etc.) may be obtained. In order for the partially-iterated CRC solutions to reproduce the results of DWBA codes, it is necessary to treat spin-orbit forces without approximation, and since spin-orbit forces almost always are those of the projectile, the asymmetric representation of channel (24) is advisable.

Identical Nuclei If one partition (κ say) is identical to another κ except that the projectile and target nuclei are exchanged, then the total wave function should be formed from (1 + πP_{κκ′}) times the above expression, where π = ±1 is the intrinsic parity of the two nuclei under exchange. A simple method of dealing with this exchange is to first form the wave function of equation (26), and then operate with (1 + πP_{κκ′}) on both the wave functions and the S-matrix elements, before cross sections are calculated. This is equivalent to the replacement

f_{α}(R) ← f_{α}(R) + c_{α,α′}f_{α′}(R)  (28)

where

c_{α,α′} = π(−1)^{L}δ_{L,L′}(−1)^{J_{T} +L−J−J′} ĴĴ′W (J_{p}L′J_{T} J_{t}; JJ′)  (29)

with α′ = |(L′J_{t})J′, J_{p}; J_{T} ⟩.

3.2 Coupled equations

The CRC equations are in many cases of the form

[E_{κpt} − T_{κL}(R_{κ}) − U_{κ}(R_{κ})] f_{α}(R_{κ}) = ∑_{α′,Γ>0} i^{L′−L} V_{α:α′}^{Γ}(R_{κ′})f_{α′}(R_{κ′})

+ ∑_{α′,κ′≠κ} i^{L′−L} ∫_{0}^{R_{m}} V_{α:α′}((R_{κ}), R_{κ′})f_{α′}(R_{κ′})dR_{κ′} (30)

where the kinetic energy term is

T_{κL}(R) = −\frac{h̄^{2}}{2μ_{κ}} (\frac{d^{2}}{dR^{2}} − \frac{L(L + 1)}{R^{2}}) ,  (31)

U_{κ}(R_{κ}) is the diagonal optical potential with nuclear and Coulomb components, and R_{m} is a radius limit larger than the ranges of U_{κ}(R_{κ}) and of the coupling terms. The equations (30) are in their most common form: they become more complicated when non-orthogonalities are included by the method of section 2.3. The V_{α:α′}^{Γ}(R_{κ′}) are the local coupling interactions of multipolarity Γ, and the V_{α:α′}(R_{κ}, R_{κ′}) are the non-local couplings between mass partitions that arise from particle transfers.

For incoming channel α_{0}, the f_{α}(R_{κ}) satisfy the boundary conditions

f_{α}(R_{κ′}) =_{R_{κ}>R_{m}} \frac{i}{2} [δ_{αα_{0}}H_{Lη_{α}}^{(−)}(K_{α}(R_{κ})) − S_{α_{0}α}H_{Lη_{α}}^{(+)}(K_{α}(R_{κ}))]  (32)

where H_{Lη}^{(−)} and H_{Lη}^{(+)} are the Coulomb functions [44] with incoming and outgoing boundary conditions respectively. The asymptotic kinetic energies are

E_{κpt} = E + Q_{κ} − ϵ_{p} − ϵ_{t}  (33)

for excited state energies ϵ_{p}, ϵ_{t} and Q-value Q_{κ} in partition κ, and

K_{α} = \sqrt{[\frac{2μ_{κ}}{h̄^{2}E_{κpt}} ]}  (34)

where μ_{κ} = A_{κp}A_{κt}/(A_{κp} + A_{κt}) is the reduced mass in the channel with partition κ, and

η_{α} = \frac{2μ_{κ}}{h̄^{2}} \frac{Z_{κp}Z_{κt}e^{2}}{2K_{α}}  (35)

is the Sommerfeld parameter for the Coulomb wave functions.

3.3 Single-nucleon states

If ϕ_{JM}(ξ) is a core+particle bound state, then for coupling order |(ℓs)j, I; JM⟩ , the wave function is

ϕ_{JM}(ξ_{c}, r) = ∑_{ℓjI} A_{ℓsj}^{jIJ} [ϕ_{I}(ξ_{c})φ_{ℓsj}(r)]_{JM}  (36)
= ∑_{ℓjI,mμm_{s}m_{ℓ}} A_{ℓsj}^{jIJ} ⟨jmIμ|JM⟩ ϕ_{Iμ}(ξ_{c}) ⟨ℓm_{ℓ}sm_{s}|jm⟩ Y_{ℓ}^{m_{ℓ}}(r̂)ϕ_{s}^{m_{s}} \frac{1}{r}u_{ℓsjI}(r)

where ξ_{c} = core internal coordinates, ϕ_{Imu}(ξ_{c}) = core internal state, ϕ_{s}^{m_{s}} = particle internal spin state, u_{ℓsjI}(r) = particle core radial wave function, and A_{ℓsj}^{jIJ} is the coefficient of fractional parentage.

3.3.1 Bound States

If the single-particle is bound at negative energy E around the core, then its wave function may be found as the eigen-solution of a given potential:

[T_{ℓ}(r) + V (r) + ϵ_{I} − E]u_{ℓsjI}(r) + ∑_{ℓ′j′I′, Γ>0} V_{ℓsjI:ℓ′sj′I′}^{Γ}(r)u_{ℓ′sj′I′}(r) = 0  (37)

with boundary conditions u_{ℓsjI}(0) = 0 and, as r ≥ R_{m}, of u_{ℓsjI}(r) ∝ W_{ℓ}(k_{I}r) where W_{ℓ}(ρ) is the Whittaker function and k_{I}^{2} = 2μ|E − ϵ_{I}|/h̄^{2} is the asymptotic wave number.

If the core cannot be excited, then these coupled equations reduce to one uncoupled equation, but solving this equation can still be regarded as a special case of the coupled bound state problem. Eigen-solutions can be found by solving either for the bound state energy E, or by varying the depth of the binding potential. In general, however, we can choose to vary any multipole of any part of the binding potentials (except the Coulomb component), so one method of solving the full coupled bound-state problem will be outlined below.

To define the phase (±1) of the overall wave function, some convention has to be adopted. One component (say that around a core I = 0 state) can be set to either positive towards the origin (r → 0), or positive towards large distances (r → ∞). The former choice is made in the FRESCO code, following the Mayer-Jensen phase convention, which is also used for harmonic oscillator wave functions in many structure calculations.

3.3.2 Solution of the Coupled-Channels Eigenvalue Problem

When, for example, the problem is to find the bound state of a particle in a deformed potential, then several channels with different angular momenta will be coupled together. There are various techniques for calculating the wave functions of these bound states: for a review see ref. [10]. The Sturmian expansion method [45] can be used, or the coupled equations can be solved iteratively. The Sturmian method has the advantage that all solutions in the deformed potential are found, where sometimes the iterative method has difficulty in converging to a particular solution if there are other permitted solutions near in energy. The iterative method has the advantage that the radial wave functions (once found) are subject only to the discretisation error for the Schrödinger’s equation, and are not dependent on the (time-consuming) diagonalisation of large matrices, often of the order of 100 or more. As they satisfy the correct boundary conditions independently of the size of a basis-state set, the radial wave functions of the iterative method therefore more accurately reflect the details of the coupling potentials and of the core excitation energies. As nuclear reactions are often confined to the surface region, it is important to satisfy the exterior boundary conditions as accurately as possible.

A method for solving the uncoupled eigenstate problem has to be included in a reaction code in any case, and since it can be generalised as described in this section to solving the coupled problem, it seems a worthwhile facility to have available. Bound states from a previous Sturmian solution can still be included as explicit linear combinations of the single particle (uncoupled) basis states used in the Sturmian expansion.

The general problem of finding eigen-solutions of a set M coupled-channels equations can be represented as the problem of finding λ such that the equations

[\frac{d^{2}}{dr^{2}} − \frac{ℓ_{i}(ℓ_{i} + 1)}{r^{2}}] ψ_{i}(r) + ∑_{j=1}^{M} [U_{ij}(r) + λV_{ij}(r)] ψ_{j}(r) = 0  (38)

with boundary conditions

ψ_{i}(R) = a_{i}W_{ℓ_{i}η_{i}}(k_{i}R)  (39)
ψ_{i}(R + δR) = a_{i}W_{ℓ_{i}η_{i}}(k_{i}(R + δR))  (40)
ψ_{i}(0) = 0  (41)

(with k_{i}^{2} ≡ κ_{i}^{2} + θλ and η_{i} ≡ nu_{i}/(2k_{i})) for given partial waves ℓ_{i}, fixed potentials U_{ij}(r), variable potentials V_{ij}(r), matching radius R, and Coulomb proportionality constants ν_{i}. The energy constants κ_{i}^{2} are the asymptotic components of the diagonal U_{ii}(r), and θ is the asymptotic component of the diagonal V_{ii}(r) (assumed all equal).

The solution begins by constructing the trial integration functions for a trial value of λ, on either side of an intermediate matching point r = ρ:

f_{i;j}^{in}(r) by integrating r from h to ρ, starting with f_{i;j}^{in}(h) = δ_{i,j} h^{ℓ_{i}+1}/(2ℓ_{i} + 1)!!, and

f_{i;j}^{out}(r) by integrating r from R in to ρ, starting with f_{i;j}^{in}(R) = δ_{i,j} W_{ℓ_{i}η_{i}}(k_{i}(R + δR)).

The intermediate point r = ρ should be chosen where the wave functions are oscillatory, to avoid having to integrate outwards in the classically forbidden region.

The solution is therefore

ψ_{i}(r) = \{ \begin{matrix}∑_{j} b_{j}f_{i;j}^{in}(r) for r < ρ \\ ∑_{j} c_{j}f_{i;j}^{out}(r) for r ≥ ρ, \end{matrix}  (42)

and the matching conditions are the equality of the two expressions and their derivatives at r = ρ. The normalisation is still arbitrary, so we may fix c_{1} = 1. In general the equations (38) have no solution as λ is not exactly an eigenvalue. The method therefore uses the discrepancy in the matching conditions to estimate how λ should be changed to λ + δλ to reduce that discrepancy, and iterates this process to reduce δλ.

Thus at each iteration we first solve as simultaneous equations the 2M-1 of the matching conditions

∑_{j} b_{j}f_{i;j}^{in}(ρ) = ∑_{j} c_{j}f_{i;j}^{out}(ρ) for all i  (43)
∑_{j} b_{j}f_{i;j}^{in}(ρ)′ = ∑_{j} c_{j}f_{i;j}^{out}(ρ) for all i ≠ 1  (44)

along with c_{1} = 1 for the 2M unknowns b_{i}, c_{i}. If the function ψ_{i}(r) is then constructed using equation (42), there will be a discrepancy as

ψ′_{in} ≡ ψ_{1}(r)|_{r<ρ} ≠ ψ′_{out} ≡ ψ_{1}(r)|_{r>ρ},  (45)

and this difference will generate δλ via

δλ ∑_{ij}^{M} ∫_{0}^{R} ψ_{i}(r)V_{ij}(r)ψ_{j}(r)dr = ψ_{1}(ρ)[ψ′_{out} − ψ′_{in}].  (46)

It is necessary while iterating in this manner to monitor the number of nodes in one or more selected components of the wave function, as in general a given potential will have different eigensolutions with different numbers of radial nodes. When the iterations have converged to some accuracy criterion on the size of δλ, the set of wave functions can be normalised in the usual manner:

∑_{i}^{M} ∫_{0}^{∞} |ψ_{i}(r)|^{2}dr = 1  (47)

and perhaps some of the components i omitted if their contribution to this norm is below some preset threshold.

3.3.3 Continuum States

If the initial and/or final single-particle states of a transfer step are unbound E − ϵ > 0, the use of single energy eigenstates ϕ_{k}(r) will result in calculations of the transfer form factors which will not converge, as the continuum wave functions do not decay to zero as r → ∞ sufficiently fast as to have square norms. One way [39], [40] of dealing with this divergence is to take continuum states not at a single energy, but averaged over a range of energies. These ‘bin’ states that result are square integrable, and if defined as

Φ(r) = \sqrt{\frac{2}{πN}} ∫_{k_{1}}^{k_{2}} w(k)ϕ_{k}(r)dk  (48)
with N = ∫_{k_{1}}^{k_{2}} |w(k)|^{2}dk  (49)

for some weight function w(k), then they are normalised ⟨Φ|Φ⟩ = 1 provided a sufficiently large maximum radius for r is taken, and that the ϕ_{k} are eigensolutions of a potential which is energy-independent. They are orthogonal to any bound states, and are orthogonal to other bin states if their energy ranges do not overlap. The construction can be easily generalised to give coupled-channels bin wave functions.

The weight function w(k) is best chosen ([40] p. 148) to include some of the effects known to be caused by the variation of ϕ_{k}(r) within the bin range k_{1} ≤ k ≤ k_{2}. If w(k) = exp(−iδ_{k}), where delta_{k} is the scattering phase shift for ϕ_{k}(r), then it includes the effects of the overall phase variations of ϕ_{k}, at least in the DWBA limit. If, however, w(k) = exp(−iδ_{k}) sin δ_{k} ≡ T_{k}^{∗}, where T_{k} is the T -matrix element for ϕ_{k}(r), then it includes in addition a scale factor which is useful if the |T_{k}| varies significantly, as it does, for example, over resonances. Both choices result in a real-valued wave function Φ(r) (for real potentials), which is computationally advantageous.

If the maximum radius (R_{m} say) is not sufficiently large, then the wave functions Φ will not be normalised to unity by the factors given in equation (48). The rms radius of a bin wave function increases as the bin width k_{2} − k_{1} decreases, approximately as 1/(k_{2} − k_{1}). These bin constructions can be used to describe the narrow resonant wave functions of say the 3^{+} state in ^{6}Li, or the 7/2^{−} state in ^{7}Li, but these states will require a large limiting radius R_{m} unless the w(k) = T_{k}^{∗} weighting factor is used to

R_{m} δE = 0.1 MeV  δE = 0.5 MeV
(fm)  e^{−iδ_{k}}  T_{k}^{∗}  e^{−iδ_{k}}  T_{k}^{∗}
10  0.105  0.923  0.108  0.977
20  0.109  0.939  0.132  0.985
40  0.112  0.951  0.258  0.993
80  0.122  0.971  0.411  0.996
160  0.142  0.986  0.614  0.998

Figure 1: Normalisations of a continuum bin state. For this 3^{+} state in ^{6}Li at 0.71 MeV, Saxon-Woods potentials were used with V = 77.05 MeV, R = R_{c} = 1.2 * 4^{1/3} fm., and a = 0.65 fm.

emphasise the increase in the interior wave function over the resonance. The 3^{+} state in^{6}Li at 0.71 MeV, for example, for which the resonance width is approximately 40 keV, yields the normalisations shown in 3.3.3. It can be seen that without a scale factor which emphasises the resonance peak, very large radii R_{m} will be needed to obtain unit normalisation.

3.4 Two-particle bound states

3.4.1 Centre-of-mass coordinates

If ϕ_{JT} (ξ_{c}, r, ρ) is a two-particle bound state with total spin J and isospin T, then for coupling order |{ L, (ℓ, (s_{1}s_{2})S)j}J_{12}, I; J⟩ we have

ϕ_{JM} = ∑ ; \begin{matrix}LℓS \\ jJ_{12}I \end{matrix}  A_{LjJ_{12}}^{J_{12}IJ} ϕ_{Iμ_{I}} (ξ_{c}).ϕ_{s_{1}}^{σ_{1}} ϕ_{s_{2}}^{σ_{2}} Y_{L}^{Λ}(r̂) Y_{ℓ}^{μ}(ρ̂) \frac{1}{rρ}u_{12}(r, ρ)

⟨J_{12}M_{12}Iμ_{I}|JM⟩⟨LΛjm_{12}|J_{12}M_{12}⟩⟨ℓμSΣ|jm_{12}⟩⟨s_{1}σ_{1}s_{2}σ_{2}|SΣ⟩ (50)

where A_{LjJ_{12}}^{J_{12}IJ} is the coefficients of fractional parentage, and ϕ_{s_{1}}^{σ_{1}} ϕ_{s_{2}}^{σ_{2}} are the intrinsic spins of the two particles.

Note that two neutron transfer can be viewed as the transfer of a ‘structured particle’ (ℓ, (s_{1}s_{2})S)j, and then becomes similar to single-particle transfers of above.

The radial wave function u_{12}(r, ρ) can be given either as a cluster product of single-particle wave functions u_{12}(r, ρ) = Φ_{L}(r)ϕ_{ℓ}(ρ), or input directly as a two-dimensional distribution e.g. from a Faddeev bound-sate calculation, or calculated from the correlated sum of products of single-particle states, as in the next section.

3.4.2 Independent Coordinates

Two-particle states from shell-model calculations or from Sturmian-basis calculations [11], and are then usually described by means of the |r_{1}, r_{2}⟩ coordinates, and then transformed internally into the centreof-mass coordinates |r, ρ⟩ of equation (50) using r_{i} = x_{i}r+y_{i}ρ. For equal mass particles, x_{1} = x_{2} = 1, and y_{1} = −y_{2} = \frac{1}{2}. The second description is as

φ_{12}(r_{1}, r_{2}) = ∑_{i} c_{i} |(ℓ_{1}(i), s_{1})j_{1}(i), (ℓ_{2}(i), s_{2})j_{2}(i); J_{12}T ⟩  (51)
→ ∑_{u} c_{i} ∑_{LℓSj} |L, (ℓ, (s_{1}s_{2})S)j; J_{12}T ⟩ ϕ_{L(ℓS)j}^{J_{12}T,i}(r, ρ)  (52)

The transformation of the i’th component in the cluster basis is

ϕ_{L(ℓS)j}^{J_{12}T,i}(r, ρ) = ⟨L, (ℓ, (s_{1}s_{2})S)j; J_{12}T | (ℓ_{1}(i), s_{1})j_{1}(i), (ℓ_{2}(i), s_{2})j_{2}(i); J_{12}T ⟩  (53)
×⟨[Y_{L}(r̂)Y_{ℓ}(ρ̂)]_{λ} | [φ_{ℓ_{1}s_{1}j_{1}}(r_{1})φ_{ℓ_{2}s_{2}j_{2}}(r_{2})]_{J_{12}T} ⟩  (54)

where (suppressing the i indices for clarity)

⟨L, (ℓ, (s_{1}s_{2})S)j; J_{12}T |(ℓ_{1}, s_{1})j_{1}, (ℓ_{2}, s_{2})j_{2}; J_{12}T ⟩ =

∑_{λ}  λ̂Ŝĵ_{1}ĵ_{2} ( \begin{matrix}ℓ_{1} & ℓ_{2} & λ \\ s_{1} & s_{2} & S \\ j_{1} & j_{2} & J_{12} \end{matrix} ) \frac{1 + (−1)^{ℓ+S+T}}{\sqrt{2(1 + δ_{ℓ_{1},ℓ_{2}}δ_{j_{1},j_{2}})}} ĵλ̂W (LℓJ_{12}S; λj)(−1)^{ℓ+L−λ}. (55)

The radial overlap integral can be derived by means of harmonic-oscillator expansions [12], with the Bayman-Kallio expansion [13] or using the Moshinsky solid-harmonic expansion[71]. This last method gives

K_{ℓL:ℓ_{1}ℓ_{2}}^{λ}(r, ρ) = ⟨[Y_{L}(r̂)Y_{ℓ}(ρ̂)]_{λ} | [φ_{ℓ_{1}}(r_{1})φ_{ℓ_{2}}(r_{2})]^{λ}⟩
= ∑_{n_{1}n_{2}} ( \begin{matrix}2ℓ_{1} + 1 \\ 2n_{1} \end{matrix} )^{\frac{1}{2}} ( \begin{matrix}2ℓ_{w} + 1 \\ 2n_{2} \end{matrix} )^{\frac{1}{2}} (x_{1}r)^{ℓ_{1}−n_{1}}(y_{1}ρ)^{n_{1}}(x_{2}r)^{n_{2}}(y_{2}ρ)^{ℓ_{2}−n_{2}}

× ∑_{Q} q_{ℓ_{1}ℓ_{2}}^{Q}(r, ρ) (2Q + 1) ℓ̂_{1}ℓ̂_{2}ℓ̂_{1} − n_{1}ℓ̂_{2} − n_{2} L̂ℓ̂

× ∑_{Λ_{1}Λ_{2}} ( \begin{matrix}ℓ_{1} − n_{1} & n_{2} & Λ_{1} \\ 0 & 0 & 0 \end{matrix} ) ( \begin{matrix}ℓ_{w} − n_{2} & n_{1} & Λ_{2} \\ 0 & 0 & 0 \end{matrix} ) ( \begin{matrix}Λ_{1} & L & Q \\ 0 & 0 & 0 \end{matrix} ) ( \begin{matrix}Λ_{2} & ℓ & Q \\ 0 & 0 & 0 \end{matrix} )

×(−1)^{ℓ_{1}+ℓ_{2}+L+Λ_{2}}(2Λ_{1} + 1)(2Λ_{2} + 1)W (Λ_{1}LΛ_{2}ℓ; Qλ) ( \begin{matrix}ℓ_{1} − n_{1} & n_{2} & Λ_{1} \\ n_{1} & ℓ_{2} − n_{2} & Λ_{2} \\ ℓ_{1} & ℓ_{2} & λ \end{matrix} ) . (56)

where ( \begin{matrix}a \\ b \end{matrix} ) is the binomial coefficient (see Appendix A).

The kernel function q_{ℓ_{1}ℓ_{2}}Q(r, ρ) which appears in this expression is the Legendre expansion of the product of the two radial wavefunctions in terms of u, the cosine of the angle between r and ρ:

q_{ℓ_{1},ℓ_{2}}^{Q}(r, ρ) = \frac{1}{2} ∫_{−1}^{+1} \frac{u_{ℓ_{1}}(r_{1})}{r_{1}}^{ℓ_{1}+1} \frac{u_{ℓ_{2}}(r_{2})}{r_{2}}^{ℓ_{2}+1} P_{Q}(u)du  (57)

3.5 Scattering Amplitudes

The Rutherford amplitude for pure Coulomb scattering (with no e^{2iσ_{0}} factor) is

F_{c}(θ) = −\frac{η}{2k} \frac{exp(−2iη ln(sin θ/2))}{sin^{2} θ/2}  (58)

The Legendre coefficients for the scattering to the projectile state J′_{p} and target state J′_{t} from initial projectile state J_{p} and target state J_{t} are given by

A_{m′M′;mM}^{L′} = ∑_{L,J,J′,J_{T}} ⟨L0J_{p}m|Jm⟩⟨JmJ_{t}M|J_{T} M_{T} ⟩

⟨L′M_{L′}J′_{p}m′|J′M_{L′} + m′⟩⟨J′M_{L′} + m′J′_{t}M′|J_{T} M_{T} ⟩

\frac{4π}{k} \sqrt{\frac{k′}{μ′} \frac{μ}{k}}e^{i(σ_{L}−σ_{0})}e^{i(σ′_{L′}−σ′_{0})}

(\frac{i}{2}) [δ_{α,α′} − S_{α,α′}^{J_{T}}] \sqrt{\frac{2L + 1}{4π}} Y_{c}(L′, M_{L′})  (59)

where Y_{c}(L, M) is the coefficient of P_{L}^{|M|}(cos θ)e^{iMϕ} in Y_{L}^{M}(θ, ϕ), σ_{L} = arg Γ(1 + L + iη) is the Coulomb phase shift, α′ refers to the primed values L′J′_{p}J′_{t}k′μ′ etc., and α refers to the unprimed values LJ_{p}J_{t}kμ.

For each outgoing channel J′_{p}, J′_{t}, we may then calculate the angular-dependent scattering amplitudes

f_{m′M′:mM}(θ) = δ_{J_{p},J′_{p}}δ_{J_{t},J′_{t}}F_{c}(θ) + ∑_{L′} A_{m′M′:mM}^{L′}P_{L′}^{m′+M′−m−M}(cos θ)  (60)

in terms of which the differential cross section is

\frac{dσ(θ)}{dΩ} = \frac{1}{(2J_{p} + 1)(2J_{t} + 1)} ∑_{m′M′mM} |f_{m′M′:mM}(θ)|^{2} .  (61)

The near-side and far-side decompositions [14] of this cross section are defined by the same process, with P_{L}^{M}(u) replaced by \frac{1}{2}[P_{L}^{M}(u) ± 2i/πQ_{L}^{M}(u)] respectively. The Coulomb scattering of equation (58) is included in the near-side component [15].

The spherical tensor analysing powers T_{kq} describe how the outgoing cross section depends on the incoming polarisation state of the projectile. If the spherical tensor τ_{kq} is an operator with matrix elements

(τ_{kq})_{mm′′} = \sqrt{2k + 1}⟨J_{p}mkq|J_{p}m′′⟩,

we have

T_{kq}(θ) = \frac{T r(fτ_{kq}f^{+})}{T r(ff^{+})}  (62)
= k̂ \frac{∑_{m′M′mM} f_{m′M′:mM}(θ)^{∗}⟨J_{p}mkq|J_{p}m′′⟩f_{m′M′:m′′M}(θ)}{∑_{m′M′mM} |f_{m′M′:mM}(θ)|^{2}}  (63)

The polarisations in the ‘transversity frame’ [16] are then

T T_{10} = \sqrt{2}iT_{11}  (64)
T T_{20} = −\frac{1}{2}(T_{20} + \sqrt{6}T_{22})  (65)
T T_{30} = −\frac{1}{2}(\sqrt{3}iT_{31} + \sqrt{5}iT_{33}).  (66)

The S-matrix elements can also be used to directly calculate the integrated cross sections

σ = ∫_{4π} \frac{dσ(θ)}{dΩ}dΩ  (67)

to give

σ = \frac{1}{k^{2}} \frac{k′}{μ′} \frac{μ}{k} \frac{4π}{(2J_{p} + 1)(2J_{t} + 1)} ∑_{J_{T} α α′}(2J_{T} + 1) |S_{α,α′}^{J_{T}} |^{2} .  (68)

The fusion cross section is defined as that amount of flux which leaves the coupled-channels set because of the imaginary parts of the optical potentials. If the incoming projectile is not spherical, then the fusion rate will depend on its orientation, and hence on the magnetic substate quantum number m. One can therefore define the fusion polarisation as the distribution σ_{m}^{fus}, which can be calculated from the S-matrix elements as

σ_{m}^{fus} = \frac{π}{k^{2}} ∑_{J_{T} ≥m}(2J_{T} + 1)

×[1 − \frac{1}{2J_{t} + 1} ∑_{J′_{p}L′Mω} |∑_{LJ}⟨J_{p}mJ − m|L0⟩⟨JmJ_{t}M|J_{T} M + m⟩e^{iσ_{L}}S_{α,α′}^{J_{T} ω}|^{2}] (69)

where ω is the parity (±1) of the coupled-channels set for each total angular momentum J_{T}

Partial Wave Interpolation: Heavy ion reactions typically involve a range of partial waves L up to several hundred or more, especially when Coulomb excitations dominate the highest partial waves. In such cases it is often advantageous to solve the coupled channels sets (30) for, say, every n’th value of J_{T}, and interpolate the intermediate values. Different values of n can be used in different reaction regions: n can be small (1 or 2) for the grazing partial waves, and up to 5 or 10 for the Coulomb-dominated peripheral processes, and can be adjusted for the required balance between speed and accuracy.

This interpolation may be performed on the S-matrix elements themselves, or on the Legendre amplitudes of equation (59) In this second method (that used in ref. [34]), cubic spline interpolations are used. The main factor limiting the accuracy of this process is that the rate of change with J_{T} of the Coulomb phase shifts exp i(σ_{L} +σ′_{L′}) does not diminish as J_{T} increases. For that reason, it is advisable to interpolate not on the A^{L′} of equation (59), but on a Ã^{L′} defined with a revised phase shift factor exp i(σ_{L} −σ′_{L′}). Since L and L′ both tend to be near J_{T}, it is only the difference the phase shifts which limits the accuracy of the interpolation. It will therefore be more accurate for smaller projectile and target spins, and incoming and outgoing channels with similar Sommerfeld parameter η (equation 35).

4 Coupling Interactions

When two nuclei interact, a variety of kinds of elastic and inelastic potentials may be needed to describe their interaction. As well as the scalar nuclear attractions and scalar Coulomb repulsions, if either of the nuclei has spin J ≠ 0, then there can be higher-order tensor interactions which couple together the spin and the orbital motion. If a nucleus has spin J ≥ \frac{1}{2}, then there can be a spin-orbit component V_{ls}(R)2l·J in the Hamiltonian H. and if its spin is one or greater (J ≥ 1), there can be tensor forces of various kinds. The most commonly used tensor force is a T_{r} potential of the form V_{T r}(R)R_{2}(R, R).S_{2}(J, J). Similar tensor forces are also generated if the projectile and target spins coupled together can reach J_{p} + J_{t} ≥ 1: such is the case with the tensor force between the proton and the neutron within the deuteron.

Inelastic potentials (4.2) arise when one or both of the nuclei have permanent deformations (as seen in their intrinsic frame), or are vibrationally deformable. The inelastic potentials which come from rotating a permanently deformed nucleus are described in the Hamiltonian by terms of the form

V_{λ} = ∑_{μ} V_{λ}(R)D_{μ0}^{λ}Y_{λ}^{μ}(R̂)  (70)

where the form factors V_{λ}(R) have both nuclear and Coulomb components for angular momentum transfers λ. Their nuclear component is approximately proportional to the derivative of the scalar potential between the two reaction partners. Simultaneous excitations of both nuclei are also possible (see e.g. [17]), but have not been included in the present code. Vibrational excitations of a nucleus have more complicated form factors in general [1], but can still be expanded in the form of equation (70). For the more intricate level schemes of strongly-deformed nuclei, it will in general be necessary for each allowed transition to have its own transition rate specified independently of a particular rotational or vibrational model.

Inelastic potentials also arise when one of the nuclei can be decomposed into a ‘core’ + ‘valence particle’ structure 4.3), such that the opposing nucleus interacts with the two components with distinct potentials acting on distinct centres-of-mass. The valence particle can be a single nucleon, as in the case of ^{17}O = 16O + n, or it can be a cluster of nucleons, as in ^{6}Li =α + ^{2}H, or ^{7}Li =α + ^{3}H. In all these cases, there arise inelastic potentials which can re-orient the ground state of the composite nucleus, or can excite the valence particle into higher-energy eigenstates.

Finally, transfer interactions (4.4) arise when the reaction brings about the transfer of a valence particle from one nucleus into a bound state around the other. As the incoming and outgoing projectiles have different centres-of-mass, with the targets likewise, the correct treatment of transfer interactions requires taking into account the effects of recoil and of the finite ranges of the binding potentials. These result in the coupling form factors becoming non-local, so that they must be specified by the two-dimensional kernel functions V_{α:α′}(R_{κ}, R_{κ′}) in equation (30). They also require that the coupled equations be solved by iteration, as will be discussed in section 5. If the effects of recoil are neglected, the ‘no-recoil’ (NR) approximation is obtained, but in general[72] this is inaccurate in ways which are difficult to predict. For that reason the NR approximation is not included in the present code. For many light-ion reactions, however, another ‘zero-range’ approximation is available, and this does remove many of the finite-range requirements. Alternatively, a first-order correction for the finite-range effects may be estimated, to give the ‘local energy approximation’. These two special cases are discussed at the end of the section.

4.1 Matrix Elements of Tensor Forces

This section presents the matrix elements for spin-orbit forces and a variety of tensor interactions. The radial form factors V_{Q}(R) which multiply these matrix elements are not specified, since these are usually determined by a fitting procedure in an optical-model search code, and a wide variety of parameterised forms have been used.

We shall use the |(LJ_{p})J_{1}, J_{t}; J_{T} ⟩ representation for the order of coupling the spins, as in equation (24).

4.1.1 Spin-orbit Interactions

For the projectile spin-orbit force L · J_{p}

⟨(LJ_{p})J_{1}, J_{t}; J_{T} |L · J_{p}|(L′J_{p})J′_{1}, J_{t}; J_{T} ⟩
= δ_{LL′}δ_{J_{1}J′_{1}} \frac{1}{2}[J_{1}(J_{1} + 1) − L(L + 1) − J_{p}(J_{p} + 1)]  (71)

This convention amounts to a 2l · s spin-orbit force, rather than one based on l · σ. These are the same for nucleons and spin \frac{1}{2} nuclei, but it means, for example, that the spin-orbit strengths for deuterons and 7Li will have to be decreased as they have s = 1 and 3/2 respectively.

For the target spin-orbit interaction L · J_{t}, we first transform

|(LJ_{p})J_{1}, J_{t}; J_{T} ⟩ = (−1)^{J_{1}−L−J_{p}}|(J_{p}L)J_{1}, J_{t}; J_{T} ⟩
= (−1)^{J_{1}−L−J_{p}} ∑_{J_{2}} |J_{p}, (LJ_{t})J_{2}; J_{T} ⟩Ĵ_{1}Ĵ_{2}W (J_{p}LJ_{T} J_{t}; J_{1}J_{2})  (72)

so

⟨(LJ_{p})J_{1}, J_{t}; J_{T} |L · J_{t}|(L′J_{p})J′_{1}, J_{t}; J_{T} ⟩

= (−1)^{J_{1}−J′_{1}+L′−L}Ĵ_{1}Ĵ′_{1} ∑_{J_{2}}(2J_{2} + 1)W (J_{p}LJ_{T} J_{t}; J_{1}J_{2})W (J_{p}LJ_{T} J_{t}; J′_{1}J_{2})

δ_{LL′} \frac{1}{2}[J_{2}(J_{2} + 1) − L(L + 1) − J_{t}(J_{t} + 1)]

4.1.2 Second-rank Tensor Forces

We use the notations of ref. [18]:

⟨S∥S_{2}∥S⟩ = \frac{1}{\sqrt{6}} \frac{3K^{2} − S(S + 1)}{⟨SK20|SK⟩} for any |K| ≤ S,  (73)

and

⟨L′∥R_{2}∥L⟩ = \sqrt{\frac{2}{3}} \frac{L̂}{L̂′} ⟨L 0 2 0 |L′ 0 ⟩  (74)

for the reduced matrix elements of the second-rank spin and radial tensors respectively. With the projectile T_{r} tensor force R_{2} · S_{2}(J_{p}J_{p}), the coupling interactions are

⟨(LJ_{p})J_{1}, J_{t}; J_{T} |R_{2} · S_{2}(J_{p}J_{p})|(L′J_{p})J′_{1}, J_{t}; J_{T} ⟩

= δ_{J_{1}J′_{1}} L̂Ĵ_{p}(−1)^{J_{1}−L−J_{p}}W (LL′J_{p}J_{p}; 2J_{1})⟨L∥R_{2}∥L′⟩⟨J_{p}∥S_{2}∥J_{p}⟩  (75)

For the target T_{r} tensor force R_{2} · S_{2}(J_{t}J_{t}) the coupling interactions are

⟨(LJ_{p})J_{1}, J_{t}; J_{T} |R_{2} · S_{2}(J_{t}J_{t})|(L′J_{p})J′_{1}, J_{t}; J_{T} ⟩

= (−1)^{J_{1}−J′_{1}+L′−L} ∑_{J_{2}} Ĵ_{1}Ĵ′_{1}(2J_{2} + 1)W (J_{p}LJ_{T} J_{t}; J_{1}J_{2})W (J_{p}L′J_{T} J_{t}; J′_{1}J_{2})

×L̂Ĵ_{t}(−1)^{J_{2}−L−J_{t}}W (LL′J_{t}J_{t}; 2J_{2})⟨L∥R_{2}∥L′⟩⟨J_{t}∥S_{2}∥J_{t}⟩  (76)

For the combined target-projectile T_{r} tensor force R_{2} · S_{2}(J_{p}J_{t}) the coupling interactions are

⟨(LJ_{p})J_{1}, J_{t}; J_{T} |R_{2} · S_{2}(J_{p}J_{t})|(L′J_{p})J′_{1}, J_{t}; J_{T} ⟩

= ∑_{SS′} Ĵ_{1}Ĵ′_{1}ŜŜ′W (LJ_{p}J_{T} J_{t}; J_{1}S)W (L′J_{p}J_{T} J_{t}; J′_{1}S′)

×L̂Ŝ (−1)^{J_{T} −L−S′}W (LL′SS′; 2J_{T})⟨L∥R_{2}∥L′⟩⟨(J_{p}J_{t})S∥S_{2}∥(J_{p}J_{t})S′⟩  (77)

where the second-rank reduced matrix element is

⟨(J_{p}J_{t})S∥S_{2}∥(J_{p}J_{t})S′⟩ = Ŝ′ˆ2Ĵ_{p}Ĵ_{t} ( \begin{matrix}S & J_{p} & J_{t} \\ S′ & J_{p} & J_{t} \\ 2 & 1 & 1 \end{matrix} ) \sqrt{J_{p}(J_{p} + 1)}\sqrt{J_{t}(J_{t} + 1)}

4.2 Inelastic Excitations

4.2.1 Nuclear Rotational Model

Consider a deformed nucleus with deformation lengths δ_{λ}. The effect of these deformations can be expressed as a change in the radius at which we evaluate the optical potentials, the change depending on the relative orientations of the radius vector to the intrinsic orientation of the nucleus. Deformation lengths are used to specify the these changes, rather than fractional deformations β_{λ}, to remove a dependence on the ‘average potential radius’ R_{U}. This is desirable because often the real and imaginary parts of the potential have different radii, and it is not clear which is to be used. It also removes a dependence on exactly how the ‘average radius’ of a potential is to be defined.

When U(R) is the potential shape to be deformed, the coupling interaction is

V(ξ, R) = U(R − δ(R̂, ξ))  (78)

where the ‘shift function’ has the multipole expansion

δ(R̂′) = ∑_{λ≠0} δ_{λ}Y_{λ}^{0}(R̂′)  (79)

(R̂′ is the vector R̂ in the body-centred frame of coordinates defined by ξ). Transforming to the spacefixed frame of reference, and projecting onto the spherical harmonics, the multipole expansion becomes

V(ξ, R) =  ∑_{λμ} V_{λ}(R)D_{μ0}^{λ}Y_{λ}^{μ}(R̂)  (80)
where V_{λ}(R) = \frac{1}{2} ∫_{−1}^{+1} U( r(R, cos θ) )Y_{λ}^{μ}(θ, 0) d(cos θ)  (81)
and r(R, u) =  R − \sqrt{\frac{2λ+1}{4π}}P_{λ}(u)δ_{λ} + ϵ  (82)
with ϵ =  ∑_{λ} δ_{λ}^{2}/(4πR_{U})  (83)

The correction ϵ is designed ([45]) to ensure that the volume integral of the monopole potential V_{0}(R) is the same as that of U(R), and is correct to second order in the {δ_{λ}}}.

When the {δ_{λ}}} are small, the above multipole functions are simply the first derivatives of the U(R) function:

V_{λ}(R) = −\frac{δ_{λ}}{\sqrt{4π}} \frac{dU(R)}{dR},  (84)

with the same shape for all multipoles λ > 0.

4.2.2 Coulomb Deformations

The deformations of the Coulomb potential can also be defined by the δ_{λ}, but unfortunately an average potential radius is again introduced. The dependence on models for average radii can be reduced by defining the Coulomb deformations in terms of a reduced matrix element such as that of Brink and Satchler [19], or that of Alder and Winther [20]. For the present purposes we adopt that of Alder and Winther, as it is hermitian upon interchanging the forward and reverse directions. We include, however, a simple phase factor to keep it real-valued. The new deformation parameter is called M(Eλ) and has units of e.fm^{λ}. In terms of the Alder and Winther reduced matrix element it is

M(Eλ) = i^{I−I′+|I−I′|} × ⟨I′∥Eλ∥I⟩  (85)

and is directly related to the observable electro-magnetic transition rate without any model-dependent parameters entering (except a sign):

M(Eλ) = ±\sqrt{(2I + 1)B(Eλ, I → I′)}.  (86)

A model dependent radius parameter R_{c} only enters in the relation to the deformation lengths of the rotational model:

M(Eλ, I → I′) = \frac{3Zδ_{λ}R_{c}^{λ−1}}{4π} i^{I−I′+|I−I′|} \sqrt{2I + 1} ⟨IKλ0|I′K⟩  (87)

for transitions from a state of spin I to one of spin I′ in a rotational band of projection K in a nucleus of charge Z. Within K=0 bands, M(Eλ, 0 → I) = M(Eλ, I → 0) have the same sign as δ_{λ}.

The only disadvantage of using reduced matrix elements as input parameters in this way is that the transitions in a rotational band do not all have the same matrix elements M(Eλ, I → I′), even when the deformation length is constant.

The radial form factors for Coulomb inelastic processes may be simply derived from the multipole expansion of |r − r′|^{−1}, giving

V_{λ}^{c}(R) = M(Eλ) \frac{\sqrt{4π}e^{2}}{2λ + 1} \{ \begin{matrix}R^{λ}/R_{c}^{2λ+1} (R ≤ R_{c}) \\ 1/R^{λ+1} (R > R_{c}) \end{matrix}  (88)

remembering that a factor δ_{λ}R_{c}^{λ−1} = β_{λ}R_{c}^{λ} is already included in the matrix element of equation (87) which appears in this form factor. This form factor is to be multiplied by the angular momentum coupling coefficients of the next section, and also by the charge of the opposing nucleus.

4.2.3 Angular Momentum Coupling Coefficients

The basic rotational coupling coefficient, with V_{λ} given by equation (70), is

X_{LI:L′I′}^{Jλ}(R) = ⟨LI; J|V_{λ}|L′I′; J⟩  (89)

The Coulomb form factors V_{λ}^{c}(R) have coupling coefficients

X_{LI:L′I′}^{Jλ}(R) = L̂L̂′(−1)^{J−I′−L+L′}W (LL′II′; λJ)⟨L0L′0|λ0⟩ V_{λ}^{c}(R)  (90)

whereas the nuclear form factors V_{λ}(R) defined for a rotational band with projection K have coupling coefficients

X_{LI:L′I′}^{Jλ}(R) = L̂L̂′(−1)^{J−I′−L+L′}W (LL′II′; λJ)⟨L0L′0|λ0⟩ V_{λ}(R)
Î′⟨I′Kλ0|IK⟩.  (91)

For projectile inelastic excitation, this coupling coefficient may be used directly as

⟨(LJ_{p})J, J_{t}; J_{T} |V_{λ}|(L′J′_{p})J′, J′_{t}; J_{T} ⟩ = δ_{J_{t},J′_{t}} δ_{J,J′} X_{LJ_{p}:L′J′_{p}}^{Jλ}(R)  (92)

whereas for target excitations,

⟨(LJ_{p})J, J_{t}; J_{T} |V_{λ}|(L′J′_{p})J′, J′_{t}; J_{T} ⟩ = δ_{J_{p},J′_{p}} (−1)^{J−J′−L+L′} ĴĴ′
× ∑_{J_{2}}(2J_{2} + 1)W (J_{p}LJ_{T} J_{t}; JJ_{2})W (J_{p}L′J_{T} J′_{t}; J′J_{2})X_{LJ_{t}:L′J′_{t}}^{J_{2}λ}(R)  (93)

4.3 Single Particle Excitations

When a nucleus consists of a single particle outside a core, the state of the particle can be disturbed by the interaction with1 another nucleus, as the force of that nucleus can act differentially on the particle and the core. If V_{cc}(R_{c}) and V_{p}(r′) are the interactions of the second nucleus with the core and particle respectively, then the excitation coupling from state |(ℓ′L′)λ⟩ to state |(ℓL)Λ⟩ is given by the singlefolding expression

X_{ℓL:ℓ′L′}^{Λ}(R) = ⟨(ℓL)Λ|V_{cc}(R_{c}) + V_{p}(r′) − U_{opt}(R)|(ℓ′L′)Λ⟩  (94)

where U_{opt}(R) is the optical potential already defined for these channels. This optical potential is subtracted to avoid double counting of either the Coulomb or the nuclear potentials, rather than disabling the potentials which have already been defined. This means that the ‘monopole’ potential V_{0}(R, r) (to be constructed) will have no long-range Coulomb component, and will not disturb the matching of the wave functions to the asymptotic Coulomb functions. It also means that if a nuclear well has already been defined, the new monopole form factor will be simply the difference between this well and that desired well calculated from the folding procedure.

If the potentials V_{cc}(R_{c})and U_{opt}(R) contain only scalar components, then the R- and r- dependent Legendre multipole potentials can be formed as

V_{K}(R, r) = \frac{1}{2} ∫_{−1}^{+1} [V_{cc}(R_{c}) + V_{p}(r′) − U_{opt}(R)] .P_{K}(u)du  (95)

where

K = the multipole moment,
u = r̂ · R̂ is the cosine of the angle between r and R,
r = aR + br is the particle-core vector,
and R_{c} = pR + qr is the core-nucleus vector.

The coupling form factor between states u_{ℓ′}(r) and u_{ℓ}(r) is then

X_{ℓL:ℓ′L′}^{Λ}(R) = \frac{1}{2} ∑_{K} ∫_{0}^{R_{m}} u_{ℓ}(r)^{∗}V_{K}(R, r)u_{ℓ′}(r)dr(−1)^{Λ+K} ℓ̂L̂ℓ̂′L̂′

× (2K + 1)W (ℓℓ′LL′; KΛ) ( \begin{matrix}K & ℓ & ℓ′ \\ 0 & 0 & 0 \end{matrix} ) ( \begin{matrix}K & L & L′ \\ 0 & 0 & 0 \end{matrix} )  (96)

4.3.1 Projectile Single-Particle Mechanisms

If the projectile has the particle - core composition, then the coupling interaction is

V_{α:α′}^{J_{T}}(R) = ⟨(LJ_{p})J, J_{t}; J_{T} |V|(L′J′_{p})J, J_{t}; J_{T} ⟩  (97)

where the initial (primed) and final (unprimed) states are

ϕ_{J′_{p}}(ξ_{p}, r) = ∑_{ℓ′sj′} A_{ℓ′sj′}^{j′I_{p}J′_{p}}|(ℓ′s)j′, I_{p}; J′_{p}⟩ and ϕ_{J_{p}}(ξ_{p}, r) = ∑_{ℓsj} A_{ℓsj}^{jI_{p}J_{p}}|(ℓs)j, I_{p}; J_{p}⟩,  (98)

respectively, and I_{p} is the (fixed) spin of the core. Then

V_{α:α′}^{J_{T}}(R) = ∑ ; \begin{matrix}F ΛI_{p} \\ jj′ℓℓ′ \end{matrix} ĵĵ′(2F + 1)Ĵ_{p}Ĵ′_{p}(2Λ + 1)W (ℓsJ_{p}I_{p}; jF)W (ℓ′sJ′_{p}I_{p}; j′F)

×A_{ℓ′sj′}^{j′I_{p}J′_{p}} A_{ℓsj}^{jI_{p}J_{p}} W (LℓJF; ΛJ_{p})W (L′ℓ′JF; ΛJ′_{p}) × X_{ℓL:ℓ′L′}^{Λ}(R)  (99)

4.3.2 Target Single-Particle Mechanisms

If the target has the particle - core composition, then the coupling interaction is

V_{α:α′}^{J_{T}}(R) = ⟨(LJ_{p})J, J′_{t}; J_{T} |V|(L′J_{p})J, J_{t}; J_{T} ⟩  (100)

where the initial (primed) and final (unprimed) states are

ϕ_{J′_{t}}(ξ_{t}, r) = ∑_{ℓ′sj′} A_{ℓ′sj′}^{j′I_{t}J′_{t}}|(ℓ′s)j′, I_{t}; J′_{t}⟩ and ϕ_{J_{t}}(ξ_{t}, r) = ∑_{ℓsj} A_{ℓsj}^{jI_{t}J_{t}}|(ℓs)j, I_{t}; J_{t}⟩,  (101)

respectively, and I_{t} is the (fixed) spin of the core in the target. Then

V_{α:α′}^{J_{T}}(R) = ∑_{I_{t}jj′ℓℓ′} A_{ℓ′sj′}^{j′I_{t}J′_{t}} A_{ℓsj}^{jI_{t}J_{t}}

× ∑_{J_{a}}(2J_{a} + 1)Ĵ′_{t}Ĵ_{t} W (JjJ_{T} I_{t}; J_{a}J_{t}) W (J′j′J_{T} I_{t}; J_{a}J′_{t})

×∑_{Λs_{a}} \{ \begin{matrix}L & ℓ & Λ \\ J_{p} & s & s_{a} \\ J & j & J_{a} \end{matrix} \} \{ \begin{matrix}L′ & ℓ′ & Λ \\ J_{p} & s & s_{a} \\ J′ & j′ & J_{a} \end{matrix} \} X_{ℓL:ℓ′L′}^{Λ}(R)  (102)

4.4 Particle Transfers

4.4.1 Finite Range Transfers

To calculate the coupling term that arises when a particle is transferred, for example from a target bound state to being bound in the projectile, we need to evaluate source terms of the form

S_{α}(R) = ∫_{0}^{∞}⟨(LJ_{p})J, J_{t}; J_{T} |V|(L′J′_{p})J′, J′_{t}; J_{T} ⟩ f_{(L′J′_{p})J′,J′_{t}}^{J_{T}}(R′)dR′  (103)

where the initial (primed) state has a composite target with internal coordinates ξ′_{t} ≡ {ξ_{t}, r′} : ϕ_{J′_{t}}(ξ_{t}, r′) = |(ℓ′s)j′, J_{t}; J′_{t}⟩ and the final (unprimed) state has a composite projectile with internal coordinates ξ_{p} ≡ {ξ_{p′}, r} : ϕ_{J_{p}}(ξ′_{p}, r) = |(ℓs)j, J′_{p}; J_{p}⟩.

The V is the interaction potential, of which the prior form is

V = V_{ℓsj}(r) + U_{cc}(R_{c}) − U_{α′}(R′)  (104)

and the post form is

V = V_{ℓ′sj′}(r′) + U_{cc}(R_{c}) − U_{α}(R)  (105)

where V_{β}(r) is the potential which binds φ_{β}(r), U_{α}(R) are the optical potentials, and U_{cc}(R_{c}) is the ‘core-core’ potential, here between the p′ and the t nuclei. The V_{β} will be real, but the U_{α} and U_{cc} will typically have both real and imaginary components.

This source function S_{α}(R) evaluates a non-local integral operator, as it operates on the function f_{α′}(R′) to produce a function of R. This section therefore derives the non-local kernel V_{α,α′}(R, R′) so that the source term, which initially involves a five dimensional integral over r and R̂, may be calculated by means of a one-dimensional integral over R′:

S_{α}(R) = ∫_{0}^{R_{m}} V_{α,α′}(R, R′)f_{α′}(R′)dR′.  (106)

Note that when the initial and final single-particle states are real, then the kernel function is symmetric

V_{α,α′}(R, R′) = V_{α′,α}(R′, R),  (107)

whereas if the states are unbound and complex-valued, then the kernel function is hermitian provided the interaction potential V is real. If the particle states and the interaction potential are complex, then both the forward and reverse kernels must be each calculated independently.

When the potential V contains only scalar potentials, the kernel calculation can be reduced to the problem of finding X_{ℓL:ℓ′L′}^{Λ}(R, R′) such that, given

⟨(LJ_{p})J, J_{t}; J_{T} |V|(L′J′_{p})J′, J′_{t}; J_{T} ⟩ = ∑_{ΛF}(−1)^{s+J′_{p}−F} ĴĴ′_{t}ĵF̂Ĵ_{p}Λ̂ \{ \begin{matrix}L′ & J′_{p} & J′ \\ ℓ′ & s′ & j′ \\ Λ & F & J \end{matrix} \}

×W (J_{t}j′J_{T} J′; J′_{t}J)W (lsJ_{p}J′_{p}; jF)W (LℓJF; ΛJ_{p})⟨ℓL; Λ|V|ℓ′L′; Λ⟩,  (108)

the integral operator ⟨ℓL; Λ|V|ℓ′L′; Λ⟩ has the kernel function X_{ℓL:ℓ′L′}^{Λ}(R, R′). Note that the F summation may be performed in an inner loop that does not evaluate the kernel function.

Now the r and r′ are linear combinations of the channel vectors R and R′: r = aR + bR′ and r′ = a′R + b′R′ where, when φ_{ℓ}(r) is the projectile bound state,

a = ν_{t}ω, b = −ω, a′ = ω, b′ = −ν_{p}ω,  (109)

with ν_{p} ≡ A_{κ′p}/A_{κp} , ν_{t} ≡ A_{κt}/A_{κ′t} , and ω = (1 − ν_{p}ν_{t})^{−1} . When φ_{ℓ}(r) is the target bound state

a = −ν_{p}ω, b = ω, a′ = −ω, b′ = ν_{t}ω,  (110)

with ν_{p} ≡ A_{κp}/A_{κ′p} , ν_{t} ≡ A_{κ′t}/A_{κt} , and ω = (1 − ν_{p}ν_{t})^{−1} . The ‘core-core’ vector is always R_{c} = r′ − r = (a′ − a)R + (b′ − b)R′.

Thus the spherical harmonics Y_{ℓ}(r̂) and Y_{ℓ′}(r̂′) can be given in terms of the spherical harmonics Y_{n}(R̂) and Y_{n′}(R̂′) by means of the Moshinsky [71] solid-harmonic expansion (see also refs. [21] and [46]

Y_{ℓ}^{m}(r̂) = \sqrt{4π} ∑_{nλ} c(ℓ, n)\frac{(aR)^{ℓ−n}(bR′)^{n}}{r^{ℓ}} Y_{ℓ−n}^{m−λ}(R̂)Y_{n}^{λ}(R̂′)⟨ℓ − nm − λnλ|ℓm⟩  (111)

where

c(ℓ, n) = \sqrt{\frac{1}{2n + 1} ( \begin{matrix}2ℓ + 1 \\ 2n \end{matrix} )},

with ( \begin{matrix}x \\ y \end{matrix} ) the binomial coefficient (Appendix A).

We now perform the Legendre expansion

V\frac{u_{ℓsj}(r)}{r^{ℓ+1}} \frac{u_{ℓ′sj′}(r′)}{r′^{ℓ′+1}} = ∑_{T}(2T + 1)q_{ℓ,ℓ′}^{T}(R, R′)P_{T} (u)  (112)

where the Legendre polynomials P_{T} (u) are functions of u, the cosine of the angle between R and R′, by using r = (a^{2}R^{2} + b^{2}R′^{2} + 2abRR′u)^{1/2} (with r′ analogously) in the numerical quadrature of the integral

q_{ℓ,ℓ′}^{T}(R, R′) = \frac{1}{2} ∫_{−1}^{+1} V\frac{u_{ℓsj}(r)}{r^{ℓ+1}} \frac{u_{ℓ′sj′}(r′)}{r′^{ℓ′+1}} P_{T} (u)du  (113)

The quadrature methods used here, and the accuracy attained, are discussed in section 5.3.

Using the Legendre expansion, the radial kernel function

X_{ℓL:ℓ′L′}^{Λ}(R, R′) = \frac{|b|^{3}}{2} ∑_{nn′} c(ℓn)c(ℓ′n′)RR′(aR)^{ℓ−n}(bR′)^{n}(a′R)^{ℓ′−n′}(b′R′)^{n′}

× ∑_{T} q_{ℓ,ℓ′}^{T}(R, R′)(2T + 1)(−1)^{Λ+T +L+L′} ℓ̂ℓ̂′ ˆ (ℓ − n) ˆ (ℓ′ − n′)n̂n̂′L̂L̂′

× ∑_{KK′}(2K + 1)(2K′ + 1) ( \begin{matrix}ℓ − n & n′ & K \\ 0 & 0 & 0 \end{matrix} ) ( \begin{matrix}ℓ′ − n′ & n & K′ \\ 0 & 0 & 0 \end{matrix} ) ( \begin{matrix}K & L & T \\ 0 & 0 & 0 \end{matrix} ) ( \begin{matrix}K′ & L′ & T \\ 0 & 0 & 0 \end{matrix} )

× ∑_{Q}(2Q + 1)W (ℓLℓ′L′; ΛQ)W (KLK′L′; T Q) ( \begin{matrix}ℓ′ & Q & ℓ \\ n′ & K & ℓ − n \\ ℓ′ − n′ & K′ & n \end{matrix} )  (114)

These formulae can also be used with V ≡ 1 to calculate the kernel functions K_{ℓL:ℓ′L′}^{Λ}(R, R′) for the wave function overlap operator K_{ij} ≡ ⟨ϕ_{i}|ϕ_{j}⟩ needed in evaluating the non-orthogonality terms of section 2.3.

One disadvantage of this method of calculating the two-dimensional radial kernels X_{ℓL:ℓ′L′}^{Λ}(R, R′) is that in the process of transforming the solid harmonics of r and r′ into those of R and R′, there appears summations containing high powers of the coefficients a, b, a′ and b′ These products will become larger than unity by several orders of magnitude, will the summed result is typically of the order of unity. This means that the summations involve large cancellations, and as the degree of cancellation gets worse for large ℓ and ℓ′, the cancellation places a limit on the maximum value ℓ + ℓ′ of the transferred angular momentum.

One way of circumventing this loss of accuracy is that proposed by Tamura and Udagawa [47], whereby solid harmonics are avoided in favour of a suitable choice of axes to render it practical to calculate mdependent form factors directly. If the z axis is not (as usual) parallel to the incident momentum, but set parallel to R, and the x′ axis set in the plane determined by R and R′, then the r and r′ vectors are also in this plane. The radial kernels may then be calculated as a sum of m-dependent integrals over cos θ = R̂ · R̂′, as before the cosine of the angle between R and R′. Although there are hence a larger number of radial integrals to be performed, there are no large cancellations between the separate terms, and there is no limit on the size of the transferred angular momentum.

A third method [22] of calculating the transfer form factors is that involving expanding the initial and final channel wave functions in terms of spherical Bessel functions:

f_{α}(R) = ∑_{n=1}^{N(L)} a_{α}(K_{n})R j_{L}(K_{n}R).  (115)

Using then the Fourier transform of the bound state wave functions u_{ℓ}(r) and u_{ℓ′}(r′), a transfer Tmatrix element may be written as a sum of a set of one-dimensional integrals over a momentum variable. Efficient codes [23] have been written for CCBA calculations of transfers induced by light ions (up to masses ∼ 10 to 15 amu).

This plane-wave expansion method has however several disadvantages when it comes to solving problems with coupled reaction channels. If transfers are to be calculated at each iteration of the coupled equations, then the expansion (115) has to be recalculated at each step. Another difficulty is that the method is not suited to heavy-ion induced transfers, as the large degree of absorption inside the nuclei in these cases requires a large number of momentum basis states K_{n} to be represented accurately. The plane-wave expansion becomes uneconomical, and sometimes the determination of the a_{α}(K_{n}) coefficients becomes numerically ill-conditioned

We will see in section 5.3.1, however, that if the cancellation which occurs in the first method is monitored, and steps taken to keep it to a minimum, a workable code [34] results which can produce accurate results for L-transfers up to around 6.

4.4.2 Zero Range Transfers

When the projectile wave functions φ_{ℓ}(r) are all s-states (ℓ = 0 and the interaction potential is of zerorange (Vφ(r) ∼ D_{0}δ(r) ), then the form factor X_{ℓL:ℓ′L′}^{Λ}(R, R′) of equation (114) can be simplified to

X_{0L:ℓ′L′}^{L}(R, R′) = D_{0} \frac{(−1)^{L′−ℓ′}}{L̂} \begin{matrix}ℓ̂′ & L̂ & L̂′ \\ √ & \overline{4π} \end{matrix} ( \begin{matrix}ℓ′ & L & L′ \\ 0 & 0 & 0 \end{matrix} ) \frac{1}{R}u_{ℓ′sj′}(R) \frac{b^{2}}{a} δ(aR + bR′). (116)

This can be made local by defining a new step size h′ = −ah/b ≡ ν_{t}h in the stripping channel α′.

4.4.3 Local Energy Approximation

If the interaction potential is of small range, though not zero, and the projectile still contains only sstates, then a first-order correction may be made to the above form factor. This correction will depend on the rate of oscillation of the source wave function f_{(L′J′_{p}),J′,J′_{t}}^{J_{T}}(R′) within a ‘finite-range effective radius’ ρ. The rate of oscillation is estimated from the local energy in the entrance and exit channels, and the result [24] is to replace u_{ℓ′sj′}(R) in the previous section by

u_{ℓ′sj′}(R) → u_{ℓ′sj′}(R) [1 + ρ^{2} \frac{2μ_{α}^{(p)}}{h̄^{2}} (U_{α′}(R) + V_{ℓ′sj′}(R) − U_{α}(R) + ϵ_{α})]  (117)

where the U(R) are the optical potentials, with V_{ℓ′sj′}(r) the single-particle binding potential in the target. The μ_{α}^{(p)} is the reduced mass of the particle in the projectile, and ϵ_{α} its binding energy.

At sub-Coulomb incident energies [25], the details of the nuclear potentials in equation (117) become invisible, and as the longer-ranged Coulomb potentials cancel by charge conservation, the form factor can be simplified to

D_{0}u_{ℓ′sj′}(R) → u_{ℓ′sj′}(R)D_{0} [1 + ρ^{2} \frac{2μ_{α}^{(p)}}{h̄^{2}} ϵ_{α}] = u_{ℓ′sj′}(R)D  (118)

where

D = D_{0} [1 + (ρk_{α}^{(p)})^{2}]  (119)

is the effective zero-range coupling constant for sub-Coulomb transfers.

The parameters D_{0} and D can be derived from the details of the projectile bound state φ_{0ss}(r). The zero-range constant D_{0} may be defined as

D_{0} = \sqrt{4π} ∫_{0}^{∞} rV (r)u_{0ss}(r)dr.  (120)

The parameter D, on the other hand, reflects the asymptotic strength of the wave function u_{0ss}(r) as r → ∞, as it is the magnitude of this tail which is important in sub-Coulomb reactions:

u_{0ss}(r) =_{r→∞} \frac{2μ_{α}^{(p)}}{h̄^{2}} \frac{1}{\sqrt{4π}} De^{−k_{α}^{(p)}r}.  (121)

It may be also found, using Schrödinger’s equation, from the integral

D = \sqrt{4π} ∫_{0}^{∞} \frac{sinh(k_{α}^{(p)}r)}{k_{α}^{(p)}} V (r)u_{0ss}(r)dr.  (122)

From this equation we can see that as the range of the potential becomes smaller, D approaches D_{0}. The ‘finite-range effective radius’ ρ of equation (119) is thus some measure of the mean radius of the potential V (r).

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_{α:α′}^{Γ}(R_{κ}) and non-local kernels V_{α:α′}(R_{κ}, R_{κ′}). 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_{i,k}(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^{3}), and if there are only local couplings. For in that case the independent solutions can be found in a single outward ‘sweep’ of m^{2} radial functions. Nonlocal 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_{κ}(R_{κ}) usually have a significant imaginary component for small R_{κ}, and hence damp the solutions f_{α}(R_{κ}) in this region. This enables lower radial cutoffs to be used for R_{κ} near zero, with little loss of accuracy.

3. The bound states u_{ℓsjI}(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 ≡ cosθ ≡ R̂ · R̂′ decreases from unity.

5.1 Integration of the Radial Equations

If the non-local interactions V_{α,α′}(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_{α,α′}^{Γ}(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 lowlying 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

\frac{d^{2}}{dR^{2}} f_{i}(R) = ∑_{j=1}^{b} C_{ij}(R)f_{j}(R) + S_{i}(R) for i = 1 ⋯ b,  (123)

and \frac{d^{2}}{dR^{2}} f_{i}(R) = C_{ii}(R)f_{i}(R) + S_{i}(R) for i = b + 1 ⋯ 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) = ∑_{j=j_{min}}^{m} C_{ij}(R)f_{j}^{(n−1)}(R)  (125)

where j_{min} = b + 1 if i ≤ b and j_{min} = 1 if i > b.

IT

1 g_{1,0}  g_{2,0}  ..  g_{m,0} inhomogeneous solns.
2  (not used)  g_{b+1,b+1}  ..  g_{m,m} uncoupled regular solns.
3 g_{1,1}  g_{2,1}  ..
4 g_{1,2}  g_{2,2}  ..  coupled regular solutions
↓  ..  ..  ..  (equations 1 → b)
b + 2 g_{1,b}  g_{2,b}  .. g_{b,b}

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.

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 ⋯ m) which is made independent of the other sets by having a distinctive starting value

g_{i,k}(R_{min} − h) = 0, g_{i,k}(R_{min}) = \frac{1}{(2L_{i} + 1)!!} (K_{i}R_{min})^{L_{i}+1}δ_{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 ⋯ 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) = ∑_{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_{0} is always unity, to match the source terms correctly. The S-matrix elements are a by-product of the linear matching equations (32).

Note that the independent solutions g_{i,k}(R) for k ≥ 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 ω, 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 ≠ 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_{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_{min} ≥ 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)_{m}in for the off-diagonal coupling terms only. This allows the absorption in the diagonal potentials to be effective at all radii outside R_{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 ψ_{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)} 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_{α} in the relevant channel. This requirement is less severe with light-ion reactions, where the wave numbers are typically ≤ 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}, ⋯ of S-matrix elements that result from iterating the CRC equations can be regarded as the successive partial sums of the polynomial

f(λ) = S_{0} + (S_{1} − S_{0})λ + (S_{2} − S_{1})λ^{2} + ⋯  (129)

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

P_{[n,m]}(λ) = \frac{p_{0} + p_{1}λ + p_{2}λ^{2} + ⋯ + p_{n}λ^{n}}{1 + q_{1}λ + q_{2}λ^{2} + ⋯ + q_{m}λ^{m}}  (130)

which agrees with the f(λ) function in the region where the latter does converge, as tested by matching the coefficients in the polynomial expansion of P_{[n,m]}(λ) up to and including the coefficient of λ^{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 ϵ-algorithm [31], which is a method of evaluating the upper right half of the Padé table at λ=1 directly in terms of the original sequence S_{0}, S_{1}, · · ·.

5.2.3 Wynn’s epsilon Algorithm

Initialising ϵ_{0}^{(j)} = S_{j} and ϵ_{−1}^{(j)} = 0, we form an array using the relation ϵ^{(}j)_{k+1} = ϵ_{k−1}^{(j+1)} + (ϵ_{k}^{(j+1)} − ϵ_{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:

ϵ_{2k}^{(j)} = 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 ϵ_{2k}^{(0)} in the ϵ table.

When accelerating a vector S-matrix elements 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

x^{−1} = (x · x^{∗})^{−1}x^{∗}  (132)

where x^{∗} is the complex conjugate of x. Otherwise there will be problems when iterating (say) a twochannel system with alternating backwards and forwards coupling, because of zero divisors in the ϵ 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 ℓ and ℓ′, this places a limit on the maximum value ℓ+ℓ′ 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_{α}(R), the bound state wave functions u_{ℓsjI}(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 q_{ℓ,ℓ′}^{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 r and r′ in terms of R and R′ be accurate to full machine precision, as also the Racah algebra coefficients in equation (114). In fact, the channel wave functions f_{α}(R) and the bound state wave functions u_{ℓsjI}(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 (θ to 180^{∘}) 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 θ increases. For a given accuracy criterion, an estimate can then be made of an adequate upper limit for the θ integration at the next R value. Typically, the upper limits of θ decrease monotonically as R increases from 0 to the upper limit R_{m}. Because the integrand is largest for θ=0, the accuracy of the angular integration for small θ is improved by a change of variable from u to x as in ref.[55]:

θ = \frac{1}{4}(3x^{2} + 1)xθ_{max}  (133)

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

q_{ℓ,ℓ′}^{T}(R, R′) = \frac{1}{2} ∫_{0}^{1} V\frac{u_{ℓsj}(r)}{r^{ℓ+1}} \frac{u_{ℓ′sj′}(r′)}{r′^{ℓ′+1}} P_{T} (u) sin(θ)\frac{dθ}{dx}dx.  (134)

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

5.3.2 Radial Grids

The methods used to calculate, store and use the non-local form factors qu_{ℓ,ℓ′}^{T}(R, R′) (equation 113) and V_{α,α′}(R, R′) (equation 108) have to be efficient in a wide variety of reactions, from light-ion reactions such as^{3}He(^{3}H,^{4}He)^{2}H or^{16}O(^{20}Ne,^{24}Mg)^{12}C to heavy-ion reactions, such as nickel on tin one-nucleon transfers. In the former cases, the radial form factors V_{α,α′}(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 ∼ 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 δR ≡ R − R′ (especially for heavy ion reactions, as the Jacobian b^{3} in equation (114) becomes large), and varies only slowly with R (if δ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 δR and R, and then to use different interpolatory intervals h_{δ} and h_{R} in the two directions respectively. Then, when nuclear masses become large compared with the mass of the transferred particle, h_{δ} 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_{δ} 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 α-particle bound states).

The h_{δ}, 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 δR, which again will be large ( ∼ nuclear radii) for light ions, and small ( ∼ 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 q_{ℓ,ℓ′}^{T}(R, δR), averaging over R. and all partial waves T, ℓ, and ℓ′.

When h_{δ} 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_{δ} 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

I = ∑_{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) = ∑_{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 ∼ i ± 2). Then I can be evaluated directly in terms of the f_{m} :

I = ∑_{j} V (x_{j}) ∑_{m} a_{m}(x_{j})f_{m}  (137)

= ∑_{m} \overline{V}_{m}f_{m}  (138)

where

\overline{V}_{m} ≡ ∑_{j} V (x_{j})a_{m}(x_{j})  (139)

is a new effective form factor This means that when h_{δ} is a submultiple of h, we do not need to store a form factor at intervals of h_{δ}, 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{q}_{ℓ,ℓ′}^{T}(R, δR) and \overline{V}_{α,α′}(R, δR) need fewer grid points in the δR direction. Less arithmetic is needed to evaluate the source functions of equation (106), which change from

S_{α}(R) = ∫_{δR_{min}}^{δR_{max}} \overline{V}_{α,α′}(R, δR) f_{α′}(R − δR)d(δR).  (140)

to

S_{α}(ih_{R}) = h ∑_{j} \overline{V}_{α,α′}(ih_{R}, jh)f_{α′}((in_{R} − j)h) where n_{R} ≡ h_{R}/h  (141)

even when the original kernel functions vary rapidly as δ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 (ℓ, (s_{1}s_{2})S)j and ρ, 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 ρ value and summed in equation (106), or the separate ρ products can be summed as early as equation (113). Because the separate ρ 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

q_{L,L′}^{T}(R, R′) = \frac{1}{2} ∫_{−1}^{+1} r^{−L−1}r′^{−L′−1} [∑_{i} w_{i}Vu_{12}(r, ρ_{i}) u′_{12}(r′, ρ_{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.

References

[1] T. Tamura, Rev. Mod. Phys. 37 (1965) 679.

[2] T. Tamura, Phys. Reports 14C (1974) 59.

[3] N. Austern, Direct Nuclear Reaction Theories (Wiley, New York, 1970).

[4] M. Rhoades-Brown, M.H. Macfarlane and S.C. Pieper, Phys. Rev. C 21 (1980) 2417.

[5] L.D. Tolsma, Phys. Rev. C 20 (1979) 592, J. Comp. Phys. 17 (1975) 384.

[6] H. Feshbach, Annals of Phys. (N.Y.) 5 (1958) 357, 19 (1962) 287.

[7] T. Udagawa, H.H. Wolter and W.R. Coker, Phys. Rev. Letts. 31 (1973) 1507.

[8] B. Imanishi, M. Ichimura and M. Kawai, Phys. Lett. 52B (1974) 267.

[9] J.S. Lilley, B.R. Fulton, D.W. Banes, T.M. Cormier, I.J. Thompson, S. Landowne, and H.H. Wolter, Phys. Letts. 128B (1983) 153, and J.S. Lilley, M.A. Nagarajan, D.W. Banes, B.R. Fulton and I.J. Thompson, Nucl. Phys. A463 (1987) 710.

[10] W. Ogle, S. Wahlborn, R. Piepenbring and S. Frederiksson, Rev. Mod. Phys. 43 (1971) 424.

[11] J.S. Vaagen, B.S. Nilsson, J. Bang and R.M. Ibarra, Nucl. Phys. A319 (1979) 143; J.M. Bang, F.G. Gareev, W.T. Pinkston and J.S. Vaagen, Phys. Reports 125 (1985) 253.

[12] T.A. Brody and M. Moshinsky, Tables of Transformation Brackets., (Monografias del Instituto de Fisica, Mexico, 1960); D.H. Feng and T. Tamura, Comput. Phys. Commun. 10 (1975) 87.

[13] B. Bayman and A. Kallio, Phys. Rev. 156 (1967) 1121.

[14] M.V. Berry, Proc. Phys. Soc. 89 (1966) 479, J. Phys. B: Atom. Molec. Phys. 2 (1969) 381; R.C. Fuller, Phys. Rev. C 12 (1975) 1561.

[15] N. Rowley, J. Phys. A: Math. Gen. 11 (1978) 1545.

[16] M. Simonius, Lecture Notes in Physics., Ed D. Fick, 30 (Springer, 1974) 38.

[17] B. Imanishi, Nucl. Phys. A125 (1969) 33.

[18] G.R. Satchler, Nucl. Phys. 21 (1960) 116.

[19] D.M. Brink and G.R. Satchler, Angular Momentum., (Clarendon Press, Oxford, 1968).

[20] K. Alder and Å. Winther, Coulomb Excitation (Academic Press, New York, 1966) and Electromagnetic Excitation (North Holland, Amsterdam, 1975).

[21] N. Austern, R.M. Drisko, E.C. Halbert and G.R. Satchler, Phys. Rev. 133 (1964) B3.

[22] D. Robson and R.D. Koshel, Phys. Rev. C6 (1972) 1125.

[23] P. Nagel and R.D. Koshel, Phys. Rev. C13 (1976) 907, Comput. Phys. Commun. 15 (1978) 193.

[24] P.J.A. Buttle and L.J.B. Goldfarb, Proc. Phys. Soc. (London) 83 (1964) 701; for a review see G.R. Satchler, [33, §6.14.1].

[25] L.J.B. Goldfarb and E. Parry, Nucl. Phys. A116 (1968) 309.

[26] T. Tamura, Oak Ridge National Laboratory Report No. ORNL-4152 (1967).

[27] M.A. Melkanoff, T. Sawada and J. Raynal, Methods in Computational Physics, Vol. 6 (1966), Academic Press (New York).

[28] J. Bang, J.-J. Benayoun, C. Gignoux, and I.J. Thompson, Nucl. Phys. A405 (1983) 126.

[29] K. Alder, F. Rösel and R. Morf, Nucl. Phys. A284 (1977) 145.

[30] M. Ichimura, M. Igarashi, S. Landowne, C.H. Dasso, B.S. Nilsson, R.A. Broglia and Å. Winter, Phys. Letts. 67B 129.

[31] P.Wynn, Numer. Math. 8 (1966) 264; see also A. Genz, p 112 in ref.[48]

[32] P. Wynn, Math. Comput. 16 (1961) 23.

[33] G.R. Satchler, Direct Nuclear Reactions. (Clarendon, Oxford, 1983).

[34] I.J. Thompson, FRESCO users’ manual and code, available from the author.

[35] F. Rösel, J.X. Saladin and K. Alder, Comput. Phys. Commun. 8 (1974) 35.

[36] G.H. Rawitscher and C.H. Rasmussen, Comput. Phys. Commun. 11 (1976) 183.

[37] L.D. Tolsma, Phys. Rev. C 35 (1987) 177.

[38] J.D. Raynal, Phys. Rev. C23 (1983) 2571.

[39] G.H. Rawitscher, Phys. Rev. C9 (1974) 2210, C11 (1975) 1152, Nucl. Phys. A241 (1975) 365. M. Yahiro and M. Kamimura, Prog. Theor. Phys. 65 (1981) 2046, 65 (1981) 2051.

[40] Y. Sakuragi, M. Yahiro and M. Kamimura, Prog. Theor. Phys. Suppl. 89 (1986).

[41] S.R. Cotanch and C.M. Vincent, Phys. Rev. C14 (1976) 1739.

[42] I.J. Thompson, in: Light Ion Reaction Mechanism., proc. of RCNP Symposi

[43] P. Nagel and R.D. Koshel, Phys. Rev. C13 (1976) 907, Comput. Phys. Commun. 15 (1978) 193.um, Osaka, 1983, ed. H. Ogata.

[44] I.J. Thompson and A.R. Barnett, Comput. Phys. Commun. 36 (1985) 363, J. Comp. Phys. 64 (1986) 490.

[45] J.M. Bang and J.S. Vaagen, Z. Physik A 297 (1980) 223.

[46] T. Ohmura, B. Imanishi, M. Ichimura and M. Kawai, Prog. Theor. Phys. 41 (1969) 391, 43 (1970) 347, 44 (1970) 1242.

[47] T. Tamura, T. Udagawa, K.E. Wood and H. Amakawa, Comput. Phys. Commun. 18 (1979) 63. T. Tamura, T. Udagawa and M.C. Mermaz, Phys. Reports 65 (1980) 345.

[48] P.R. Graves-Morris, Padé Approximants, (Institute of Physics, London and Bristol, 1973)

[49] I.J. Thompson, Ph.D. report (University of Auckland, 1979, unpublished).

[50] M. Kawai, M. Kamimura, Y. Mito, K. Takesako, Prog. Theor. Phys. 59 (1978) 674 and 676. M. Kamimura, Prog. Theor. Phys. Suppl. No. 62 (1977) 236.

[51] J. Raynal, Computing As a Language of Physics, (IAEA, Vienna, 1972) 281.

[52] M. Rhoades-Brown, M.H. Macfarlane and S.C. Pieper, Phys. Rev. C 21 (1980) 2436.

[53] B. Buck, A.P. Stamp and P.E. Hodgson, Phil. Mag., 8 (1963) 1805.

[54] L.D. Tolsma and G.W. Veltkamp, Comput. Phys. Commun. 40 (1986) 233.

[55] P.J.A. Buttle, Comput. Phys. Commun. 14 (1978) 133.

[56] T. Tamura, Comput. Phys. Commun. 8 (1974) 349.

[57] S. Yoshida, Proc. Phys. Soc. (London) A69 (1956) 668.

[58] A.P. Stamp, Nucl. Phys. B3 (1966) 232.

[59] G.H. Rawitscher, Phys. Rev. 163 (1967) 1223.

G.H. Rawitscher and S.N. Mukherjee, Phys. Rev. 181 (1969) 1518.

[60] J.M. Bang and S. Wollesen, Phys. Lett. 33B (1970) 395.

[61] R.J. Ascuitto, N.K. Glendinning and B. Sørensen, Nucl. Phys. A170 (1971) 65.

[62] R.S. Mackintosh, Nucl. Phys. A209 (1973) 91.

[63] R. Schaeffer and G.F. Bertsch, Phys. Lett. 38B (1972) 159.

[64] M. Toyama, Phys. Lett. 38B (1972) 147.

[65] P.D. Kunz, Program CHUCK3 (University of Colorado, 1978, unpublished),

[66] Extended version of [65] by J.R. Comfort, (University of Pittsburgh, 1980, unpublished).

[67] A.J. Baltz, Phys. Rev. C25 (1982) 240.

[68] B. Imanishi and W. von-Oertzen, Phys. Lett. 87B (1979) 188.

[69] M. Toyama and M. Igarashi, 1972 and 1977, unpublished.

[70] N.M. Clarke, J. Phys. G: Nucl. Phys. 10 (1984) 1219.

[71] M. Moshinsky, Nucl. Phys. 13 (1959) 104.

[72] R.M. DeVries, Phys. Rev. C8 (1973) 951.

A Notation and Phase Conventions

Spherical Harmonics The phase convention used here is

Y_{L}^{m}(θ, ϕ) = \sqrt{\frac{2L + 1}{4π} \frac{(L − m)!}{(L + m)!}}(−1)^{m}e^{imϕ}P_{L}^{m}(cos θ)

for m ≥ 0, and Y_{L}^{−m} = (−1)^{m}Y_{L}^{m∗} to give negative m values. Angular Momentum Coupling Coefficients The notation ⟨ℓ_{1}m_{1}ℓ_{2}m_{2}|LM⟩ has been used for the Clebsch-Gordon coupling coefficient for coupling states ℓ_{1}m_{1} and ℓ_{2}m_{2} together to form LM. The

( \begin{matrix}a & b & c \\ α & β & γ \end{matrix} ) ≡ \frac{(−1)^{a−b−γ}}{ĉ} ⟨aαbβ|c − γ⟩

represents the Wigner 3-j symbol, and

x̂ ≡ \sqrt{2x + 1}.

The 9-j coupling coefficient is used in two forms related by

\{ \begin{matrix}a & b & c \\ d & e & f \\ g & h & i \end{matrix} \} ≡ ĉf̂ĝĥ ( \begin{matrix}a & b & c \\ d & e & f \\ g & h & i \end{matrix} ) .

The binomial coefficient is

( \begin{matrix}x \\ y \end{matrix} ) = \frac{x!}{y!(x − y)!}.

B Coupled Channels Codes in Nuclear Physics

There is a natural progression of complexity in the codes being considered here:

1. One-step (DWBA) codes Inelastic excitations Zero-range transfers (ZR) No-recoil transfers (NR)

2. Coupled channels (CC) codes with local form factors Inelastic excitations (CC) Zero-range transfers (ZR-CC, sometimes included in CCBA) No-recoil transfers (NR-CC, sometimes included in CCBA)

3. One-step (DWBA) codes for exact finite-range transfers (EFR-DWBA)

4. Coupled channels Born approximation (CCBA): a coupled set of channels followed by a finiterange transfer (sometimes called EFR-CCBA).

5. Two-step DWBA Zero-range transfers (2-step ZR-DWBA) No-recoil transfers (2-step NR-DWBA)

6. One-step DWBA codes for exact finite-range transfers (EFR-DWBA)

7. Two-step DWBA for exact finite-range transfers (2-step EFR-DWBA)

8. Coupled reaction channels (CRC), allowing finite-range transfers.

The following is a summary of the more widely known coupled channels codes (codes which can only perform one-step DWBA calculations have been excluded).

1. Yoshida [57]: Inelastic CC with δ-function interactions

2. Buck, Stamp and Hodgson [53] and Satchler (see ref.[53]): Inelastic CC

3. Tamura [1]: General purpose inelastic CC

4. Stamp [58]: ZR-CC

5. Rawitscher [59]: ZR-CC using iterated Green functions

6. Tamura and Low [56] and [47]: Saturn-Mars - NR-DWBA and EFR-DWBA

7. Ohmura et al. [46]: CRC for deuterons

8. Ascuitto et al. [61]: ZR-CCBA using source terms

9. Mackintosh [62]: ZR-CRC for deuterons and protons

10. Bang and Wollesen [60]: two-step ZR-DWBA.

11. Toyama [64]: two-step ZR-DWBA

12. Schaeffer and Bertsch [63]: two-step ZR-DWBA

13. Rösel et al. [35], extended by Rawitscher [36]: AROSA - CC for Coulomb excitations

14. Cotanch and Vincent [41]: CRC for deuterons

15. Raynal [51] and [38]: General purpose inelastic CC and ZR-CC

16. Kunz [65], and later Comfort [66]: CHUCK - General purpose CC (inelastic and ZR)

17. Nagel and Koshel [43]: OUKID - EFR-CCBA for light ions

18. Kawai [50]: CRC for deuterons

19. Baltz [67]: QUICC - Inelastic CC for heavy ions 20. Imanishi [68]: CRC for^{12}C+^{13}C reactions

21. Tolsma [5] and [37]: PIECANSOL - Inelastic CC for heavy ions

22. Toyama and Igarashi [69]: TWOSTP - 2-step ZR- and EFR-DWBA for light ions, and Igarashi: TWOFNR/PTFF - 2-step EFR-DWBA for sequential and simultaneous transfers of two nucleons for light ion reactions.

23. Thompson [49]: CRC for deuterons

24. MacFarlane, Pieper and Rhoades-Brown [52] PTOLEMY/1 - Inelastic DWBA and EFR-DWBA for heavy and light ions PTOLEMY/2 - General purpose inelastic CC for heavy ions

25. Kunz : CHORK - ZR-CC, and NR-CC for heavy ions

26. Thompson [34]: FRESCO - General purpose CRC for light and heavy ions

27. Clarke [70]: a new ‘zero-angle’ approximation for CC finite-range transfers (i.e. ‘ZA-CC’)

Note that it is sometimes difficult to put these developments in a definite chronological order.