Theory of deuteron stripping and pick-up reactions for nuclear structure studies

Deuteron stripping and pick-up experiments (d, p) and (p, d) have been used for a long time to study the structure of nuclei. Today these experiments are often carried out in inverse kinematics in state-of-the-art radioactive beams facilities around the world, extending the boundaries of our knowledge of the nuclear chart. The nuclear structure information obtained from these experiments relies entirely on transfer reaction theory. We review the theory of (d, p) and (p, d) reactions starting from early formulations and ending with the most recent developments. In particular, we describe the recent progress made in the understanding of the three-body dynamics associated with the deuteron breakup degrees of freedom, including effects of nonlocality, and discuss the role of many-body degrees of freedom within the three-body context. We also review advances in structure model calculations of one-nucleon overlap functions — an important structure input to (d, p) and (p, d) reaction calculations. We emphasize the physics missing in widely-used standard approaches available to experimentalists and review ideas and efforts aimed at including this physics, formulating the crucial tasks for further development of deuteron stripping and pick-up reaction theory. ©2019 TheAuthors. Published by Elsevier B.V. This is an open access article under the CCBY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Removing or adding one neutron from a nucleus using a (d, p) or (p, d) transfer reaction has been a popular choice over a half-century of nuclear structure studies. It originates from a suggestion made by Butler [1] and Bhatia et al [2] that useful information concerning the structure of the target and residual nuclei could be separated from the details of the reaction mechanism. According to this idea, the shapes of one-neutron transfer angular distributions are determined by the quantum numbers of the orbits populated or vacated by the neutron while their absolute values are strongly related to the single-particle (or spectroscopic) strength of final or initial target states through occupancies of orbits with these quantum numbers. This makes (d, p) and (p, d) reactions an obvious and excellent tool for testing the shell-model picture of atomic nuclei. On the other hand, the spectroscopic strength information they provide also has valuable applications to predicting the cross sections for radiative nucleon capture reactions in stars [3] or placing restrictions on the matrix elements for double-beta decay [4]. Most recently, these reactions have become an important tool in the study of neutronand proton-rich nuclei away from β-stability using radioactive beam facilities [5]. In addition, established (d, p) cross sections are demanded for data libraries by on-going strategic research programs such as ITER, IFMIF, SPIRAL2-NFS [6].
The spectroscopic strength of nuclear states is derived from (d, p) and (p, d) reactions by comparing experimental and theoretical cross sections. Experimentalists traditionally use the Distorted-Wave-Approximation theory and easily available computer codes, such as TWOFNR [7], Ptolemy [8] and DWUCK4 or 5 [9], FRESCO [10], written a long time ago. However, transfer reactions theory enjoys constant developments, which are not always reflected in the analysis of modern experiments but could still be included by a clever use of these computational tools (especially FRESCO) by treating quantities generated using state-of-the-art theories as external input. Here we offer a comprehensive review of recent progress made with these theories, with the aim of familiarizing experimentalists with new theoretical developments and thus eventually improving the quality of spectroscopic and other nuclear structure information extracted from experiments.
Another aim of this review is to provide a guide for nuclear structure theorists working to extend the application of nuclear structure theories to nuclear reactions, especially in the context of ab-initio calculations. Ab-initio approaches to the description of many-body systems are booming in all areas of physics and chemistry due to significant improvements in computing power and huge progress in high-performance computing. In nuclear physics, these approaches are now used not only for predictions of binding energies but also for more complicated quantities, including spectroscopic factors -the main focal point of interest of nucleon transfer experiments (see [11,12], for example). One nucleon overlap functions -an essential element of (d, p) and (p, d) reaction models -have now been calculated for selected nuclei in several different ab-initio approaches [13][14][15] and the first ab-initio calculations of nucleon optical potentials have been published [16][17][18]. However, providing an input from ab-initio approaches to a transfer reaction amplitude based on an oversimplified distorted-wave approximation does not make the reaction description truly microscopic.
To date only four truly ab-initio calculations of one-nucleon transfer have been published. The first two, d(d, p)t and 3 H(p, d)d, involve the simplest possible targets. They have been calculated in a four-body multichannel microscopic cluster model [19] and in an Alt-Grassberger-Sandhas approach [20,21]. These calculations describe the experimental angular distributions well, thus validating existing nucleon-nucleon (NN) potential models. The low-energy nuclear reaction t(d, n) 4 He, relevant to fusion in Tokamaks and in primordial Big-Bang nucleosynthesis, has been calculated in the no-core shell model (NCSM) combined with the resonating group method (RGM) [22], keeping, however, only the most important configurations. A more recent and more advanced calculation of this reaction, including polarization observables, can be found in [23]. A very strong basis cut-off was required to describe the low-energy 7 Li(d, p) 8 Li reaction in the same NCSM-RGM approach [24]. In fact, these calculations employed an ansatz for the nine-body wave function which means that they can only provide the best outcome for a wave function with the structure assumed in the ansatz.
To understand the nature of the limitations of these approaches we consider the widely studied reaction 40 Ca(d, p) 41 Ca at a typical incident deuteron energy of 10 MeV. This reaction involves 42 nucleons in the compound nucleus 42 Sc at an excitation energy of 20 MeV where the following binary decay channels are open: 41 Ca * + p, 41 Sc * + n, 39 Ca * + 3 He, 38 K * +α, 36 Ar * + 6 Li. Three-body decay channels are also open, the most important of which are 40 Ca * + p + n. The total number of channels is measured in hundreds and is unmanageable. Moreover, 40 Ca(d, p) 41 Ca cross sections are determined by the 42 Sc wave function in the 40 Ca-d and 41 Ca-p channels at a (macroscopic) fragment separations equal to the distance between the detector and the target. Ab-initio and other nuclear structure theories are good in predicting the wave functions at the (microscopic) femtometre scale but require explicit inclusion of binary channels to go far beyond the nuclear interior. A possible alternative, the time-dependent Hartree-Fock approach to nuclear reactions [25], has not yet been developed to include (d, p) reactions. The general theory of transfer reactions was developed in the 1960-70s and can be found in books such as [26][27][28][29]. Theoretical modelling of these reactions, including (d, p) and (p, d), necessarily involves selecting the most relevant physics associated with the reaction mechanism. Our review summarizes developments over the last decade, much of which was aimed at providing a correct description of the n-p continuum populated in these reactions. We discuss the theory of (d, p) reactions only, as (p, d) cross sections can be obtained using time reversal ideas. We concentrate on direct reactions populating discrete bound and unbound final states lying in areas of low-level density. Readers interested in compound nucleus contributions are referred to [28] while those interested in inclusive A(d, p)X reactions can find the latest developments in [30][31][32]. We do not discuss polarization and analysing power measurements where little new has been done recently on the experimental side here. On the theoretical side, a recent paper [33], based on the adiabatic theory discussed in Section 3 predicts a new type of spin-dependent force in the deuteron channel but no numerical calculations assessing their importance have been published.
We start our review in Section 2 by recalling the origins of the widely-used Distorted-Wave Approximation, highlighting outstanding issues associated with this model. In Section 3 we discuss the contribution of deuteron breakup as treated in three-body models of (d, p) reactions. Emphasis is placed on the correct account of many-body degrees of freedom within the three-body context. Both exact and approximate methods of solving the three-body Schrödinger equation are discussed. Section 4 is devoted to the Faddeev formalism for solving three-body problem in the context of (d, p) and (p, d) reactions. The extension of three-body methods to include target excitations is reviewed in Section 5. Section 6 is devoted to progress in extending three-body models to include nonlocality of the nucleon optical potentials while Section 7 discusses how the energy-dependence of nucleon optical potential can be treated in three-body calculations. Coupled-reaction-channel calculations are discussed in Section 8 and nuclear structure models for one-nucleon removal overlap functions are reviewed in Section 9. Section 10 concentrates on (d, p) reactions populating unbound states in the final nucleus. Conclusion and outlook are given in Section 11.

Main features of the distorted-wave theory
The huge number of experimental papers published over the last 50-60 years that base their (d, p) analysis on the Distorted Wave Born Approximation (DWBA) usually have a theory section in which they quote the following expression for the DWBA amplitude T DWBA dp as the starting point: where χ d is the distorted wave describing deuteron-target elastic scattering in the entrance channel, χ p is the distorted wave describing the relative motion between the proton and the final nucleus, ϕ is the single-particle wave function of the neutron transferred from the deuteron to one of the mean-field target states, V np is the neutron-proton interaction potential in the incident deuteron, ψ d is the deuteron ground state wave function and S is the spectroscopic factor corresponding to the state of the final nucleus. The definition of the coordinates r, R d , R p and r n is shown in Fig. 1 and in the caption. In this Section we recall the origin of Eq. (1) within a complete theory of the (d, p) reaction and discuss the approximations involved. We emphasize some open problems and summarize discussion in the literature on whether the spectroscopic factor S is indeed a quantity that is determined by the (d, p) reaction.

Exact formulation and antisymmetrization
The state of the A + 2-body system in a coordinate system in which total momentum is zero is described by a wave function that depends on A + 1 internal linear-independent coordinates, for example, Jacobi coordinates ξ i = √ i/(i + 1)( ∑ i j=1 r j /i − r i+1 ), where r i is individual coordinate of nucleon i in arbitrary system. For a system of fermions this wave function must be antisymmetrized with respect to nucleon permutations. Permutations can either be performed between neutrons and protons separately or, in the case of isospin formalism that treats neutron and proton as two different projection of the same isospin- 1 that is antisymmetric in the A target nucleons but treats the two extra nucleons that form the incident deuteron as distinguishable from the target nucleons. We will use the short notation (n, p, A) to denote coordinates of Ψ (+) k d throughout the paper for the case of different neutron and proton identities where A stands for all internal coordinates of nucleus A and we will make comments about the isospin formalism where appropriate.
We first discuss the partially antisymmetrised theory before returning later to the fully antisymmetrised case. Ψ (+) k d (n, p, A) satisfies the many-body Schrödinger equation where T is the kinetic energy operator in the centre-of-mass system and E is the total energy of the (A+2)-nucleon system, and satisfies the boundary conditions corresponding to a plane wave in the incident deuteron channel and outgoing spherical waves in the elastic and all other open channels. All distorted-wave theories assume that nucleons interact via two-body potentials V ij although it has been known for long time that three-body nucleon potentials are needed to explain structure and reaction observables in few-body systems. We will consider traditional reaction theory based on two-nucleon interactions only and discuss a recent work on non-trivial effects from the three-nucleon force in Section 2.6. Differential cross section measurements use proton detectors placed at values of R p , the distance of the proton from the target, that are very large on a nuclear scale. For R p → ∞ in the direction (θ , φ), the projection of the partially antisymmetrised state Ψ (+) k d (n, p, A) onto a particular intrinsic state ψ B of nucleus B, made of nucleons in A and the neutron n belonging the incoming deuteron, has the form where φ p is the spin wave function of outgoing proton p and the integration is over a complete set of internal variables τ B of B. Standard textbooks e.g., [26], Section 1.2 [34], Ch.3, [35], Vol. I, Ch. X, Section4, [36], Ch.11, use a time-dependent wave-packet treatment to show that the differential cross section for the reaction A(d, p)B to a particular state B is related to the amplitude f dp (θ , φ) by dσ dp where v dA and v pB are the deuteron and proton velocities relative to the targets in the entrance and exit channels, respectively.
Except in the case of explicit calculations involving very few target nucleons it is not practical to obtain the amplitude f dp directly from the definition Eq. (3). The following integral representation derived, for example, in [27], Ch.2, is used instead: where µ pB is the reduce mass of the p + B system. This amplitude contains a distorted wave χ kp (R p ) obtained from a two-body Schrödinger equation with an auxiliary potential U p . The expression (5) still requires knowledge of the total many-body function, but only in a restricted region of the A + 2-nucleon configuration space. It allows reaction models to be developed based on physical ideas about the most important contributions to the reaction mechanism. For example, although Ψ (+) k d (n, p, A) contains components corresponding to many open (and closed) channels, in the direct reaction regime the deuteron elastic scattering channel component is expected to produce the largest overlap with the final state in the matrix element (5).
The scattering state needed for a fully antisymmetrised theory can be obtained from Ψ (+) k d (n, p, A) and is given by [27] Ψ (as)(+) k d where N n = (N A + 1) −1/2 and N p = (Z A + 1) −1/2 . The operator P nα exchanges the neutron in the deuteron with one of the N A target neutrons labelled α and P pβ exchanges the proton in the deuteron with one of the Z A target protons labelled β.
The factors N n and N p in Eq. (6) guarantee that a wave-packet constructed from this state will correspond to an incident flux of deuterons in the incident channel that is the same as the flux calculated using the non-antisymmetrised state when the deuteron and the target are far from each other. The neutron and proton exchange terms in Eq. (6) have very different physical effects. In the first place, in the calculation of the projection needed in Eq. (3), because the state ψ B is antisymmetric in the target neutron coordinates, all the N A neutron exchange terms give the same result as the first, partially symmetrized, term but with an additional minus sign, and the neutron exchange terms can be fully accounted for by multiplying the partially symmetrized term by N A + 1. Combining this factor N A + 1 together with the factor N n gives a net factor √ N A + 1 in the (d, p) amplitude and completely takes into account neutron antisymmetrization. The effect of proton antisymmetrisation is quite different and cannot be taken into account by incorporating a numerical factor into the partially antisymmetrised amplitude. All Z A terms coming from the P pβ proton exchanges in Eq. (6) do give the same contribution to the projection (3) but they are all different from the partially antisymmetrised term. They describe a distinct physical process in which one of the target nucleons interacts with all the remaining Z A + N A + 1 nucleons and then appears in the outgoing channel, while proton p that was in the incident deuteron ends up in nucleus B. The complete antisymmetrised transition matrix is where N A + 1 is the number of neutrons in B and T dp is the direct matrix element calculated with the partially antisymmetrised scattering wave function Ψ (+) k d (n, p, A) as in Eq. (5). An explicit expression for T ex can be found in [34], page 837, Eqs. (553) and (554). Note that in calculating the (d, p) cross section the factor N p from Eq. (6) is cancelled by the factor (Z A + 1) that arises because the outgoing proton can be any one of the (Z A + 1) protons in the system. It has become standard to neglect the proton exchange term T ex in practical calculations. For medium and high energy reactions this can be justified by the high momentum transfer involved in converting a target nucleon in a shell model orbit into an outgoing proton with the required momentum. Such a process is expected to be suppressed compared with the direct transition to final nuclear state with a large spectroscopic factor, but this argument cannot be considered a very strong one at low incident deuteron energies and/or very light targets. Proton exchange can also be enhanced for proton-rich weakly-bound nuclides such as halo nuclei at the proton drip-line.
Some light is shone on these effects by a completely different approach to antisymmetry presented in a series of papers based on a model in which the deuteron is incident on a target that is represented by a single determinant of filled single particle states [37][38][39][40]. Antisymmetrization of the nucleons in the deuteron with the target nucleons leads to a three-body model, n+p+A, but with an effective, density dependent, n−p interaction renormalized by Pauli blocking. All other antisymmetrisation effects are contained in the nucleon optical potentials that describe the nucleon-A interactions. These V np modifications were found to produce small effects on elastic scattering and stripping cross sections, mainly because absorption tends to suppress contributions from the nuclear interior where Pauli blocking occurs most strongly. Note, however, that in these calculations the crucial absorption has been inserted by simply adding phenomenological imaginary parts to the real mean fields that define the model used. It is not clear whether this procedure is consistent with many-body theory. Pauli blocking also predicts novel spin-dependent effects in the deuteron channel that have been searched for in experiments with polarized deuterons, see [41] and references therein.
Many microscopic nuclear structure and reaction calculations employ the isospin formalism where neutron and proton are two different projections of the same isospin 1 2 -particle, the nucleon. Following reasoning similar to that used to derive Eq. (7), one can show that the (d, p) amplitude in this formalism is also presented by two terms, The direct term T dp is given by Eq. (5) involving the wave function Ψ (+) k d (1, 2, . . . , A; A + 1, A + 2) while T ex is given by the same equation but with Ψ (+) k d (1, 2, . . . , A − 1, A + 2; A + 1, A). Here a semi-column separates antisymmterised nucleon groups. In both cases Eq. (5) involves wave function ψ B (1, 2, . . . , A, A + 1). No attempts to calculate the contribution from T ex in the isospin formalism are known. It is always neglected. The factor of √ A + 1 is usually combined with the overlap integral ⟨ψ B |ψ A ⟩ that appears in T dp and is discussed below. The factor √ 2 is cancelled by the Clebsch-Gordan coefficient that couples neutron and proton isospin states into the isospin-zero deuteron.
Before proceeding any further we have to note that antisymmetrization is fully taken into account in all ab-initio calculations of (d, p) reactions on very light targets [19][20][21][22]24]. The relative contribution from direct and exchange terms in these papers is not quantified. Such quantification would be of a great help for further developments of (d, p) reaction theory.

Separating out the DWBA term
Practical calculations start with a rearrangement of the exact amplitude (5) where ψ d and ψ A are the internal wave functions of the deuteron and A, respectively, and H is the A+2-nucleon Hamiltonian The modelling of (d, p) reactions starts by partitioning the many-body Hamiltonian H into intrinsic Hamiltonians H d and H A of the deuteron and the target A and introducing an auxiliary potential U d that depends only on the relative coordinate R d of d and A. For simplicity we will ignore the possibility that U d can also depend on the intrinsic spins of the deuteron and A, although in practice these issues are important. The Green operator in Eq. (9) satisfies the identity ıϵ where T d is the deuteron kinetic energy operator in the centre of mass system, where p d is the momentum operator associated with coordinate R d and µ dA is reduced mass in the d − A channel. Using this, Ψ where is the distorted wave generated by U d and defined by Here E d is the incident kinetic energy in the centre-of-mass system and φ k d = exp(ık d .R d ) is the incident deuteron plane wave. The potential V d that appears in the second term on the right of Eq. (13) is the sum of the interactions between the nucleons in the deuteron with the nucleons in the target A, while G (ϵ) d is the Green's function [27], Eq. (13) should be carefully distinguished from what is often called the Lippmann-Schwinger equation and has the form Ψ (+) where It is well known [42] that Eq. (16) may not have a unique solution because of the singular nature of the limit ϵ → 0 in the second term on the right-hand-side of Eq. (13). The unique solution of Eq. (13) for ϵ > 0 is the well-defined state given by Eq. (9).
Inserting (13) into (5) and taking limit ϵ → 0 we get a new expression for the (d, p) amplitude, Repeated substitution of the right hand side of (13) for Ψ (18) produces a series with higher and higher powers of , which is known as the Born series for the (d, p) amplitude. The question of whether or not this Born series convergences is discussed in [43,44] and in [27] and references therein. In practice, it is the first term in Eq. (18) which is normally retained and the corresponding reaction amplitude is called the Distorted Wave Born Approximation (DWBA). Note that the convergence of the Born series is irrelevant. All that matters for the validity of the Born approximation is that the second term on the right of Eq. (18) is much smaller than the first term for some choice of U d .

The DWBA transition operator
The DWBA provides an opportunity to model (d, p) reactions in terms of deuteron and proton distorted waves and the structure of the wave functions ψ A and ψ B of the target and daughter nuclei. The transition operator in the DWBA term of Eq. (18) can be rewritten as The contribution to the DWBA amplitude coming from the V np term in Eq. (19) contains nuclear structure information through the projection ⟨ψ B |ψ A ⟩ of ψ B into the target wave function ψ A . This projection depends on the degrees of freedom, r n , of a single neutron. This attractive feature is not shared with the ∆V term because, unlike V np , ∆V depends on the coordinates of the target nucleons.
It is often stated that the second term in (19) where V pA = ∑ i∈A V ip and m stands for an excited state of A and 0 corresponds to its ground state. The optical potential U p , according to Feshbach [45], is determined by the matrix element where V pB = ∑ i∈B V ip , Q B is the operator that projects the total p + B wave function onto excited states of B and is determined by the Hamiltonian H B of nucleus B, the kinetic energy operator T pB in the p-B system and the nucleus B binding energy E B . The term ⟨ψ B |∆V |ψ A ⟩ contains a contribution from ⟨ψ A |V pA |ψ A ⟩−⟨ψ B |V pB |ψ B ⟩, which is the difference between the p-A and p-B folding potentials. These potentials are expected to be similar for mediumto-heavy neighbouring nuclei so that neglecting the corresponding term is somewhat justified. However, there will be two more terms in ⟨ψ B |∆V |ψ A ⟩, one of which is related to the dynamical part of U p , determined by the propagator (e p − Q B V pB Q B ) −1 , and another one being ⟨ψ B |Q A V pA |φ A ⟩. Both of them provide a flux into the proton channel via excited core states A * of nucleus B. A quantitative investigation of this flux was reported in [46] using the 40 Ca(d, p) 41 Ca reaction as an example and assuming a double-magic structure for 40 Ca within a harmonic oscillator shell model. It was found out that in most cases the contribution of ∆V did not exceed 6%, but for some choices of the deuteron optical potential it was much higher. It was also pointed out that the crudity of the model would not allow generalization of the conclusions to other nuclei and a more realistic nuclear structure model should be used for a better idea of the contribution of core excitations in B to the ∆V term. Such a task may become feasible in the near future due to the rapid progress of modern microscopic approaches aimed at linking nuclear structure and nuclear reactions.
The above considerations were based on Hamiltonian H with bare NN interactions assuming that exact nuclear wave functions could be obtained in a many-body approach. In practice, modelling of (d, p) reactions is based on a Hamiltonian H P which is a projection of H onto a truncated space, typically the one that contains the initial d + A and/or final p + B partitions. Effective Hamiltonian should be based on effective interactions arising due to the coupling to neglected parts of configuration space [47]. Such interactions are expected to be complex and thus one can expect that the complex sum ∑ i∈A V ip (r ip ) could be approximately cancelled by a proper choice of a complex potential U p (R p ). In practice, the ∑ i∈A V ip (r ip ) in most cases is modelled by some potential V pA (r p ). This assumption allows ∆V to be calculated on a routine basis, for example, using the reaction code FRESCO. While U p is usually chosen as a proton optical potential describing elastic scattering in the final channel, the choice of V pA is uncertain since the p-A relative energy in the A + p + n system is poorly defined. It is not surprising that choosing for V pA an optical potential describing p-A scattering at an energy equal to E p makes the contribution ∆V to the DWBA cross sections very small, typically of the order or 1%-2%, in the most important angular region around the first maximum of the (d, p) angular distribution. However, the validity of such a choice is unknown.
A more advanced recent treatment of ∆V in [48], where noncentral terms in the V pA interaction were accounted for in a collective model, concluded that the contribution from the flux into the proton channel via the first excited state of A depends on the excitation energy E x of this state. According to [48] the effect of the core-excitation mechanism in enhancing the transfer cross section increases with decreasing E x . Its contribution is about 2%-4% and 11% for E x ∼ 3 MeV and 0.8 MeV, respectively. The introduction of an effective Hamiltonian and effective interactions requires the use of an effective interaction V eff np in the transition amplitude, including in the deuteron ground state wave function ψ d . Usually V np is chosen to be a real n-p potential that gives the correct deuteron binding energy. In fact, the fine details of this interaction are not important because the short-range nature of the strong interaction means that the (d, p) amplitude usually depends on the properties of V np mainly through the volume integral All deuteron models, from the simplified to the most modern realistic models, provide a D 0 value close to −126 MeV fm 3/2 [27,49]. It was suggested in [47] that the effective interaction V np should be complex, energy-dependent and density-dependent. In [47] the density-dependence of V np , was assumed to be the same as that of the Brueckner-Goldstone NN G-matrix used to produce the correct saturation properties in nuclear matter calculations [50]. It was then found to have a noticeable affect on the shape of the angular distribution of the 54 Fe(d, p) 55 Fe reaction and gave a better description of the experimental data within the DWBA [47]. The basic origin of the density dependence of the Brueckner-Goldstone NN G-matrix is the role of the Pauli principle in suppressing excitations of the Fermi sea through the NN tensor force [50]. However, it was pointed out in [39] that when the Pauli principle is properly taken into account in an antisymmetrized version of the (d, p) reaction the effective interaction that appears in the (d, p) transition matrix differs from the Brueckner-Goldstone G-matrix and detailed calculations produce a significantly different density dependence associated with Pauli blocking effects. These modified effects were studied in detail for both deuteron stripping and elastic scattering in [38][39][40] and failed to confirm the work of [47]. It confirmed that the strength of V eff np is reduced in the nuclear interior and may affect reactions with large Q -value corresponding to poor momentum mismatching where the contributions from the nuclear interior are normally significant. For high bombarding energies, the modifications of V np were found to be less important.
The calculation of the remnant term ∆V can be avoided if U p is chosen to be exactly equal to the ∑ i∈A V pi . The consequence of this choice in the case of infinitely heavy target A is that the distorted wave χ p is the solution of the p − B scattering problem but with the p − A optical potential rather than with the p − B one [34]. For finite nuclei this effect has been considered in [51] within a three-body model of (d, p) reactions. A detailed discussion of the amplitude based on the choice U p = ∑ i∈A V pi is given in Section 3.2 where we review the three-body model of (d, p) reactions within the many-body context.

Overlap functions, spectroscopic factors and asymptotic normalization coefficients and their observability
In the absence of the ∆V term information about the wave functions of A and B enters the DWBA amplitude via the quantity ⟨ψ A |ψ B ⟩ -the overlap function or overlap integral -which is a function of r n , the coordinate of the transferred neutron with respect to A. The definition of the overlap integral often includes the factor √ N A + 1 to incorporate the factor of N A + 1 that arises in the cross section from antisymmetrization, as discussed in Section 2.1: As in (3) and (5), the internal coordinates τ are the A − 1 and A Jacobi coordinates for A and A + 1 nucleons, respectively.
The equivalent definition of the overlap function in isospin formalism is where the τ now include an isospin coordinate as well as appropriate space and spin coordinates. The definition (24) also incorporates the factor √ A + 1 discussed after Eq. (8). One advantage of including the factor √ N A + 1 or √ A + 1 into the definition of I(r n ) is that in a simple model in which ψ A is a closed shell determinant of shell model orbitals and MeV on the radius r 0 of the Woods-Saxon potential. Global (left) and fitted (right) deuteron optical potentials were used in the calculations. Source: Figures are taken from [56]. ψ B is a closed shell plus one nucleon determinant, then I(r n ) is precisely the wave function of the extra nucleon and the corresponding cross section is indeed the 'single-particle cross section'.
In the second-quantization formulation both Eqs. (23) and (24) can be written as where ψ(τ n ) is a neutron destruction operator. This formula needs careful interpretation. See [52] for a definition of the overlap (25) in the second-quantization formalism that respects translation invariance. The overlap functions depend on angular momenta J A and J B and their projections M A and M B of nuclei A and B, respectively, and should carry these indices as well. These functions, I ℓ,j (r n ) and angular Y ℓm ℓ (r n ) parts: where φ σ is the spin function of the nucleon in B not belonging to A and the Clebsch-Gordan coefficients couple all the angular momenta. The norm of the radial part of the overlap function is called the spectroscopic factor, Interest in spectroscopic factors has been triggered by their shell model interpretation proposed in [53]. Spectroscopic factors can be expressed as reduced matrix elements of particle creation operators (see Eq. (25)). This suggests the interpretation of S J A J B ℓj as a measure of the occupancy of nucleon orbits ℓ, j in A and B. That is why the modelling of overlap functions in the form in which the function ϕ ℓj (r) is normalized to unity, became popular. This model leads to the representation of the (d, p) amplitude as given in by Eq. (1), that factorizes the cross section as the product of a spectroscopic factor and a ''single particle cross section''. The DWBA amplitude (1) is widely used in the analysis of experimental data for extracting the spectroscopic factor by dividing the experimental cross sections by theoretical ''single particle cross section'' calculated with S J A J B ℓj = 1. These ratios are then compared to the shell model predictions made with the help of widely available computer codes such, for example, as Oxbash [54] or NuShell [55].
The spectroscopic factor extracted from experimental cross sections in this way depends on the assumption made about the function ϕ ℓj (r). The latter is assumed to be a single-particle wave function of a nucleon moving in a mean field of the target A, usually described by a Woods-Saxon shape. The widely used phenomenological shell model, that gives rise to the occupancy-interpretation of the spectroscopic factor, does not provide any single-particle wave functions. All it involves are single-particle energies and two-body matrix elements of the NN interactions fitted to reproduce nuclear spectra [57]. Any uncertainties in modelling ϕ ℓj (r) propagate into S exp ℓj determined from Eq. (1). Such uncertainties are large. For example, using a standard radius r 0 and diffuseness a of the Woods-Saxon potential well to model the ϕ ℓj (r) results in cross sections that are about 25% higher than those obtained with potential geometries constrained by the Hartree-Fock calculations [58]. Fig. 2 demonstrates the dependence of the experimental spectroscopic factors on the choice of r 0 , taken from one of the most recent calculations in [56]. One sees that uncertainty in this choice can affect spectroscopic factors by up to 50%. To understand the large uncertainty of the spectroscopic factors due to the choice of r 0 it is useful to examine their exact definition (27). Unless B is a weakly-bound nucleus for which the corresponding overlap I ℓ,j (r) has a very long tail, the main contribution to the spectroscopic factor comes from the interior region of A. Fig. 3 displays the integrand of (27) and the integral from some R max to infinity for two cases, an s-wave halo nucleus 11 Be and the double-magic 16 O, taken from [59]. It shows that for halo nucleus the important contribution to S ℓj comes from the surface area of the nucleus supplemented by a noticeable contribution from the asymptotic classically forbidden region while for doublemagic nucleus 16 O the main contribution comes almost entirely via the internal nuclear region. For most other nuclei, the spectroscopic factor is determined mainly by internal contributions.
Some (d, p) reactions, e.g., sub-Coulomb reactions, often have enhanced sensitivity to the tail of the overlap function which has the form I ℓ,j (r) ≈ C ℓj W −η,l+1/2 (2κr)/r, (29) where W is the Whittaker function, κ = √ 2µS n /h, S n is the separation energy of the neutron from B and C ℓj is the asymptotic normalization coefficient (ANC). This leads to the factorization of the reaction amplitude in terms of squared ANCs rather than spectroscopic factors. If Eq. (28) is used to model the overlaps then C 2 ℓj = S ℓj b 2 ℓj where b ℓj is the ANC of the single-particle wave function ϕ ℓj (r). Slight changes in model parameters can strongly affect b ℓj and, consequently, the spectroscopic factor S ℓj determined from experiment, as indeed seen in Fig. 2. When the internal region does contribute to the (d, p) amplitude, S ℓj can be deduced more precisely if the ANC C ℓj is determined independently from a different reaction. A method to do this was originally proposed in [60] and used in [61] and more recently in [62][63][64] but the error bars associated with this method can be large. In the isospin formalism nuclei A and B have well-defined isospins T A and T B and their projection M T A and M T B and the overlap will also include the isospin Clebsch-Gordan coefficient (T A M T A 1 2 τ n |T B M T B ). In the radial expansion (26) this coefficient is embedded in I ℓ,j (r), S lj and C lj . Shell model definitions based on this formalism often make this factor explicit. This can be convenient if overlaps for isobar analog states are needed. For such cases it is sufficient to carry out a shell model calculation once and then multiply the result by an isospin Clebsch-Gordan coefficient, denoted as C , so that the norm of the overlap function in this definition is C 2 S ℓj rather than S ℓj . There could potentially be a confusion in interpretation of C as it could mean either the isospin Clebsch-Gordan coefficient (as in C 2 S) or an ANC (as in C 2 = Sb 2 ) so one should always check the context.
Over the last decade the observability of spectroscopic factors has been questioned in a series of publications [65][66][67]. It has been pointed out that spectroscopic factors are not invariant with respect to unitary transformations and, therefore, being ill-determined quantities, are not observables. Although non-observability follows directly from the exact (d, p) amplitude (5), which does not contain any spectroscopic factors at all, the claims of their ill-defined nature are exaggerated. Indeed, the many-body Schrödinger equation with correct boundary conditions does have a well-defined solution (up to an arbitrary common phase factor). Therefore, the overlap of two well-defined quantities is a well-defined quantity together with its norm, the spectroscopic factor, irrespective of its observability. Such a quantity can be sensitive to the choice of NN interactions in the Schrödinger equation and to the choice of the model space cutoff in the same manner that any other nuclear property does such as spectra, radii, quadrupole and magnetic moments, etc. This gives rise to apparent non-invariance of spectroscopic factors with respect to unitary transformations. The question of whether unitary transformations of the exact Hamiltonian conserve invariance of spectroscopic factors could now in principle be answered for the lightest nuclei through direct ab-initio calculations.  [79] where all details are given.

Optical potential uncertainty in the DWBA
The main ingredients of the distorted-wave theory are the optical potentials describing elastic scattering in the entrance and exit channels. A few decades ago (d, p) and/or (p, d) experiments also reported elastic (p, p) and/or (d, d) measurements that were used to determine optical potentials at the energies needed. As experimental nuclear physics research is now mainly concentrated at radioactive beams facilities elastic scattering is not always measured or measured over a small angular range. In this situation global systematics of local optical potentials are often used to analyse the transfer data, the most popular of which are shown in Table 1. The general theory of optical potentials tells us that they should be nonlocal. A review of (d, p) reaction studies with nonlocal potentials will be given in Section 6. Before that only local potentials are considered.
There exist different families of optical potentials that give the same description of elastic scattering, which means that they all generate the same phase shifts and the corresponding asymptotic parts of the distorted wave functions. However, the internal and surface parts of distorted waves do depend on the optical potential choice, which affects the calculated neutron transfer cross sections. The dependence of (d, p) and (p, d) cross sections on the choice of optical potentials is almost always investigated in every experimental paper devoted to these reactions. In many cases uncertainties of up to 30% in the angular region of the main forward peak of the (d, p) or (p, d) cross section have been reported. For nuclei beyond the β-stability line only very limited elastic scattering data are available and optical potentials are less well known. As a result uncertainties in the nuclear structure information obtained from (d, p) reactions performed in inverse kinematics with radioactive beams can be much larger (see [76]). A new way of analysing these uncertainties, based on the correlated χ 2 matrix has been proposed and applied in [77,78]. Constraining (d, p) cross sections using Bayes' theorem has been developed in [79]. As an example, Fig. 4 shows the 95% confidence interval for elastic 48 Ca(d, d) and transfer 48 Ca(d, p) cross sections obtained in Bayesian analysis assuming 10% and 5% errors in experimental data.
All global potentials shown in Table 1, except the one labelled DOM, are obtained by fitting the real and imaginary parts of the optical potentials independently. However, it follows from the causality principle stating that a particle cannot be emitted before the incident particle is absorbed, that these two parts are related [80]. As an example of the application of this principle to proton scattering by nucleus B, the optical potential U p given by Eq. (21) is split into two parts, the folding potential ⟨ψ B |V pB |ψ B ⟩ and the energy-dependent dynamical part, which we denote here as V dyn (E). Both the real part and the imaginary part, W (E), of V dyn (E) are linked by the dispersion relation [80] where P denotes the Cauchy principal value. The dispersive optical model (DOM) potential shown in Table 1 has been fitted in [73] assuming relation (30) between V dyn (E) and W (E). The first application of this potential to (d, p) reactions in [81] demonstrated its advantages. Both the spectroscopic factors and ANCs extracted from experimental data showed less dependence on the energy of the projectile, as it should be as they are intrinsic structure properties of the nuclear state populated in the reactions. When other global optical potentials were used in the reaction analysis, these quantities were found to be energy-dependent.
Optical potentials can be calculated using microscopic models. We refer the reader to the most recent review [82] of developments in this area. Here we note that the microscopic model used most frequently in applications to (d, p) and (p, d) reactions is based on the approach of Jeukenne, Lejeune and Mahaux (JLM) [83,84]. In this approach the optical potential is constructed from the energy-and density-dependent potential U NM (ρ, E) obtained in nuclear matter calculations with realistic NN interactions. Parametrization of U NM in terms of ρ and E can be found in [85,86]. Using a local density ρ(r) in U NM (ρ, E) gives an r-dependent optical potential U(r). Examples of applications of the JLM approach to (d, p) and (p, d) can be found in [58,87,88].
At low energies ab-initio calculations of nucleon optical potentials within the coupled-cluster approach [16,17] and self-consistent Green's function theory [18] methods have been reported recently. These optical potentials lack absorption leading to a considerable overestimation of elastic cross sections. In [17] the absorption was increased artificially by setting the value of iϵ in e p in Eq. (21) to a large energy. This corresponds to averaging the cross sections over the corresponding large energy range, which is difficult to reconcile with the experimental situation. Moreover, the ab-initio optical potentials from [16][17][18] are not translation-invariant [52]. A translation-invariant theory of the nucleon optical potential has been recently proposed in [89], but no numerical implementation is yet available.

The effect of a three-nucleon force on (d, p) reaction cross sections
The (d, p) theory outlined above is based on a Hamiltonian (2) with two-body NN interactions only. However, it has been known for a long time that a three-nucleon (3N) force is vital for a correct understanding of the binding energies of the lightest nuclei [90]. A 3N force was also found to be necessary for the analysis of reactions with three-and four nucleons [91]. Since the 2000's [92,93] the inclusion of 3N forces has become standard in ab-initio nuclear structure calculations of light nuclei and in calculations of the nuclear equation of state at high density relevant to the physics of neutron stars [94]. It is clearly necessary to also include 3N forces in the Hamiltonian (2) when modelling (d, p) and (p, d) reactions.
Following the reasoning of [26,27] one can show that when the A + 2-body Hamiltonian (2) contains 3N interactions, ∑ i<j<k W ijk , the DWBA amplitude T dp becomes [95] T dp = ⟨ψ B φ p χ kp (R p where U p is an auxiliary potential discussed in Section 2.3. Following the usual practice, we assume that U p is chosen to generate the distorted wave χ kp and that This assumption should work better than the standard assumption phenomenological optical potentials are fitted to experimental elastic scattering data in which 3N effects are always present. From (32) and (31) we obtain T dp = T 2N dp + T 3N dp , where There is an important difference between T 2N dp and T 3N dp . In the former case V np does not depend on the internal coordinates of A and, therefore, T 2N dp contains the nuclear structure information via the overlap function ⟨ψ B |ψ A ⟩. In contrast, the 3N force involves the internal coordinates of A so that T 3N dp contains the form factor ⟨ψ B |  one state in A only. This introduces a new theoretical uncertainty into extracting spectroscopic factors from experimental data if the contribution from the 3N force is important. The T 3N dp contribution has been estimated in [95] assuming a zero-range 3N force where the coordinate vectors r ip and r in are depicted in Fig. 1 and τ i is the Pauli matrix in isospin space associated with nucleon i. This force is similar to the simplest part of the 3N interaction derived from chiral effective field theory (EFT) in next-to-next-to-leading-order (N2LO) [96,97]. The strength I 3 of this interaction has been calculated in [95] from the low-energy constant chosen in [97] to simultaneously reproduce light nuclei (A = 3, 4, 5) binding energies, neutron-α scattering, and neutron matter properties within chiral EFT at N2LO. It was shown in [95] that if nucleus A is double-magic then using the contact interaction (36) results in a standard DWBA calculation with a modified overlap function where ρ A and ρ (n) A are the matter and neutron density of the target A, respectively, D 0 is given by Eq. (22) and ψ d (0) is the deuteron wave function at zero separation between proton and neutron, also taken for consistency from the chiral N2LO interaction. Apart from density-dependence, the function I mod (r) also depends on the thickness of the neutron skin, with different results when adding or removing nucleons to neutron-or proton-rich nuclei [95]. This dependence comes from nucleon exchange in the Hartree-Fock model of double-magic nuclei that represent all possible core excitations The DWBA cross sections calculated for the 40 Ca(d, p) 41 Ca reaction at two deuteron incident energies with and without the contribution from the 3N force are shown in Fig. 5 and indicate that the 3N effects are noticeable. It should be noted that 3N force has other components, for example, those arising from two-pion exchange. Future work will be able to assess their significance for transfer reactions. The results discussed above are based on a nuclear Hamiltonian with a bare 3N interaction. Working with truncated Hamiltonians usually results in induced effective 3N interactions which may be stronger than bare interactions. Their contribution to the (d, p) will need further clarification. It is not obvious that such interactions will have a contact form. In addition the contributions from 3N terms may have a different structure, in which case a simple renormalization of the overlap function of the type (37) will not be valid. The calculation of contributions from more realistic 3N forces will require new expressions for the partial wave decomposition of the DWBA amplitude and new computer codes will have to be written.

Introduction
The deuteron is an n-p system bound by only 2.22 MeV. There is a large probability of finding the neutron and proton in the classically-forbidden region with n-p separations outside the range of V np . Any interaction with a third body can easily disrupt them. That is why the deuteron easily breaks up in nuclear reactions and the neutron can be transferred to the target from the n-p continuum states as well as from the deuteron ground state. This physics can be taken into account by introducing a three-body model of (d, p) reactions.
A three-body model of d + A elastic scattering was first suggested by Watanabe [98] who proposed the model where U p and U n are nucleon-A complex optical potentials and T 3 is the three-body kinetic energy operator in the overall centre-of-mass system. This Hamiltonian can only describe a process in which the target state is unchanged, but it does couple bound and different continuum states of the n − p system because of tidal forces associated with U p + U n . As applied to elastic deuteron scattering the additional assumption made in [98] was that the deuteron itself was not excited (broken-up) in the reaction. This leads to a simple prediction for the deuteron optical potential, U d , the Watanabe model The model was generalized in [99] to include energy-independent nonlocal nucleon optical potentials and gave some justification for the prescription that when local equivalent nucleon potentials are used in Eqs. (38)-(39) they should be taken at half the incident deuteron kinetic energy E d .
Going beyond the no-breakup assumption of Watanabe and using the E d /2 prescription for potentials U p and U n , Johnson and Soper [100] developed an approximate adiabatic description of the coupling to breakup degrees of freedom showing that they play important role in elastic deuteron scattering. The same approach was shown to give an improved account of (d, p) reactions on a medium mass targets [101]. Johnson and Tandy [102] developed a formalism for how this approach could be systematically improved. The simplest form of this theory, known as the Adiabatic Distorted Wave Approximation (ADWA) [103], continues to be the main tool, together with the DWBA, for the analysis of (d, p) experiments [5]. The ADWA is discussed in detail in Section Section 3.4 below. Given a three-body Hamiltonian the corresponding three-body wave function can be calculated rigorously for all values of the three-body coordinates using, for example, the Faddeev equation formalism. The first work of this kind was published in 2007 [104]. This method can be applied to the Watanabe Hamiltonian, but the transfer amplitude calculated from the Faddeev wave function by using, for example, the analog of Eq. (3) is not the same as the one discussed in Sections 2.1 and 2.3. In the subsections below we discuss three-body models of (d, p) reactions from the point of view of how they can be embedded in a full many-body theory and what are the main approximations involved in reducing the theory to a practical form. We will see that the effective three-body Hamiltonian certainly cannot have the Watanabe form.

Relation between three-body and many-body models of the reaction A(d, p)B
We consider a theory in which the identity of the neutron n and proton p in the incident deuteron are treated as distinguishable from the nucleons in the target nucleus A and we will use a schematic notation through all these sections that does not show angular momentum couplings between states of B or A with those of p or d. We start with the exact expression (5) and choose U p equal to ∑ i∈A V ip . Following the reasoning of [51] and [105], the A(d, p)B transition matrix becomes This equation involves two different scattering wave functions satisfying different equations and boundary conditions: describing the many-body scattering state corresponding to a deuteron with momentum k d incident on nucleus A in its ground state. It satisfies Eq. (2), which we write here more explicitly as where H A is the intrinsic Hamiltonian of the target and v nA and v pA are the sums of two-body interactions between n and p and the nucleons in A. We define the three-body state |Ψ (+) k d (p, n)⟩ as the projection of the many body state |Ψ (+) The three-body wave function Ψ (+) k d (n, p) has outgoing waves that describe elastic deuteron scattering and elastic deuteron breakup exactly. It is determined by an effective Hamiltonian, H eff , which will be discussed in the subsection below. The of the A + 2 many-body scattering state is, where Q A = 1 − P A projects onto excited states of the target. This component contains excited target states and is determined by off-diagonal matrix elements of the operator U defined in subsections below together with the propagator e. The complete many-body scattering wave function is given by The (d, p) transition matrix element now reads The amplitude (46) is different from the usual DWBA in the sense that it includes recoil excitation and breakup (REB) kp,B (p, B) contains components in which the relative motion of A and n in the nucleus B is in its ground state or excited to one of its bound or continuum breakup states. The latter components are mixed in by the p − A interaction, i.e., by the recoil of A being transmitted to n through V nA . In a similar way to Ψ (+) the particular final state appropriate to the A(d, p)B reaction of interest, and orthogonal components in which B is in all other states. An operator that generates these two components can be defined, but we do not need the details here. We use the term ''three-body model'' to mean the approximation to the transition matrix T (1) dp obtained by retaining only the (40). In this way we obtain T 3body dp This amplitude contains a single-nucleon overlap function ⟨ψ B |ψ A ⟩ as defined in Section 2.4. This is precisely the situation that leads to a proportionality between the (d, p) cross section and the spectroscopic factor associated with this overlap. The term T (2) dp , which involves explicitly paths to the final state in which the target is first excited before the deuteron is stripped, will be discussed in Section 5. We will discuss in the next subsection how excited states of A also contribute in an important way to T (1) dp .
In the rest of this section we will address work based on the assumption (48) that neglects the REB effects. See, however, [51] where important contributions from coupling to excited continuum states of B are found for a halo nucleus 11 Be populated in the 10 Be(d, p) reaction. More detailed investigation of this effect has been carried out in [106]. Note that the amplitude (48) includes multistep population of B via exciting states (bound or unbound) because Ψ (+) k d (n, p) has non-zero projections to these states.

The effective three-body Hamiltonian
The three-body wave function Ψ (+) k d (n, p), introduced above, is determined by an effective Hamiltonian, H eff . Exact expressions for H eff are derived from Eq. (41) in [107], Appendix 1, using projection operator techniques originally described by Feshbach [45]. It is shown there that where the bra-ket notation around the operator U, which is an operator in all A + 2 coordinates of n, p and A, implies integration over the target nucleus co-ordinates to leave an operator in three-body coordinates only. The many-body operator U in Eq. (49) accounts for all the effects of the target nucleus degrees of freedom. Using for simplicity a version of the theory of Ref. [107] in which certain arbitrary auxiliary one-body nucleon potentials have been set to zero, gives an explicit expression for U [107] where N is n or p and the operator e in the energy denominator in (50) is given by in which E A is the ground state energy of the target A. The three-body energy E 3 is related to the incident centre of mass kinetic energy E d and deuteron binding energy ε d by The latter sums up all processes via excited target states and the deuteron ground and breakup states that are coupled to the target ground state by v NA and which begin and end on the target ground state. This can be seen from the multiple scattering expansion of the operator U [107]: where U (i) contains Q A /e in the ith order and The leading order contribution U (0) is associated with the operator U NA which is strongly reminiscent of the Feshbach operator for the two-body optical potential given by Eq. (21). However, U NA has a different denominator, e, which explicitly contains information about the position of and interaction with the third particle making ⟨ψ A |U NA |ψ A ⟩ an operator in the three-body coordinate space so that the N-A optical potential in a three-body system is not the same as in an isolated N-A system. In addition, U (i̸ =0) contain physical processes that correspond to excitation of the target A by n followed by de-excitation by p, that are certainly not included in the nucleon-A optical potentials. These processes have the nature of a three-body force in contrast to the two-body nucleon optical potentials which sum up processes that begin and end on the ground state and where A is excited and de-excited by the same nucleon to all orders.
In almost all practical applications the effective interaction ⟨ψ A |U|ψ A ⟩ is replaced by the sum of the n-A and p−A optical potentials, as in the Watanabe model. It is clear from the explicit expression for U given in Eq. (50) that this cannot be generally correct. Estimates of contributions to U (0) and U (1) due to three-body effects induced by target excitations are given in [107] and [108], respectively, and they will be discussed in Sections 7.2 and 7.3. These estimates were made using the adiabatic approximation. Below, we introduce this approximation and consider methods that can be used to go beyond it.

Formulation of the ADWA
The matrix element (48) requires knowledge of the three-body scattering wave function Ψ (+) k d (n, p) in a specific region of configuration space where the neutron and proton are within the range of V np . The adiabatic approximation to the calculation of this component was first introduced in [100]. The adiabatic idea was applied to simplify the treatment of the deuteron breakup spectrum by assuming that the important deuteron excitation energies involved in the reaction, where ε k is the energy of an n-p scattering state, are a small fraction of the total energy E. Arguments in support for this assumption where given in [100]. This assumption, together with the assumption of a zero-range V np force leads to the approximation where the adiabatic distorted wave χ k d (R d ) is a coherent sum of the deuteron elastic scattering and breakup components of Ψ (+) k d (n, p) with all possible energies ε k of the n-p pair. This distorted wave satisfies a one-channel two-body Schrödinger equation with the Johnson-Soper optical potential U n (R d ) + U p (R d ). Note that these arguments are not designed to provide a complete solution to the three-body problem, n + p + A.
The ADWA emphasizes the evaluation of the product V np |Ψ (+) k d (n, p)⟩, as required in Eq. (48). Formally, it can be derived from an expansion of the three-body wave function Ψ [102,109,110].
These were originally defined by Weinberg as eigenstates of the kernel Under very general conditions these states form a complete set of functions of r. They are orthogonal in the sense that where the normalization −1 for i = j is convenient for an attractive V np . The eigenvalue equation (55) features a fixed deuteron energy ε d and n-p kinetic energy operator T r . As i increases the α i increase monotonically and the state φ W i becomes increasingly oscillatory with φ W i possessing i nodes within the range of V np . At large distances all the φ W i decay exponentially like the deuteron ground state wave function. Illustrative explicit examples of the Weinberg states corresponding to a Hulthén model of V np can be found in [111]. The Weinberg basis is a natural framework for searching for solutions in a coordinate space limited by the short range of V np . The properties of the basis suggest that a deuteron-nucleus scattering state that involves only low deuteron breakup energies ε k , and hence does not oscillate rapidly within the range of V np , could be well described in terms of a small number of Weinberg states within this range. The Weinberg basis thus provides a link with the adiabatic approximation as applied in Ref. [100] to the evaluation of the component of d − A scattering state that appears in the (d, p) matrix element (48). In this basis the deuteron-nucleus scattering state Ψ while Daenick global potential is used for the deuteron distortion potentials in the DWBA case. Standard geometry was employed for the bound-state potential well and all spectroscopic factors for 56 Ni are set equal to five whole those for 48 Ca are set to one. Experimental data for 56 Ni and 48 Ca are from [112] and [113], respectively.
where the Weinberg components χ W (+) i are functions of R d . According to the orthogonality property (56) they are related From the Schrödinger equation satisfied by Ψ (+) k d a set of coupled differential equations for the χ W (+) i can be deduced in a standard fashion [102]. These differ from those frequently found in the nuclear physics literature by having constant coupling potentials that do not vanish at large distances. Methods for treating such unusual couplings are discussed in [102] and [111]. The calculations reported in [111] go beyond the ADWA and include a detailed comparison with experiment. They give an analysis of the validity of the ADWA in a particular case and appear to be the only calculations of this kind performed so far for a (d, p) reaction.
The ADWA ignores all i ̸ = 0 components in the expansion (57). As we will see below, this is a very good assumption. However, to calculate the first i = 0 Weinberg component all couplings to higher Weinberg components are also neglected in the ADWA. This gives a simple model for Ψ (+) k d (n, p), described by the expression (54) but with a finite-range V np and the adiabatic distorted wave χ ADWA(+) Here T R d is the kinetic energy operator associated with the n-p centre-of-mass coordinate R d and the adiabatic potential Here, the matrix element in the nominator implies integration over all internal degrees of freedom of A and over n and p spin coordinates together with the relative n-p coordinate r. In the zero-range limit of V np the Johnson-Tandy potential To demonstrate the role of deuteron breakup, Fig. 6 compares our DWBA and ADWA calculations for two cases: removal of a strongly bound neutron (S n = 16.643 MeV) from the double-magic radioactive nucleus 56 Ni in a (p, d) reaction and the addition of a neutron bound by S n = 5.146 MeV to a stable neutron-rich isotope 48 Ca in the (d, p) reaction. The projectile energies are 37 MeV and 10 MeV, respectively. In the first case, momentum mismatch enhances the contribution from the nuclear interior and makes the difference between the DWBA and ADWA calculations more significant. In addition, employing the Johnson-Tandy potential instead of the Johnson-Soper one also makes a difference. In the second case, where the Q-value is small, the difference between DWBA and ADWA is smaller, although it grows with the scattering angle. This difference would translate into a systematic error in spectroscopic factors and ANCs extracted using one or another theory. It was noted in [101] that contributions from low angular momentum partial waves are suppressed in the ADWA. This often leads to the dominance of contributions from the nuclear periphery thus emphasizing the role of ANCs in the (d, p) reaction. For further discussion see [114].

Continuum-discretized coupled channel approach
The continuum-discretized coupled channel (CDCC) approach was first introduced and applied to (d, p) reaction calculations in [115,116], although the earlier calculations of [100] can be regarded as a simplified CDCC calculation in which the deuteron breakup continuum is represented by a single effective channel. A comprehensive review of early research in this area has been published in [117]. A review of more recent applications of the CDCC to various nuclear reactions is available in [118]. The most recent CDCC calculations of (d, p) reactions and comparisons with DWBA and/or ADWA calculations can be found in [56,63,119,120].
In the CDCC, the three-body wave function Ψ (+) k d (n, p) is found by solving the equation where the orthonormal basis of n-p continuum bin functions φ bin where φ k (r) is the n-p scattering wave function with the relative momentum k corresponding to the energy ε k , w(k) is the weight function and the normalization factor N is defined as In a more complete specification each bin i is specified by a set of orbital, spin and total angular momentum quantum with coupling potentials U ij = ⟨φ i |U A (n, p)|φ j ⟩.  small incident deuteron energies this means including closed channels i defined by E 3 − ε i < 0. The role of the closed channels has been discussed in details in a recent work [120]. It was found that at a rather high energy of E d = 40 MeV their contribution is in the range of 0%-5%, in rare cases reaching 9%. For low incident deuteron energies (∼5 MeV) the closed channels can affect the angular distributions of (d, p) reactions populating strongly-bound states, requiring renormalization of the cross sections by up to 40%, although this effect decreases with increasing electric charge of the nucleus. The population of weakly-bound state is much less affected. The closed-channels effects are also more prominent for low orbital momentum transfer.
In [120] the transition matrix element (48) was split into its elastic transfer (ET) and breakup transfer (BT) contributions as determined by the elastic and breakup components of (62). The elastic transfer was found to be dominant, with the cross sections corresponding to BT being one or two orders of magnitudes smaller. Despite this, their interference was found to be important. The CDCC results for the ET and BT contributions were compared to ADWA calculations, which were performed by setting ε i = ε 0 everywhere in the coupled equations (65). While the ADWA reproduces the ET cross sections, populating weakly-bound states, reasonably well, it often overestimates BT cross sections, and this affects the ET-BT interference. At low energies (∼5 MeV) and large separation energies the ADWA breaks down and cannot reproduce the ET cross sections as well. A large number of plots of CDCC calculations comparing the ET and BT contributions as well as contributions from closed channels and comparisons to ADWA can be found in [121]. In many cases a typical difference between the CDCC and ADWA is less than 10%-20% but there are some cases where the ADWA differs significantly from the CDCC. An example of comparison between the CDCC, ADWA and DWBA cross sections is presented in Fig. 7 for two selected reactions, taken from Ref. [63].
The dependence of CDCC (d, p) cross sections on the choice of the geometry of the overlap functions ⟨ψ B |ψ A ⟩ used in (48) has been studied in [56,63]. It was found that, as in the ADWA, at low deuteron incident energies (up to 15-20 MeV) CDCC cross sections are determined mainly by the external parts of the T 3body dp integrand, which implies that the cross sections are proportional to the square of ANCs. The situation is different at higher energies, around 60 MeV, where internal contributions are significant. This provides an opportunity to determine spectroscopic factors consistently by using information about ANCs obtained at lower energies. Examples of such an analysis are presented in [63].

Connection between the CDCC and Weinberg basis wave functions
It is known from experience with CDCC calculations that the energy range of n-p continuum states that are coupled to the incident deuteron channel is limited to tens of MeV. For these energies we expect that inside the range of V np the wave function Ψ (+) k d (n, p) will not be a strongly oscillatory function of r and only a few terms of the Weinberg expansion (57) will be needed to evaluate the (d, p) matrix element. Note that this has nothing to do with the strength of the coupling between Weinberg components in Ψ (+) k d (n, p) or how rapidly the Weinberg expansion for Ψ (+) k d (n, p) itself converges, but rather how rapid is the convergence of the sequence of contributions to the (d, p) amplitude from the different Weinberg components. These components can be constructed from the CDCC wave functions.
It was shown in [122] that the Weinberg χ W i and CDCC χ CDCC The transformation coefficients C ij are determined entirely by the bound and scattering states of V np in the energy range of the relevant bin states: They do not depend on any other details of the reaction, such as the deuteron incident energy, the transferred angular momentum, or the structure of the target nuclei. The same transformation coefficients also appear in the transformation that quantify the contribution of each Weinberg state to a particular CDCC bin state, in the presence of the weight factor V np .
The values of C ij do depend on how the CDCC bin states were constructed, through the bin sizes ∆k i , see [122] for details. In Fig. 8 they are shown for i ≤ 5 and j < 14, respectively. The point with j = 0 shows the C 10 that connects with the deuteron ground state. One can see that the bin states in the relevant energy range, 0-40 MeV, are dominated by the first Weinberg component with only small contributions from Weinberg states i = 2 − 5. This dominance is particularly marked for the low energy continuum that is the most strongly coupled to the deuteron ground state by the breakup mechanism. At the higher continuum energies the bin states are mixtures of several Weinberg states, as was expected.
Referring to Eq. (66), the dominance of the C 1j coefficients for low n-p continuum energies will make χ W 1 the dominant Weinberg distorted wave provided the contributions from continuum bins with energies greater than of order 30 MeV are not large. This is illustrated in Fig. 9 where the Weinberg distorted waves χ W 1 , χ W 2 and χ W 3 generated for the 132 Sn(d, p) 133 Sn reaction at two different deuteron incident energies are shown for the most important partial waves.
Distorted wave cross sections that retain only one Weinberg component χ W i , the DWχ i A, are plotted in the same figure where they are compared to the CDCC calculations. As was expected the DWχ 1 A and CCDC cross sections are very similar for both deuteron energies.
The dominant role of the first Weinberg component gives a precise realization of the basic idea used in [100] that only the three-body scattering state inside the range of V np is needed in (48). By focusing on the first Weinberg component of the three-body scattering wave function, this finding opens up a way of going beyond the ADWA approximation to  the (d, p) transition matrix without solving coupled equations. This could be done perturbatively following ideas of [105] where an exact expression for the effective potential, V eff , that drives χ W (+) 1 is derived for an arbitrary U A (n, p), but with a separable model for V np . This idea is reviewed in the next section.

Perturbative corrections to the ADWA distorting potential
Given that only the first Weinberg component of V np Ψ (+) k d (n, p) is really needed to calculate the (d, p) cross sections and that the Johnson-Tandy potential V ADWA frequently generates a good first approximation to χ (+) it is natural to examine a perturbation treatment of the effective distorted potential V eff which goes beyond V ADWA given in Eq. (60). It is shown in [105] that this can be achieved relatively easily for the case of a rank-1 separable n-p potential A special property of the rank-1 separable potential is that only one Weinberg state exists and it corresponds to the that gives a complete account of the product V np |Ψ (+) k d (n, p)⟩ that appears in T 3body dp in Eq. (48). In [105] an exact expression for χ (+) k d is given that involves the product of two nonlocal operators in R d space that are themselves constructed from matrix elements in r space of the three-body Green's operator G and U A (n, p). Replacing G by an average operator reduces V eff to Johnson-Soper and Johnson-Tandy results (see [105] for details).
It is shown in [105] that to first order in where the operator has a separable denominator in coordinates R d and r. The kinetic energy operator T r that appears in G ADWA can be replaced by an average valueε that is defined so that to first order in T r −ε First order corrections to this approximation vanish ifε is chosen to bē With this choice Eq. (71) for V eff reduces to The second term on the right involves the Green function G ADWA (E + 3 −ε). As a result the leading correction to the ADWA potential has a very different character depending on whether the incident three-body energy E 3 is above or below a threshold valueε. For E 3 <ε the Green function is an exponentially decaying function of R d outside the range of ⟨f np |T r ∆U A |ψ d ⟩ and V ADWA (R d ). For E 3 >ε the Green function will behave like an oscillating outgoing wave of wave number determined by E 3 −ε and V eff (R d ) will have long range components very different from V ADWA (R d ).
The thresholdε is the energy in the n − p continuum at which the corresponding n-p scattering states first show an oscillatory dependence on r within the range of V np and differ significantly from that of the deuteron ground state wave function for the same range of r. For a general V np these states would be associated with higher Weinberg states whose coupling to the lowest Weinberg state are neglected in the ADWA. How important these corrections to the ADWA are depends on the magnitude of the coupling produced by the matrix element in Eq. (75). This coupling can be calculated for any given |f np ⟩ and U A . An explicit evaluation is provided in [105] for the case of a Yamaguchi V np [123] and the choice U A = V nA (n)+V pA (p). The connection between these semi-quantitative insights into the limitations of the ADWA and the results of detailed Weinberg expansion calculations of A(d, p)B [111] and comparisons with Faddeev methods reported in [124] can also be found in [105], Sections VI, E and F. At present, no numerical assessment of the quality of the perturbative corrections to V ADWA has been reported. Such an assessment is important for future development of (d, p) reaction theory. Before proceeding further we note that for the case of the separable V np as given by Eq. (69) the threshold energyε The latter is associated with the n-p kinetic energy averaged over the volume of space where the interaction V sep np is strong. We will meet the same quantity for a more general potential V np later in Section 6.2 and Section 7.2 where it appears in the ADWA theory with nonlocal nucleon-target optical potentials.

Faddeev formalism for (d, p) reactions
The Faddeev equations were developed originally as a way of obtaining scattering state solutions of the threeparticle Schrödinger equation that did not suffer from some of the non-uniqueness problems of the Lippmann-Schwinger equation, as discussed following Eq. (16). In the Faddeev method the scattering wave function Ψ of the three-body system corresponding to an incident channel with particle 1 with energy E 1 =¯h 2 2µ 1(23) k 2 1 in the centre-of-mass system incident on a bound state of 2 and 3 with (negative) binding energy ε 23 and wave function ψ 23 , is represented as Ψ = Ψ 1 + Ψ 2 + Ψ 3 . The three components Ψ 1 , Ψ 2 , Ψ 3 are defined as the solutions of the coupled equations [125] where E = E 1 + ε 23 and Φ(R, r) = exp(ık 1 R)ψ 23 where r is the separation of particles 2 and 3 and R is the separation of their centre-of-mass from particle 1. Note that Φ satisfies where H 0 is the operator for the total kinetic energy of the system in the centre-of-mass system and V 23 is the interaction potential between 2 and 3. The operator G 0 (E + ı0) satisfies and the notation ''E + ı0'' means The remaining quantities in (76), the operators T ij (E + ı0), satisfy The uniqueness of the solution of Eqs. (76) were first proved in [125]. A distorted wave form for the Faddeev equations for the three-body scattering problem was discussed in [126,127]. A general distorted wave Faddeev method for three-body systems was introduced in [128]. A multipole expansion was used to select parts of the interactions in rearrangement channels to serve as three-body distorting potentials in the first Faddeev equation. The consequence is a greatly reduced role for the second and third Faddeev equations. Truncating the multipole expansion and discarding the second and third equations reduces the theory to the CDCC method described above.
The Faddeev method provides the amplitude for transition between binary channels, including transfer reaction amplitudes. In binary channels the wave function has the form of a normalized two-body bound state wave function times a plane wave or a spherical outgoing way when the third particle is far away. For a (d, p) reaction, the Faddeev (d, p) amplitude T Fad dp connects the bound n-p initial state, incident on a structureless target A, to the final A-n bound state.
When written in a form familiar in (d, p) reaction theory, this amplitude is given by Eq. (5) in which Ψ (+) k d is replaced by the Faddeev wave function Ψ = Ψ 1 + Ψ 2 + Ψ 3 , ψ B is replaced by a two-body bound state wave function corresponding to a real V An two-body potential and ∑ i∈A V ip is replaced by the two-body interaction V Ap that appears in the Faddeev equations. It was shown in [129] that the Faddeev amplitude can also be rewritten in the form given by Eq. (48), that does not contain the remnant term, in which once again a normalized two-body A-n wave function is used. These alternative representations of T Fad dp show that the Faddeev method lacks any nuclear structure input and provides no justification for extracting spectroscopic factors as a ratio of experimental (d, p) cross sections and cross sections calculated from Faddeev amplitudes. However, if the Fadeev method is used as an alternative to the CDCC or ADWA to provide the wave function Ψ (+) k d (n, p) in Eq. (48) then the ratio between the experimental and "Faddeev" (d, p) cross sections could be more reasonably associated with a spectroscopic factor and compared with the results obtained with the CDCC and ADWA methods.
Perhaps the biggest problem associated with Faddeev description of (d, p) reactions is the uncertainty in choosing interaction potentials. Originally, this method has been formulated for short-range pairwise potentials V ij that reproduce bound-state and scattering data in the (ij) subsystem and has been successfully used for many years in ab-initio descriptions of various three-body systems. The case of the A(d, p)B reaction is more complicated. Any attempt to reduce the A +n+p many-fermion problem to an equivalent problem involving only three-particle degrees of freedom inevitably produces a Hamiltonian with energy-dependent, nonlocal, two-and three-body forces and involves terms that cannot be simply identified with n-A and p-A optical potentials and a bare n-p interaction potential (see discussion in Section 3.3). Even if all these terms except two-body potentials are discarded the uncertainty of choosing A-n and A-p interactions remains. In the incident channel the neutron and proton in the deuteron asymptotically have energies peaked at half of the deuteron energy and should be described by complex optical potentials. However, in the exit channel the neutron is bound to A with a real potential while the outgoing proton has an asymptotic energy very different from that in the incident channel. In fact, the three-body Faddeev theory requires information about subsystem T matrices for all subsystem energies from the three-body energy to minus infinity.  [130] (the AGS equations). The AGS equations are formulated in terms of three-particle transition operators that are smoother functions than Ψ 1 , Ψ 2 and Ψ 3 . In [104] the AGS equations were first used to calculate the cross sections for the 12 C(d, p) 13 C reaction. The problem of the interaction choice as well as the importance of the internal degrees of freedom were recognized. The excitation of the 12 C(2 + ) state was taken into account using semi-separable channel interactions. In most cases the calculated (d, p) cross sections were underestimated. In a series of papers on the Faddeev-AGS description of (d, p) reactions involving halo nuclei, published in [131][132][133], three models of the n-A and p-A interactions were suggested: local energy-independent, local energydependent and nonlocal energy-independent. The latter two models will be reviewed in the sections below. Here we concentrate on local energy-independent models where a comparison with calculations performed using other three-body methods such as the CDCC and the ADWA is possible.
To describe reactions involving a one-neutron halo nucleus a real n-A potential that reproduces the bound state spectrum of nucleus (An) is needed. This is a standard choice in considering elastic scattering p + (An) reactions where the p-A potential is complex and is taken at the proton laboratory energy. The reactions d+A → p+ (An) and p+ (An) → d+A are related by time reversal provided the energy in the centre of mass system is the same. Therefore one can calculate the cross sections of the latter using the Hamiltonian with the nucleus (An) in its ground or excited state and apply detailed balance to obtain the observables for the pick-up reaction. This is equivalent to calculating the d + A → p + (An) transfer with a real n-A potential and complex p-A interaction, which is a non-standard choice. Nevertheless, it provides quite a reasonable description of the d + A elastic scattering, as demonstrated in Refs. [132,133]. A detailed comparison between these non-standard Faddeev-AGS and the corresponding CDCC calculations has been reported in [119]. One example is shown in Fig. 10 where these calculations are labelled FAGS2 and CDCC2, respectively. In general, good agreement is obtained for low incident deuteron energies but at high energies significant discrepancies appear between the CDCC and these potentials in all partial waves but one were also performed, keeping a real n-A potential in the partial wave for which the bound state exists. Once again these two methods give very similar predictions at low incident energies but disagree at higher energies. The FAGS1 and FAGS2 give similar predictions at low incident energies as do the CDCC and CDCC2 calculations. The difference between them grows with the incident energy. A similar trend has been seen in [124] where Faddeev-AGS was compared to the ADWA. Using the E d /2 and E p proton optical potential in the proton channel within the ADWA gave the same results at low incident energies but led to significant discrepancies at higher energies. Overall, the agreement between ADWA and Faddeev-AGS for reactions with beam energies around 10 MeV/u, typical for several radioactive beam facilities, was found to be better than 10% [124]. The deviation of the ADWA from the Faddeev model increases with the angular momentum of the neutron-bound state as well as with the separation energy of the neutron. The dependence on the orbital momentum l of the transferred neutron was found to be stronger than that on the separation energy. More comparison between ADWA and Faddeev methods can be found in [134].

Extended three-body models. Core excitations
The lack of reference to internal degrees of freedom in the Faddeev description of (d, p) and (p, d) reactions can be remedied to some extent by including excited states of nucleus A explicitly [104,127,135]. In this approach, two-body matrices describing subsystems A + n and A + p are determined from coupled-channel equations with multichannel twobody interactions and potential parameters obtained by fitting scattering and bound-state data in all channels considered. If all excited states were taken into account the optical operator U from Eq. (52) would be real as the projector Q A would be zero. Given that the total (bound and continuum) number of excited states of any nucleus is infinite it is essential to simplify the problem by taking explicitly only the most important channels, others being either ignored or lumped together in a complex effective interaction.
The first numerical three-body calculations with one excited core state were presented in [104] where AGS matrix three-body integral equations were derived and real separable nucleon optical potentials were used. Low-lying excited states of nucleus A are usually a consequence of deformation. Faddeev-AGS calculations with deformed n-A and p-A potentials based on a local phenomenological optical potential CH89 were first published in [135]. However, using this phenomenological potential, which has been fitted to the experimental data assuming that the Feshbach projection P A contains only the ground state of nucleus A, is inconsistent with an extended model space that includes the first excited state, 1, in addition to the ground state, 0. To account for this, a method was proposed in [136] the essential idea of which is to subtract from the elastic amplitude the contributions that are explicitly generated in the scattering equations by the coupling to the core excited state. To achieve this, the diagonal optical potential U 00 in the elastic N-A channel was modified using the formula where g 0 = (E + i0 − H NA ) −1 , H NA is the extended N-A Hamiltonian that includes the excited state of A, V 01 and V 10 are off-diagonal coupling potentials and T 1 is the transition matrix for the core in its excited state but without coupling to the ground state. T 1 is obtained with an n-A potential that is not constrained by the data but is usually taken from a standard optical potential parameterization. The resulting potential U mod 00 is nonlocal. The application of this method to the 10 Be(d, p) 11 Be and 11 Be(p, d) 10 Be reactions showed that the subtraction mostly affects the population of the excited state in the (p, d) reaction. For the (d, p) reaction to the ground state of 11 Be the subtraction procedure has little influence on (d, p) cross sections for E d < 22 MeV but its effects become noticeable for higher energies. In the case of the (d, p) reaction populating the first excited state 11 Be(1/2 − ) the subtraction effect increases from 6% at very low E d to 14% at E d = 40 MeV.
In the case of the 10 Be(d, p) 11 Be reaction the Faddeev-AGS calculations with core excitation have been compared to the CDCC predictions in [137] for E d = 20 MeV and 60 MeV. The CDCC method allows different contributions to the transfer cross section to be traced. Four different paths can lead to the population of 11 Be: from the (I) d+ 10 Be(g.s.), (II) d+ 10 Be * , (III) d * + 10 Be(g.s.) and (IV) d * + 10 Be * configurations. For both deuteron energies the main contribution comes from the d+ 10 Be(g.s.) path and in both cases the path (IV) is the weakest. The relative role of the paths (II) and (III) depends on the incident energy (see Fig. 11). Comparing complete cross section calculations, the CDCC agrees with the Faddeev-AGS at low deuteron incident energy but a noticeable discrepancy appears at high energies (see Fig. 11), similar to that found in calculations that ignored explicit core excitations. The reason for this is not understood.
Including excitation of the 10 Be(2 + ) core in the 10 Be(d, p) 11 Be reaction reduces the cross sections [137][138][139]. As the consequence, the normalization factor needed to match the predictions (dσ /dΩ) 0 without excited core to those, (dσ /dΩ) CX , obtained with the core excitation included, is smaller than unity. In general, this renormalization factor is not Source: All calculations are taken from [138] and [137].
the same as the norm S of the 10 Be(0 + )+n component in the 11 Be(1/2 + ) state. Moreover, it was found that the ratio calculated using cross sections in the main forward peak, depends on the deuteron incident energy. Fig. 12 shows that it decreases from about 0.9 to 0.5. Such behaviour is obtained both in the Faddeev-AGS and the CDCC calculations as well as in the ADWA with excited core. A detailed CDCC analysis suggests that the deviation of R x from unity stems mostly from the interference of the direct (I) and multistep (II) transfer amplitudes with the deuteron remaining in the ground state. At low energies the path (II) is suppressed thus explaining why R x ∼ 1. This ratio is weakly affected by the transfer via the deuteron breakup states, and their effect seems to be rather well simulated by the adiabatic three-body wave function. At larger angles and for high incident energies, the effect of breakup channels becomes more significant.
The ratio R x also depends on the neutron separation energy in the final state and on the core excitation energy [138]. This is illustrated in Fig. 12 where the neutron separation energy in 11 Be has been increased from the experimental value of 0.5 MeV to 5 MeV and the excitation energy of the core decreased from the experimental value of 3.7 MeV to 0.5 MeV. A smaller core excitation energy would be expected to lead to an increased contribution from the path (II), populating 11 Be via the 10 Be * state, and lead to a stronger deviation of R x from unity as indeed seen in the right panel of Fig. 12. In Ref. [139], R x has been studied for a different reaction, 20 O(d, p) 21 O, where the first excited state 2 + in the core 20 O is located at 1.67 MeV. The energy-dependence of R x was found to be similar to that observed in the 11 Be case.
It is pointed out in [139] that by examining carefully a three-body model A + n + p that includes a core excitation partition A * + n + p it is possible to get an insight into separate core-excitation contributions of a two-and three-body nature. The latter arises from neutron (proton) exciting the target A which then interacts with the proton (neutron) and plays a role of an energy-dependent effective three-body force (E3BF). This force can be directly associated with the multiple scattering terms discussed in Section 3.3 but includes the contribution from one excited state only.
Two different sectors of the Hilbert space, corresponding to the ground and the first excited states of A and coupled by nucleon-core transition potentials, are represented in the AGS equations by two-body and three-body transition operators. By excluding terms in these operators that correspond to three-body paths it is possible to get an idea of the role of the E3BF in (d, p) reactions. The separation of the two-body and E3BF contributions has been studied in [139] for the case of the 20 O(d, p) 21 O reaction at two deuteron energies, 21 and 63 MeV. The contribution from the three-body force was found to be significant, as illustrated in Fig. 13, being mainly responsible for additional absorption. At a lower energy, 21 MeV, calculations without core excitation, scaled with some spectroscopic factor, and those with core excitation, including E3BF, give very similar results. At the higher energy of 63 MeV this does not happen, most likely due to an increased contribution from the path via 20 O * +d, as happens in the 10 Be(d, p) 11 Be case discussed above. Ref. [139] finds an analogy between core excitation in the (d, p) reactions and delta-isobar excitation in nucleons in the n-d scattering.
All the calculations discussed above have been done with only one excited core state included. In reality, most nuclei of interest have many core states which give significant contributions to the spectroscopic factors sum rules. In particular case of 11 Be, discussed above, the spectroscopic factors corresponding to the 10 Be * states near the one-neutron decay threshold are comparable in magnitude to the ground state spectroscopic factor, which results in strong population of these sestates in the 11 Be(p,d) 10 Be * reaction [140]. In 21 O, the particle-bound 4 + 1 core-state at 3.57 MeV has a spectroscopic factor which is 7 to 8 times larger than the 0 + ground state spectroscopic factor, depending on whether it is obtained in the shell model or a microscopic cluster model [141]. Including more excited core states either into the CDCC of Faddeev calculations is numerically challenging. However, it is important to clarify their role, especially for a region of nuclear chart that hosts deformed nuclei.

The (d, p) reaction with nonlocal optical potentials
Most (d, p) reaction studies use phenomenological local optical potentials that depend on the relative distance between projectile and target. However, the general theory of the optical potential, as given by Feshbach [45] and expressed by Eq. (21), tells us that these potentials are nonlocal. The origin of this nonlocality is the exclusion of all reaction channels except the elastic one from the model space describing the many-body wave function. During the scattering process, the system can leave the elastic channel at some point r to an intermediate point in the excluded model space and then come back to the elastic channel at a different spatial coordinate r ′ . A nonlocal two-body N-A potential generates a distorted wave χ NL satisfying the equation where r = r A − r N is the radius-vector between A and N, T is the corresponding kinetic energy operator and E is the energy of the N-A system in the centre of mass. The Coulomb potential V coul should be included if N is a proton. Today we see a growth of interest in nonlocal potentials, both from a phenomenological and a microscopic point of view. Until five year ago only three systematic studies of nonlocal nucleon optical potentials were known [142][143][144]. Recently, new nonlocal systematics have been developed [145][146][147][148]. Significant progress has been made in the theoretical understanding of nonlocal nucleon optical potentials from a many-body perspective. A comprehensive review of these developments can be found in [82]. Below, we review the work aimed at understanding the manifestation of nonlocality in (d, p) reactions.

A localized version of the nonlocal two-body problem and the DWBA
An approximate way of accounting for nonlocality in the DWBA was developed in [149] for the particular class of optical potentials proposed by Perey and Buck [142] and it has beed widely used in analysis of experimental data over the last 50 years. The Perey-Buck potentials have the form and are characterized by a nonlocal factor H of range β, and a complex formfactor U, usually assumed to be of Woods-Saxon shape. For this class of potentials an approximation χ 0 to the exact two-body wave function χ NL is obtained from the Schrödinger equation with a local potential U 0 loc which is the solution of the transcendental equation [142] U 0 where µ is the reduced mass of the N +A system. For proton scattering, the centre-of-mass energy E in the r.h.s. of Eq. (91) must be reduced by the local Coulomb interaction V coul (r), which can be approximated by a constant, for example, the one given in [143]. The quality of the local approximation is demonstrated in Fig. 14 for the case of p+ 41 Ca elastic scattering at E p = 20 and 50 MeV, respectively. This figure shows the cross section for the nonlocal model (84) with the Perey-Buck potential [142] calculated by us using the method of [150]. The local calculations with U 0 loc provide a good agreement with the exact solution for both energies for a wide range of angles.
Corrections to the local-equivalent model (87)- (88) have been derived in [151] where χ NL is approximated by the solution χ 1 of the equation which includes the Perey factor ) , (90) and the correction to the lowest-order local potential U 0 loc is The correction terms involve the first and second derivatives (U 0 loc ) ′ and (U 0 loc ) ′′ . In this case χ 1 (r) = f (r)ϕ(r), where ϕ is found from the local Schrödinger equation If the Schrödinger equation (84) Fig. 14 shows that the improvements obtained in elastic scattering predictions by including ∆U are more important for higher energy. The link between χ 1 and ϕ gives rise to a simple method to account for exit channel nonlocality in the DWBA description of (d, p) reactions by multiplication of the proton distorted wave obtained from a local optical model by the Perey factor f (r) in the DWBA amplitude [149]. The option to include this factor is available in reaction codes DWUCK4 (5) and TWOFNR and has been used routinely. Both χ 1 and ϕ are identical in the asymptotic region giving the same elastic cross sections obtained with χ 1 (r) are reduced in the main peak with respect to those calculated with ϕ(r) by 9% and 3%, respectively. The difference in the maximum between the local calculations with χ 1 and exact nonlocal calculations is 3% and 8%. A comparison between these calculations for the whole angular range is shown in Fig. 15. In these calculations the deuteron distorted waves were generated by the global optical potential of [74], the neutron potential well was a standard choice and spectroscopic factors were set equal to one.
Recently, a different formulation of a local-equivalent potential for Perey-Buck potentials was proposed in [152]. There, the nonlocal wave function is also approximated by a product of a Perey factor F (r) and a solution of the Schrödinger equation with the local potential U LE (r). The F (r) and U LE (r) are related by equations This system of coupled transcendental equations was not solved in [152]. Instead U LE (r) was chosen to describe the cross sections from the exact nonlocal model and then the F (r) was deduced. It is not clear from the comparison of expressions (93)-(94) with (88), (89), (90) and (91) if the methods of [151] and [152] should give the same results. No comparisons between these two methods have yet been published. Comparison between the (d, p) calculations, carried out with a U LE that fits the Perey-Buck nonelastic proton scattering and uncorrected by the Perey factor, and those using exact nonlocal distorted wave shows that the difference between them, depending on the choice of reaction system and energy of the projectile, can reach 17% (we estimated this number from Figs. 5, 6, 7 and 8 of [152]). Ref. [152] also studies the influence of nonlocality in the two-body Schrödinger equation that generates the neutron bound state wave function ϕ lj (r) for Eq. (28), needed for transfer calculations. Using the Perey factor F (r) suppresses this wave function in the nuclear interior so that the nonlocal wave function has a larger root-mean-square radius. Renormalization of the bound state corrected by F (r) increases its magnitude in the asymptotic region relevant to transfer reaction and increases the (d, p) cross section. The combined effect of including nonlocality in both the proton channel and in the neutron bound state can affect the (d, p) cross sections by up to 20% [152]. Nucleon optical potentials generated in ab-initio approaches are nonlocal. An analysis of their form in [16] lead to the conclusion that they cannot be described by the gaussian Perey-Buck form. However, in general any potential can be approximated by a sum of gaussians with different ranges. Indeed, the nonlocal DOM potential that fits neutron and proton elastic scattering on 40 Ca [145] is represented as a sum of seven terms of the Perey-Buck type with nonlocality ranges of from 0.5 to 2.0 fm. Constructing a local-equivalent potential in this case is discussed in [153]. The combined effect of nonlocality in the proton channel and the neutron bound state for this potential has also been studied in [154] for the 40 Ca(p, d) 39 Ca reaction.

Nonlocal potentials in the ADWA
The treatment of nonlocality in the proton channel discussed above can be extended to the deuteron channel in the DWBA formalism of (d, p) reactions. However, the importance of deuteron breakup means that three-body models should be used to describe the wave function Ψ (+) k d in the deuteron channel. For nonlocal nucleon optical potentials the three-body Schrödinger equation reads [155]: where J = 8 ( A+1 A+2 ) 3 . Note that [155] uses a different definition of the direction of R d . Eq. (95) uses the definition of coordinates in Fig. 1. In this subsection, we discuss the solution of this equation when the ADWA is appropriate.
In the ADWA the state Ψ , r) can be expressed in terms of the distorted wave χ (+) k d which satisfies the nonlocal two-body equation where U C is the d-A Coulomb potential and U dA is given by As in the two-body N-A case, it is possible to construct an approximate local-equivalent model of Eq. (96). According to [155,156] for nonlocal nucleon optical potentials of the Perey-Buck type, the leading-order local potential U 0 loc equivalent to U dA is the solution of the transcendental equation where µ dA is the deuteron-target reduced mass and α = (A + 2)/(A + 1), β d = √ M 2 /(2M 0 ). The quantities M 2n , n = 0, 2, that appear in Eq. (99) and define β d are given by Eq. (99) is very similar to (88). It also has the same structure as the local-equivalent d-A potential obtained in the Watanabe model with nonlocal optical potentials [99]. The new feature that appears in (99) is the numerical value of the coefficient M 0 . It has been shown in [99] that with the value M 0 ≈ 1 associated with the Watanabe model, the solution of (99) is the sum of the proton and neutron optical potentials taken at half the deuteron incident energy. The smaller value of M 0 ≈ 0.5-0.8 in the ADWA case has a consequence that ln M 0 ̸ = 0. A non-zero value of the logarithm of M 0 is associated with a shift in the energy at which the local nucleon optical potentials must be evaluated if all the energy dependence of the optical potentials was due to nonlocality [155,156]. In leading order, this shift is related to half ⟨T np ⟩ V where i.e., the kinetic energy of n-p relative motion in the deuteron ground state averaged over a spatial volume with size determined by the range of V np . The same quantity appears in the expression for the effective deuteron distorting potential in uniform nuclear matter (see discussion in Section A3 of [155]) and in a perturbative correction to the ADWA potential in a model with a separable V np [105] discussed in Section 3.7. Depending on the n-p interaction model, ⟨T np ⟩ V ranges from 88 to 247 MeV and gives rise to an energy shift in the range 31-74 MeV [49]. This means that the local-equivalent d-A potential U 0 loc should be constructed from local nucleon optical potentials evaluated at a much higher energy than the E d /2 value traditionally used in the ADWA. The optical potential depths decrease with increasing energy so that U 0 loc calculated with NN potentials that have larger ⟨T np ⟩ V are shallower. This is demonstrated in Fig. 16a for the case of the Hulthén [157] and AV18 potential [158] with ⟨T np ⟩ V = 107 and 218 MeV, respectively. The difference is mainly due to the contribution from the deuteron d-state [150]. Although this state contributes only 5%-6% to the norm of the deuteron wave function, its contribution to the deuteron potential energy is about 40% [49,150], which leads to a strong sensitivity of U 0 loc and the corresponding (d, p) cross sections to the NN model choice. Fig. 16b shows this effect for the 40 Ca(d, p) 41 Ca cross sections.
The leading-order local-equivalent d-A potential U 0 loc does not differ much from a trivially-equivalent local-equivalent , which is illustrated in Fig. 17a for the d+ 40 [150] for E d = 11.8 MeV. In (b) full nonlocal calculations are compared with the local-equivalent models in leading order (LO), using U 0 loc , and in next-to-leading order (NLO). Exact nonlocal functions are used in the proton channel. Nonlocal optical potentials are from [143] and the n-p potential is the Hulthén model. an additional local surface potential ∆U and a Perey factor. Their exact expressions can be found in [155]. Fig. 17 shows that the next-to-leading order local-equivalent 40 Ca(d, p) 41 Ca cross sections are very close to those obtained in the exact nonlocal ADWA. In the main peak they differ by less than 3%.
Comparison of (d, p) nonlocal calculations with a local ADWA that uses N-A optical potentials taken at half the deuteron incident energy can be found in [156] and [159]. Table 2 shows the percentage differences, calculated in [159], of the two cross sections in the main peak when the n-A bound and exit p-B potentials were the same in both calculations. These results are labelled ''Nonlocal Entrance rel to Local''. The combined effect of including nonlocality in all channels, also studied in [159], is shown in Table 2, labelled ''Nonlocal Relative to Local''. One sees that the changes in the cross section in the main peak can reach 43%, depending on the reaction system. Since the main peak is always used to determine spectroscopic factors and/or ANCs this table gives an idea of their uncertainty due to nonlocal effects.

Velocity-dependent potentials
Nonlocality can arise from a position-dependence of the effective mass of a nucleon moving in a nuclear medium [160]. It results in a velocity-dependent optical potential. Such potentials have been used to describe a wide variety of physics problems, including pion-nucleon, nucleon-nucleon and electron-atom scattering, and in nuclear matter calculations, as well as in describing the dynamics of electrons in semiconductors and quantum dots (see [160] for references). Recently, several papers [160][161][162][163][164] have been published where new systematics were proposed that include nucleon velocity-dependence. Table 2 Percent difference of the (d, p) cross sections in the first peak with respect to calculations with local neutron and proton potentials for a range of targets at two proton energies calculated in [159]. Comparison is done for calculations with nonlocal wave functions in both the entrance and exit channels (first and third columns) and for a nonlocal wave function in the entrance channel only (second and fourth columns). The simplest form of the Schrödinger equation describing N-A scattering and accounting for change in mass of nucleon N due to interactions with nucleons in A takes the form [160]: where V (r) is the local optical potential that can include spin-orbit interaction and ρ(r) is an isotropic function. There exists an equivalent form of this equation, obtained by dividing Eq. (102) by 1 − ρ(r): where the new potentialŨ and the velocity-dependent force F are given by equations ln(1 − ρ(r)).
Using velocity-dependent nucleon potentials to describe the deuteron channel within the A + n + p three-body model means that two representations are possible.
• Case I is based on the model of (102) with a local interaction V nA (r), V pA (r) and both first (∇ n , ∇ p ) and second (∇ 2 n , ∇ 2 p ) derivatives with respect to variables r n = R d + r/2 and r p = R d − r/2: • Case II is based on the model of (103) with local effective interactionsŨ nA (r),Ũ pA (r) and first derivatives with respect to r n and r p only: Both cases were considered within the ADWA in [165] for the reaction 40 Ca(d, p) 41 Ca at E d = 20 MeV with a ρ(r) that has the surface shape proposed in [160]. It was found that the two alternative representations gave significantly different descriptions of the (d, p) angular distributions. This difference is directly attributed to the same quantity, ⟨T np ⟩ V , that gives rise to the strong n-p high-momentum dependence in the nonlocal ADWA. The difference between cases I and II disappears when the Watanabe model is used instead of the ADWA to solve Eqs. (106) and (107) (see [165] for details). It is worth noting that treating a general nonlocal potential in the next-to-leading order using Eq. (89) is equivalent to solving a velocity-dependent problem (103) with a complex function ρ(r). It has been shown in [165] that in case II the three-body wave function Ψ , r) can be represented as , r), (108) where P N is the nucleon Perey factor and the local three-body function ϕ (+) r) satisfies a three-body Schrödinger equation with effective local potentialsŨ nA andŨ nA and an additional three-body recoil-induced potential whose contribution is small. Solving this equation in the ADWA gives for ϕ (+) k d

(R d , r) a result very similar to the standard local
Johnson-Tandy model. Corrections from nonlocality come mainly from the product P n (r n )P p (r p ) of neutron and proton Perey factors. As these factors are to be evaluated at r n ≈ r p ≈ R d the product P n (r n )P p (r p ) is very close to the Perey factor f d (R d ) for the two-body deuteron scattering problem with the Johnson-Soper potential. This justifies the application of the deuteron Perey factor in traditional adiabatic description of (d, p) reactions with local nucleon potentials such as employed, for example, in the spectroscopic factor survey in [166].

Nonlocality in the CDCC approach
Although methods for solving coupled set of differential equations for generic nonlocal potentials are available, no methods to evaluate the nonlocal matrix elements arising in the CDCC have yet been developed. However, a leading-order local-equivalent approximation has been derived in [167] for nonlocal optical potentials of the Perey-Buck type (85). In this approximation, and assuming an s-wave deuteron, the CDCC channel functions χ CDCC i (R) are found from a coupled set of differential equations with local-equivalent coupling potentials U loc ii ′ that satisfy a system of the transcendental matrix equations written for ii ′ : ii ′ X kl X li ′ + · · · = 0.
Here f The coupling potentials are given by They contain the usual bin functions φ bin i and the modified-by-nonlocality continuum bin functions Eq. (110) have been solved in [167] using Newton's method. The U loc ii ′ obtained were read into the CDCC reaction code FRESCO. Although the left-hand side of (110) contains an infinite sum of terms it was found that keeping n max = 3 was sufficient to get converged solutions U loc ii ′ with a good accuracy. The local-equivalent model (109) can be generalized to account for the deuteron d-state. This involves solving transcendental matrix equations of the type (110) but taking couplings between all angular momenta into account. Including the d-state changes the ADWA (d, p) cross sections significantly, as discussed in Section 6.2, however, it does not make much difference to the CDCC results. Because of the long range of bin functions φ i ′ andφ (n) i the coupling potentials U (n) ii ′ (R d ) are not sensitive to any short-range behaviour determined by fine details of the NN models. This is demonstrated in Fig. 18 where the CDCC calculations are compared to the ADWA for the Hulthén potential, which generates only an s-wave component in the deuteron, and the realistic Reid-soft-core (RSC) potential [168] which includes both s-and d-wave states. Both potentials give the same behaviour of the deuteron s-wave bound state and continuum bins outside the range of the NN interactions and, therefore, sensitivity to the short-range part of the deuteron wave function is significantly reduced beyond the adiabatic approximation. Fig. 18 also show that non-adiabatic effects in (d, p) cross sections calculated with nonlocal potentials are very important and larger than they are for local optical potentials. Therefore, the development of exact methods for solving nonlocal CDCC equations is very important.
An alternative approximate method to account for nonlocality in the CDCC is to start from a next-to leading order representation (89) of the nonlocal N-A scattering problem (84). Such a representation has the form of a Schrödinger equation with velocity-dependent optical potential (103) in which the force F (r) is complex. The three-body Schrödinger equation is then given by Eq. (107) and its solution is given by Eq. (108) discussed above. The CDCC expansion (62) can now be applied to ϕ (+) [169] and solved using standard techniques, for example, with the computer code FRESCO. Since for (d, p) reactions the main contribution comes from r p ≈ r n ≈ R d , the product P n (r n )P p (r p ) could be applied to the overlap function rather than to the deuteron scattering waves, which allows available reaction codes to be used without  any modifications. Fig. 19 demonstrates the role of the Perey factor for the case of 12 C(d, p) 13 The cross sections were calculated in [169] using a local equivalents of the nonlocal nucleon optical potentials from [143] and the Hulthén model for the n-p interaction. Including the Perey-effect affects the cross sections by 19% and 24% in the first maximum for populating 13 C(1/2 − ) and 13 C(1/2 + ), respectively. In the latter case, the Perey-effect changes ratios between the first and second maximum as well. The leading-order local-equivalent CDCC, defined by Eqs. (109), (110), (112) calculated with the same optical potentials predicts cross sections that are close to those obtained with traditional local CDCC with local-equivalent potentials. Therefore, it is reasonable to hope that next-to-leading order local-equivalent CDCC would be close to the standard local CDCC with the Perey factor included. The latter significantly simplifies the treatment of deuteron breakup when nonlocal nucleon optical potentials are involved.

Faddeev calculations with nonlocal nucleon optical potentials
Ref. [132] reports the first results for several A(d, p)B calculations with nonlocal n-A and p-A potentials within the Faddeev-AGS approach. In this work energy-independent nonlocal potentials of the Perey-Buck type [143] were used. This avoided the dilemma of choosing optical potentials at different energies in the entrance and exit channels. At the same time, since nonlocality induces energy-dependence in local-equivalent optical potentials, using the potential of [143] allowed the energy-dependence of optical potentials to be taken into account implicitly. The calculations with nonlocal potentials were compared to those that used local-equivalent optical potentials. However, choosing the energy at which local-equivalent potentials should be used in Faddeev equations is ambiguous. In [132] a real energy-independent potential was used for the n-A interaction to provide the correct binding of nucleus B in the exit channel and the p-A potential has been evaluated at the laboratory proton energy in the exit channel. Comparison between nonlocal and local calculations is shown in Fig. 20 for the 12 C(d, p) 13 C reaction populating the ground and first 1/2 + excited states in 13 C. Qualitatively, the difference between nonlocal and local-equivalent calculations is similar to that obtained in the CDCC calculations with and without the Perey effect discussed in the section above. It is interesting that the CDCC calculations with velocity-dependent optical potentials reproduce Fadeev-AGS predictions in the first maximum reasonably well. Fig. 20 shows that including nonlocality mostly affects large angles for the case where the final bound state 1/2 − does not have nodes and redistributes the cross section between the first and second maxima when populating the 1s 1/2 state. Similar observations can be made by examining the Faddeev-AGS results obtained in [133] for light halo nuclei 11 Be and 15 C.
Calculations with angular-momentum and parity dependent nonlocal optical potentials were introduced in [170]. They were applied to the description of the differential cross sections and deuteron analysing powers of the 16 O(d, p) 17 O reaction. It was found that angular-momentum or parity-dependence of optical potentials is not important for this particular reaction whereas nonlocality is essential for a successful description of the (d, p) differential cross section data.
Most recently, the nonlocal Faddeev-AGS approach was used to calculate the differential cross sections for the 26 Al(d, p) 27 Al reaction using a number of realistic hard-and soft-core n-p potentials [171]. The aim was to check if strong sensitivity of the 26 Al(d, p) 27 Al cross sections to the n-p interaction model, arising when nonlocal nucleon optical potentials are used in the ADWA [49], persists beyond this approximation. Only a weak dependence on the n-p force model was found in [171], typically one order of magnitude lower than in the ADWA study. This study encourages the use of the simplest models of the n-p interactions in (d, p) reaction theories that include three-body dynamics explicitly. However, the question of sensitivity of the (d, p) cross sections to the n-p model remains open since the n-A and p-A optical potentials were fixed in [171]. In reality, their dependence on the NN interactions may induce a sensitivity of (d, p) cross sections to the NN model chosen.

Explicitly energy-dependent optical potentials in three-body models of (d, p) reactions
All numerical results shown above as well as all published ADWA and CDCC calculations assume that the threebody Hamiltonian describing the deuteron channel contains energy-independent optical potentials only. These optical potentials are chosen from phenomenological optical potentials evaluated at half the deuteron incident energy. However, Feshbach theory shows that optical potentials are complicated energy-dependent operators and indeed phenomenological fits support this view. Two rather different ideas have been proposed to treat energy-dependence of optical potentials in the three-body approach. We will describe both of them.

Faddeev-Alt-Grassberger-Sandhas approach
This approach involves evaluating the two-body T -matrices T nA and T pA given by Eq. (81). Even if the N-A potential is energy-independent, the Faddeev-Alt-Grassberger-Sandhas equations require that these T -matrices have to be calculated at the two-body energies E N = E −h 2 q 2 N /2µ N , where q N is the relative momentum between nucleon N and the c.m. of the remaining nucleon-target A pair. Here, µ N is the corresponding particle-pair reduced mass, and E is the three-body energy in the over-all c.m. system. Solving the Faddeev/AGS equation requires an integration over q N . Hence, in three-body calculations the particles in all pairs scatter at two-body energies between E and −∞. If optical potentials used are chosen at a fixed energy, then they do not describe two-body scattering over the range that is relevant for the solution of the three-body Faddeev/AGS equation.
It was suggested in [131] at each value of q N the pair potential corresponding to the two-body energy E N = E − h 2 q 2 N /2µ N should be used to calculate the pair T -matrix. This idea was implemented in the Fadeev-AGS calculations in [131]. It was also pointed out in there that a similar treatment is possible within the CDCC since an alternative set of scattering equations can be derived, starting from the AGS equations, that represent the integral form of CDCC. They could be solved by discretizing the n-p continuum in the momentum space while accounting for energy-dependence of optical potentials in the same way as it could be done within the Faddeev-AGS approach. Note, however, that it is not obvious how the modified set of AGS equations that follow from the procedure suggested in [131] can be derived from a three-body Schrödinger equation that refers only to three-body degrees of freedom. It is possible that this procedure can be derived within some other framework, but the usual justification for using the AGS equations based on their ability to give an exact solution of a three-body problem does not seem to be possible for the modified AGS equations suggested in [131].
The Faddeev-AGS calculations with local energy-dependent potentials are described in [131]. For negative nucleon energies, the corresponding potentials were taken as real and energy-independent. They were tuned to support a number of bound states that correspond to the ground and excited states of the A+N nucleus. Any Pauli forbidden states were removed. Comparison have been made with a model in which only potentials in the lowest partial waves (s, p and d) were energy-dependent and tuned to reproduce the low-lying spectrum of the final nucleus. In all other partial waves the optical potentials were energy-independent and fixed at the laboratory energy of the outgoing proton. Such a treatment could be considered as energy-independent since the lowest partial waves usually do not contribute (or contribute very little) to transfer reaction cross sections.
A comparison between energy-dependent and energy-independent calculations reveals that there is not much difference between the main peaks of these cross sections. This is an encouraging outcome for spectroscopic studies since both the spectroscopic factors and ANCs are determined mainly from comparison of theoretical and experimental angular distributions in this particular angular range. However for larger angles, θ > 40 • , the Faddeev-AGS approach with energydependent potentials predicts much larger cross sections, which could be explained by the small absorption at the low N-A energies that appear in the calculations. When compared to the experimental cross sections, the calculations with energyindependent cross sections perform better. As an example, the cross sections of 13 C(p, d) 12 C and 17 O(p, d) 16 O reactions are shown in Fig. 21. The curves M1 and M2 correspond to energy-independent and energy-dependent calculations, respectively. A hybrid calculation, in which energy-dependence is included only in the s, p and d partial waves of the N-A system and chosen to provide the correct energy level spectrum, is represented by the curve denoted as M3.
The fact that a better physically motivated approach produces worse results means that something is missing in the three-body Hamiltonian. From a three-body theory point of view using energy-dependent potentials creates complications. Solutions of the Schrödinger equation for different energy eigenvalues are not orthogonal and completeness relations and unitarity could be compromised. However, we recall that describing three-body dynamics within a manybody system can be formulated as a projection of the many-body wave function onto a three-body space as described in Section 3. In the case of the two-body scattering these projections satisfy a modified completeness relation, as given for example in a microscopic many-body theory such as the self-consistent Green's function method [172]. Although analogous relations have never been formulated for projections into three-body channels, one can imagine that some modified relations should exist since many-body wave functions obtained from a many-body Schrödinger equation with energy-independent NN potentials do form an orthonormal basis.
It was already pointed out in [131] that the elimination of excitations, multiconfiguration mixing and the breakup of nucleus A from the Hilbert space leads to effective complex and energy-dependent three-body potential. This is similar to the well-known cases of scattering within three-and four-nucleon systems in which nucleons are explicitly excited to a ∆ isobar, which gives consistent effective energy-dependent NN and 3N forces [173]. However, no effective three-body force was included in [131]. Some contribution from this force is recovered in calculations with explicit core excitations [135,136,138,139] but in general it is not obvious how this can be done consistently. We will show in Section 7.3 how this can be done approximately in the ADWA.

Energy-dependent interactions in the ADWA
Treating the energy-dependence of optical potentials within the three-body problem formulated as the projection of the many-body scattering state Ψ (+) k d (n, p, A) of Eq. (41) onto a three-body A + n + p space, starts by recognizing that the two-body N-A optical potential is a matrix element of the optical operator as defined by Eq. (21) in which B is replaced by A. It is important that the optical operator is determined by the propagator which contains the nucleon energy E N and the kinetic energy associated with the relative N-A coordinate. As shown in Section 3.3, the optical operator U arising from projection of a many-body system into the A + n + p channel has a much more complicated structure than the two-body optical potential does. The simplest term U (0) in this operator includes the n-A or p-A interactions only and thus could be considered as an analog of the sum of two-body optical potentials. However, the propagator in each of these potentials contains the three-body energy E 3 , the three-body kinetic energy operator T 3 and includes the n-p interaction potential V np as well. Therefore, unlike in the two-body case, where the optical operator acts on one variable only, the two-body optical operator in the three-body system is an operator on both the neutron and proton coordinates. This observation on its own questions the application of the two-body optical potentials in the three-body model space.
When the three-body problem with a complicated three-body optical potential is solved by expanding the three-body scattering wave function Ψ (+) k d (n, p) onto a complete basis in the space of variable r, the resulting coupled equations involve an averaging of the three-body optical operator U in the matrix elements between the basis states. It was shown in [107] that choosing a Weinberg state basis and retaining only the first component of this basis (which is equivalent to the ADWA approximation) allows the averaging procedure to be done in a straightforward way that leads to a simple result.
As discussed in Section 3.4, in the ADWA the first Weinberg component of Ψ (+) k d (n, p) is found by determining the distorted wave, χ ADWA(+) k d , from the solution of the two-body Schrödinger equation (59) with the adiabatic potential V ADWA given by Eq. (60). If only the U 0 part of the optical operator is retained then V ADWA = ⟨φ 1 ψ A |U (0) |ψ d ψ A ⟩, where φ 1 is given by Eq. (98). It is shown in [107] that because of the short-range nature of φ 1 the averaging procedure results in where U opt NA is the optical model operator that describes two-body N-A scattering at an effective energy This differs from the commonly used value of half the deuteron incident energy by half of the same quantity ⟨T np ⟩ V -the n-p kinetic energy, averaged over the short range of the n-p interactions and given by (101) -that we have already met in Section 6.2 in discussing ADWA calculations with nonlocal potentials. Thus treating the energy-dependence of optical potentials requires evaluating them at a significantly higher energy than E d /2 prescription. The value ⟨T np ⟩ V depends strongly on the model chosen for the n-p interaction. The values of ⟨T np ⟩ V given in [49] give E eff values between 57 MeV and 120 MeV larger than E d /2. We have learned in Sections 6.2 and 6.4 that a strong sensitivity to the n-p model choice is the drawback of the ADWA because it emphasizes the role of the short-range part of the n-p interaction. Therefore, we must be aware of significant changes that could happen if non-adiabatic corrections are introduced within the context of a treatment of energy-dependent optical potentials. At the moment it is not clear if such corrections are possible. One possible idea, originally proposed in [105] has been discussed in Section 3.7. The Source: The calculations are from [107] and [153].
formalism proposed does show that non-adiabatic corrections depend on the same quantity ⟨T np ⟩ V , which gives the hope that including non-adiabatic corrections may reduce the n-p sensitivity of the optical potentials. Numerical ADWA calculations with explicitly energy-dependent nonlocal potentials were reported in [107] and [153] using ⟨T np ⟩ V /2 = 57 MeV, obtained from the Hulthén model. This model does not contain high n-p momenta and in comparison to other n-p models it provides a better agreement between nonlocal ADWA and the leading-order local-equivalent CDCC cross sections [167]. Nonlocal calculations for 40 Ca(d, p) 41 Ca at E d = 11.8 MeV are shown in Fig. 22 for two optical potentials, the Giannini-Ricco-Zucchiatti (GRZ) [144] and the Nonlocal Dispertive Optical Model (NLDOM) potential [145], respectively. In both cases calculations employing optical potentials evaluated at E eff gave similar differential cross sections. They were compared to Johnson-Soper calculations that used local-equivalent n-A and p-A potentials given by (88). The outcome of these two calculations is very different. The GRZ potential has an energyindependent real part and an energy-dependent imaginary one with the depth of W = 17.5(1 − exp(−0.005E N )) MeV. The (d, p) calculations employing E eff with GRZ, presented in [107] and performed for E d ∼ 11.8 MeV, require optical potentials with W ∼ 60 MeV. However, the standard Johnson-Soper E d /2 calculations involve a much smaller absorption, W ∼ 0.5 MeV, thus significantly increasing the (d, p) cross sections, as seen in Fig. 22a. However, with NLDOM, where the absorption of optical potentials is different, the situation is opposite (Fig. 22b), leading to lower Johnson-Soper cross sections. Compared to the Faddeev-AGS results, where including explicit energy-dependence leads to increased cross sections at higher angles only, the ADWA treatment of energy-dependence produces different results depending on the choice of the optical potential. For both choices a uniform increase or decrease of the (d, p) cross sections was observed. For both nonlocal optical potentials the ADWA cross sections shown in Fig. 22 overestimate the experimental cross section. We know that going beyond the ADWA could reduce these predictions. However, other terms in the optical operator U can also reduce cross sections, as indeed happens when an effective energy-dependent three-body force due to core excitation is included in the Faddeev-AGS formalism. In the next subsection we show how the contributions from all excitations of the core can be estimated in the ADWA.

An estimate of the induced three-body force contributions to (d, p) reactions in the ADWA approximation
In Eq. (52) the optical model operator U that defines the projection Ψ (+) k d (n, p) of the many-body wave function into the A + n + p channel is expressed as an infinite sum of multiple scattering terms. We review the work of [108] where it is shown how the contribution of the first two terms in U, namely U (0) + U (1) , can be estimated in the context of the ADWA.
The exact expressions for U (0) and U (1) are and where the U NA are given in Eq. (53). The term U (1) is the leading contribution to the effective energy-dependent three-body interaction generated by all the excited target states included in Q A . It was shown in [108] that because of the short-range of the state φ 1 that appears in the ADWA potential, Eq. (60), the contribution from U (0) + U (1) can be approximated by As a next step, the result given by Eq. (114) of the last subsection can be used in (119). The quantity ⟨ψ A |U opt NA (E eff )|ψ A ⟩ can be identified with the optical potential V opt NA , written as taken at effective energy E eff . In the last expression V HF NA is the Hartree-Fock term and the dynamic term V dyn NA is given by the dispersion relation (30). The result of [108] is thus that in the ADWA the effect of including the induced three-body force in the leading order is to the double the dynamical part of the optical potential if ⟨ψ A |v NA |ψ A ⟩ is identified with V HF NA , i.e., This equation provides a practical way of estimating the effect of the first term in the induced three-body force if both the separate contributions V HF NA and V dyn NA to the nucleon optical potential are known. Fig. 23 shows the effect of doubling the dynamical parts of the GRZ and NLDOM potentials on predictions of the 40 Ca(d, p) 41 Ca cross sections at E d = 11.8 MeV using the Hulthén value ⟨T np ⟩ V /2 = 57 MeV in E eff . In both cases the cross sections are reduced, by 37% and 28% for the GRZ and NLDOM cases, respectively. In this particular case the introduction of the three-body induced force brings theoretical predictions in agreement with experimental data. However, the contributions from other terms in the three-body force and from non-adiabatic effects have been neglected. The calculations of Ref. [108] show that neglecting the induced three-body force induced by target excitations in three-body calculations of (d, p) is not justified, but to account for these effects properly will necessitate going beyond the methods used here.

The Coupled Reaction Channels (CRC) method
It has been shown in Section 3 how a three-body model of the n + p + A many-nucleon system can be derived using the Feshbach projection operator technique with the operator P = |ψ A ⟩⟨ψ A | that projects onto the target ground state. However, using this model does not recover all contributions to the (d, p) reaction associated with the many-body structure of the target. Inserting the three-body channel function Ψ (46) gives only the first part T (1) dp of the transition amplitude. The second part T (2) dp , which is given by (47), involves the state Q A Ψ (+) k d (n, p, A) that describes all processes that begin with the target A in its ground state and end with the target in an excited state.
Calculating (d, p) cross sections in the Fadeev-AGS method assumes all bound states of the n+A system are normalized to unity. In most cases this contradicts the true many-body nature of the system where overlaps ⟨ψ B |ψ A ⟩ are normalized to a spectroscopic factor. Adding a single excited state of A to the model space in the Fadeev-AGS approach does not change the situation because the sum of spectroscopic factors corresponding to the overlap of a particular state of B with the ground and one excited state of A does not necessarily sum to unity, as happens in the extended Faddeev-AGS scheme.
Another way of extending the model space is to use a trial wave function of the form Ψ (+) Projection operators for rearrangement collisions that handle complications that arise because of the non-orthogonality of the two terms in Eq. (122) have been given in [174]. It was proved there that such operators exist. Explicit expressions for operators of three kinds were derived and coupled equations for channel functions were given in [174]. However, practical applications do not use the methods of [174] and instead proceed as follows. It is assumed that Ψ (+) CRC is obtained from a model Hamiltonian H eff that can be written in terms of coordinates appropriate to either the dA or the pB channels as The requirement that (E − H eff )Ψ (+) CRC = 0 within the model space leads to a pair of coupled equations for the channel functions χ dA and χ pB , where the kernels K ij are given by In this equation J ij is Jacobian for transformation from the set of coordinates in the channel i to the one in channel j.
The channel target wave function ψ i is either ψ A or ψ B depending on whether i is p or d, respectively. The integration is done over internal coordinates τ A and r in the channel d and over τ A+1 in the channel p. The kernels K ij can be split into two parts, one of which has a similar structure to the DWBA amplitude, and thus carries information about the spectroscopic strength of the state of B populated in the reaction, and another one arises due to the fact the p and d channels are non-orthogonal. Both kernels depend on the choice of the channel potentials ⟨ψ d ψ A |V d |ψ d ψ A ⟩ and ⟨φ p ψ B |V p |φ p ψ B ⟩. The non-orthogonal terms are long-ranged and very nonlocal [175]. Their contribution, although small, is not entirely negligible [176,177], being within overall uncertainties associated with the models used. However, the cumulative effect of these terms with increasing number of coupled channels can reach 50%-80% for light targets [178].
In practice Eqs. (124) are solved using phenomenologically determined potentials This requires the availability of cross sections in the (d, d), (d, p) and (p, p) channels. Without these experimental cross sections no CRC calculations can be done since they rely on the channel optical potentials that differ from traditional proton and deuteron optical potentials. In other words, these calculations do not have any predictive power. The only thing they can say about the underlying nuclear structure is whether the nuclear input is the best consistent with the wave function of a particular structure given by (122). An exception may occur in multichannel CRC calculations of lowenergy reactions on light targets, such as those performed for the 6 He(p, d) 5 He reaction in [179]. In this case a small total number of channels, manageable for CRC calculations, is involved and including all of them justifies using real channel potentials that could be derived, for example, from a folding model. Note, however, that the trial wave function used in these calculations is not fully antisymmetrized in the nucleons in the incident deuteron and the nucleons in the target. As examples, we refer to a few recent CRC calculations of (d, p) reactions published in [180][181][182]. In [181], the angular distribution of the 7 Li(d, p) 8 Li reaction, measured in inverse kinematics, has been analysed together with available data on deuteron elastic scattering while the overlap between the 8 Li and 7 Li wave functions was fixed by using spectroscopic amplitudes from the literature. A set of optical potential parameters was found that reproduced the angular distributions of both the (d, p) and the (d, d) reactions. It is not clear if an equally good fit could be obtained with any other set of overlap functions. The overlap functions were varied in the CRC calculations for (d, p) reactions on 10 Be, 12 C [180] and 11 B [182] targets populating several excited states in the final nuclei. A very good description of the angular distributions has been achieved. However, the parameters of the potential well binding the transferred nucleon where very unusual for all states except the 12 B and 13 C ground states and the first excited state of 12 B: the radius r 0 was between 1.75 fm and 1.9 fm (as compared to the standard values of 1.2-1.3 fm) while the diffuseness a could reach 0.9 fm. There is no independent confirmations from modern microscopic structure theories to such a non-standard geometry. This example shows that the results of CRC analyses should be treated with caution.

Model calculations of overlap integrals
The common feature of the DWBA, ADWA, CRC and coupled-channel approaches is the overlap function I ℓ,j (r), defined by Eqs. (23)- (26), which is a projection of the many-body wave function of B onto the ground state of A. An important feature of this function is its asymptotic behaviour, given by Eq. (29) and already discussed in Section 2.4. This asymptotic behaviour follows from the inhomogeneous equation that I(r n ) satisfies [183], where T n is the kinetic energy operator associated with variable r n and S n is the (positive) neutron separation energy. Because of the short range of the interaction V ni between the neutron n and nucleon i the right hand side of Eq. (126) goes to zero at large r n inducing exponential decrease in I ℓ,j (r n ) with a decay constant determined by S n . It is important to use the experimental value of S n to give the correct asymptotic behaviour of I ℓ,j (r n ) in calculating (d, p) reaction amplitudes because contributions from outside the nuclear interior often dominate. A popular way of modelling I ℓ,j (r) is to use Eq. (28) and to fit the depth of the Woods-Saxon potential well to reproduce the experimental neutron separation energy S n of the single-particle wave function ϕ ℓ,j (r). This procedure guarantees the correct asymptotic decrease of I ℓ,j (r). However, its shape as well as its single-particle ANC strongly depend on assumptions about the geometry of the Wood-Saxon potential well. Below, we review available calculations of overlaps based on other models.

Potential model with excited core
The wave function of the many-body nucleus B can be always represented by an expansion over a complete set of the wave functions of its many-body subsystem A: Here the notation is rather schematic since some or all states of A have nonzero spin which should be coupled with the angular momentum of I m (r n ) to get the correct spin of B. Using expansion (127) in the inhomogeneous equation (126) we obtained coupled equation for the overlaps I m (r), where S In several papers, referenced in this subsection below, the expansion (127) was truncated to include only first excited state of A apart from its ground state while matrix elements v mm ′ (r n ) were modelled by some effective n-A potentials.
It should be noted that any truncation of Eq. (127) leads to a model wave function ψ model B that is not antisymmetrized in the extra particle. Moreover, the coupled system (128) provides information only about relative normalization of the overlaps with different m and tells nothing about their absolute normalization. It is usually assumed that normalization factors of individual overlaps I m add up to unity. However, this is true only when all m are present in the expansion (127) and the factor √ B is not included in the overlap integral definition (with this factor the normalization factors add up to B). The absolute normalization can be found by comparing the transfer cross sections calculated with these overlaps to the experimental data in the same way as the spectroscopic factors are determined.
Originally, potential models with core excitations became popular in connection with nuclear structure studies of the halo nucleus 11 Be where parity inversion of the first levels 1/2 − and 1/2 + is experimentally observed. It has been found that this parity inversion can be explained by the deformation and excitation of the 10 Be core [184][185][186]. In these works, the diagonal potentials v mm were modelled as deformed Woods-Saxon potentials and the off-diagonal v mm ′ couplings were taken from a rotational model. The wave functions from [186] were used in theoretical calculations of the first transfer reaction measured in a radioactive beams facility 11 Be(p, d) 10 Be [187] and the predicted admixture of 10 Be(2 + ) in the 11 Be wave function was found to be consistent with experimental data analysed with the DWBA. The latest application of this model is to the CDCC calculations with core excitation discussed in Section 5. See also [188] for a more detailed analysis of the halo nucleus 11 Be in a coupled-channel model.
The coupling scheme (128) has been extended in [189] to include more core excited states and to account for the Pauli principle by using orthogonalizing pseudopotentials. This model, named the multichannel algebraic scattering theory, is aimed at providing both bound and continuum spectrum of B. It produces a good description of the low-energy n-A cross section. Extended coupled-channel potential models (more than one excited core state) have not yet been used in (d, p) calculations.

Ab-initio calculations of overlap integrals
Full ab-initio calculations with hard-core realistic NN interactions including a 3N force have been reported only for the lightest nuclei with A ≤ 7. One-nucleon experimental separation energies for these nuclei are very well reproduced. The overlap functions ⟨d|t⟩, ⟨d| 3 He⟩, ⟨t| 4 He⟩ and ⟨ 3 He| 4 He⟩ have been calculated in the hyperspherical harmonics (HH) approach in [190,191] and the correct asymptotic behaviour of these functions have been achieved by keeping a large number of states in the HH basis expansion. The HH results have been reproduced by the Green's Function Monte Carlo (GFMC) method in [14]. The GFMC predictions are also available for ⟨ 6 Li| 7 Li⟩, ⟨ 6 He| 7 Li⟩ and ⟨ 6 Li| 7 Be⟩. However, the calculated overlaps are accurate only up to 6 fm, after which their statistical uncertainties are large. To be used in applications involving nucleon removal reactions, these overlap integrals were approximated with a good accuracy by Eq. (28) with ϕ ℓj (r n ) fitted at 0 ≤ r n ≤ 5 fm by the solution of a two-body potential model with the separation energy set equal to that obtained in the GFMC calculation. This allows a reliable extrapolation of I ℓ,j (r n ) into the asymptotic region important for transfer reactions. The availability of both the GFMC spectroscopic factors and geometry parameters of the nucleon potential well makes these overlap functions accessible to a large number of theoretical studies within standard reaction models such as DWBA, ADWA and CDCC.
The GFMC method uses trial functions from variational Monte Carlo (VMC) calculations. Constructing the overlaps with these functions is easier. A larger variety of overlap functions obtained with the VMC method for A ≤ 12 is available at [192]. As in the case of GFMC, they all have significant statistical errors at large r n so they cannot directly be used in reaction calculations. A smoothing procedure must be applied to these functions to give a reliable approximation into the asymptotic region. An example of the application of the VMC overlaps to single-nucleon removal reactions can be found in [193]. Over the last decade, the ab-initio approach was extended to treat heavier nuclei, aiming to describe the mediummass region of the nuclear Segre chart [194][195][196]. This has become possible because of the development of the unitary-transformations-based similarity renormalization group (SRG) approach, which suppresses off-diagonal NN matrix elements leading to greatly improved convergence properties while preserving observables [197,198]. Chiral effective field interactions (χ EFT) [199] softened by the SRG have become especially popular in many-body calculations. Several applications of these interactions to the spectroscopic factors calculated within ab-initio coupled-cluster and selfconsistent Green's function (SCGF) methods have been published [11,12]. These calculations, as well as earlier SCGF calculations of [200], show that for double-magic nuclei the spectroscopic factors are reduced with respect to those given by the independent particle model that normally represents these nuclei, due to the NN correlations [11,12,200], and they are not very sensitive to the NN choice.
An application of the SCGF to the overlap functions ⟨ 13 O| 14 O⟩ and ⟨ 15 O| 16 O⟩ can be found in [15] where they are obtained as the solutions of the Dyson equation in the oscillator representation [201]. Here the self-energy Σ is the analog of the (nonlocal) optical potential at a negative neutron energy. This equation is an eigenvalue problem for S n that provides the eigenfunction I(r n ). It does not determine the norm of the overlap function (the spectroscopic factor). The latter is obtained from the first derivative of the selfenergy at E = −S n [172]. The oscillator representation has a disadvantage of insufficient convergence, within a model space accessible to ab-initio calculations, of the calculated overlap integrals in the asymptotic region important for transfer reactions. To guarantee the correct asymptotic decrease of the SCGM overlap functions the Dyson equation for the overlap function I(r), should be solved exactly. If S n is consistent with experiment then the Dyson equation together with the spectroscopic factor should provide ANCs that could be compared to experimental values obtained from peripheral transfer reactions. An important quantity related to the surface region of the overlap functions is the I(r) root-mean-square radius, which in turn is related to the radii of nuclei A and B. The χEFT interactions from [199] significantly underpredict nuclear radii leading to nuclear collapse [202]. Using the corresponding overlap in transfer calculations would significantly affects the cross sections. Therefore, all the overlap functions in [15] were rescaled in coordinate space by the same factor (i.e., by introducing only one phenomenological correction) to account for the difference between the predicted [202] and experimental r.m.s radius of 16 O. More recently, a new version of the χEFT NN interaction has been proposed in [203] where NN and 3N forces are optimized simultaneously to low-energy NN scattering data, as well as binding energies and radii of few-nucleon systems and selected isotopes of carbon and oxygen. Such a potential, referred to as NNLO sat , is appropriate for ab-initio overlap functions studies and one such example will be referred to in the subsection below. Current ab-initio model calculations based on the Dyson equation (either for the overlap or for single-particle Green's functions) neglect the centre-of-mass motion problem. In particular, it was shown in [52] that definition of the projection of a scattering state of the many-body system B onto the ground state of A, used in the SCGF method, does not satisfy translation invariance. A similar statement applies to I(r n ), 1 which is the same quantity but taken at a negative value of the last nucleon's energy.
Translation invariance of the overlaps can be easily dealt with in the NCSM [13]. This model can span a large harmonic oscillator space which, however, is still not sufficient to provide a correct asymptotic behaviour of I(r) n without a special treatment of the single-nucleon cluster channels. We will consider this treatment in the next subsection.

Explicit treatment of single-nucleon channels. Microscopic cluster model
When the asymptotic behaviour in the A + n channel is important the wave function of B = A + 1 can be represented by the ansatz where A is the antisymmetrization operator that permutes the nucleon outside A with those in A while m run over selected excited states of A. The notation in (131) is schematic as it ignores spin variables but it can be easily modified to include all the spin couplings. The important point is that it is possible to insist on correct asymptotic behaviour of the relative function g (m) AB (r n ) at large r n using either the RGM [204] or the microscopic R-matrix approach [205]. If the asymptotics in other cluster partitions of B are important then antisymmetrized products of internal wave functions and relative functions in these partitions should be added to Eq. (131). This approach, called the microscopic cluster model, has been extensively used for description of low-energy nuclear reactions of astrophysical interest [206].
With the anzatz (131) the overlap ⟨ψ A |ψ B ⟩ is given by a sum of two functions, I (m) AB dir (r n ) + I (m) AB ex (r n ), the first of which is I (m) AB (r n ) and the other one, I (m) AB ex (r n ), arises from nucleon exchange and has a short range [207]. Earlier work assumed a simple 0hω harmonic oscillator shell model for the core wave functions ψ (m) A with some effective NN interactions. Examples can be found in [208,209]. Such calculations rarely reproduce the experimental separation energy of nucleon, and therefore give the incorrect asymptotic form for I(r n ). The NN interaction is normally tuned to reproduce this energy, for example, by changing slightly its Majorana part. In general, the overlap functions are sensitive to the choice of the effective NN interaction which can affect its r.m.s. radii and ANC [207,209]. The spectroscopic factor is less sensitive to the NN interaction when the number of core-excited states in (131) is fixed. However, the number of core excitations included can also affect the spectroscopic factors.
In modern versions of the microscopic cluster model, much more complicated structure models for the core wave in the ansatz (131) are employed. At present, using the NCSM is the most advanced way to link the overlap function with the structure of the core and the latest NN and 3N forces [210]. The most recent calculations of this type are reported in [211] for the ⟨ 10 Be| 11 Be⟩ overlap, where an additional term, representing the compound 11 Be system in the NCSM model space, was added to (131). These calculations used a range of the χEFT NN+3N forces SRG-softened to accelerate the convergence of the NCSM. The calculations unavoidably involve a cutoff in the 3N regularization, which is not precisely defined and whose variation can slightly change some calculated quantities such as the neutron separation energy. In calculating the overlap ⟨ 10 Be| 11 Be⟩ that could provide the ANC directly comparable to experimental data the cutoff has been adjusted in [211] to ensure that the experimental neutron separation energy in 11 Be is reproduced, which is similar in spirit to the treatment of energy thresholds in all older microscopic cluster model calculations [206][207][208][209]. In [212] it was found that the NCSM ANC is consistent with the one obtained from the latest ADWA analysis of the 10 Be(d, p) 11 Be reaction. It should be noted that that the 11 Be spectrum is extremely sensitive to the details of the nuclear NN+3N interactions and constitutes an important benchmark for future forces. It would be interesting to know how the choice of the interaction model affects the predicted ANCs as well.

Modelling overlaps within the source term approach
Eq. (126) tells us that the overlap function I(r n ) is completely determined by the source term ⟨ψ A | ∑ i∈A V ni |ψ B ⟩. If the experimental separation energy is used in this equation then the correct asymptotic behaviour of I(r n ) is guaranteed. A useful feature of calculations based on the source term method is the possibility of using model functions, approximating ψ A and ψ B in the internal region with a good accuracy, without worrying about their behaviour outside the nucleus.
Because of the short range of the NN force any inaccurate contributions from this region will be cut off without significantly affecting the calculated overlaps. This idea works also for removing protons. In this case the Coulomb potential of the point charge should be added to both sides of Eq. (126) to cancel long-range contributions from the Coulomb part of V ij [213]. Early applications of the inhomogeneous equation idea for finding the overlap functions are reviewed in [27,214]. In these applications the potential ∑ i∈A V ij was decomposed into a mean-field part and a residual interaction (see, for example, [215]), while the mean-field wave functions were used for ψ A and ψ B assuming a valence-nucleon model space only. It was noticed that the overlap I(r n ) depends strongly on the choice of the residual interactions and the latter were suggested to be tuned to give the same norm (or the spectroscopic factor) of I(r n ) that the direct overlap (see [216] for a detailed discussion). This is a direct consequence of the NN correlations which reduce spectroscopic factors with respect to S D .
An interesting feature of the source term approach using the 0hω oscillator shell model where all nucleons are active is that once the effective interaction V eff ij is chosen it gives a very similar ratio S IE /S D for all the double-magic nuclei [216,217], suggesting a universal nature for the effective interactions for one-nucleon removal overlap functions. When applied to the open shell A < 16 light nuclei, the source-term approach (STA) predicts the asymmetry of S IE /S D with respect to looselyand strongly-bound nucleon removal [216][217][218][219]. Asymmetry in the cross section ratio σ th /σ exp , associated with reduction of the spectroscopic strength with respect to the shell model prediction is known from the eikonal model analysis of a wide range of one-nucleon knockout reactions [220]. However, analyses of (d, p) and (p, d) neutron transfer experiments do not support this finding [88,134,166]. Ab-initio calculations predict an asymmetry in the ratio of the spectroscopic factors S ab-initio to S IPM as obtained in the independent-particle model, but much weaker than that suggested by knockout reactions [11]. The asymmetry in S IE /S D obtained from STA [217,219] is stronger than the one in S ab-initio /S IPM from an-initio calculations.
It should be noted that until now the STA used only two-body effective NN interactions. Including effective 3N force may modify these results. Incorporating this force together with a further development of the STA to include open shell A > 16 nuclei could give a possible way to account for NN correlations and explore the reduction of spectroscopic strength for all nuclei accessible to shell model studies.
It should be noted that both the overlap I(r n ) and the source term ⟨ψ A | ∑ i∈A V ij |ψ B ⟩ are defined in terms of internal variables and do not depend on the position of their centres of mass. Expanding the wave functions ψ A and ψ B in the translation-invariant oscillator basis offers the easiest way of dealing with this problem when computing the matrix elements needed for the source term. The contribution from the centre-of-mass motion has been estimated in [217] by taking the A → ∞ limit. It was found there that removing the centre-of-mass motion increases both the ANCs and spectroscopic factors. The centre-of-mass contribution to the squared ANCs was found to be 20% for 16 O and 13% for 208 Pb, while spectroscopic factors were affected by the centre-of-mass motion by 12% and 3%, respectively, which is larger than the 1/A contribution naively expected. It is worth noticing that estimates from the Hartree-Fock model also predict larger than 1/A centre-of-mass effects in spectroscopic factors [221].

Pre-asymptotic abnormalities in overlap functions near drip-lines
The asymptotic behaviour given by Eq. (29) is reached at those r n where the interaction between the removed nucleon with the rest is negligible. For light nuclei this takes place around r n > 6 fm. However if the one-neutron decay threshold in nucleus B is located close to the two-neutron decay threshold then the exact asymptotic behaviour is reached only at larger r n . Physically, this situation occurs when one neutron is removed from a weakly-bound neutron-rich nucleus which has two neutrons outside the core, such as 12 Be, 15 B, 20 C. After removing one nucleon the daughter nucleus, 11 Be, 14 B, or 19 C, has a small neutron separation energy, which means that its last neutron has a large probability to be in a classically forbidden region outside its core, 10 Be, 13 B, or 18 C. However, being far away from the core the halo neutron can still be close to and interact with the first removed nucleon. This interaction can produce a long-range contribution to the source term of the inhomogeneous equation (126) leading to a slow convergence of I(r n ) to its asymptotic form. An additional contribution to the overlap integral due to the NN correlation of two weakly-bound neutrons can be obtained using the Feynman diagram approach [222].
The abnormally slow convergence to the exact asymptotic form has been confirmed by three-body calculations of the ⟨ 11 Be| 12 Be⟩ overlap in [223]. The ratio I ℓ,j (r)/W −η,l+1/2 (2κr), shown in Fig. 24a, becomes constant only at 9 fm for 11 Be(1/2 + ) as compared to 6 fm expected from a standard n-11 Be two-body potential model. For the case of the first excited state in 11 Be(1/2 − ), where the neutron separation energy is only 177 keV, and for the overlap with unbound 11 Be(5/2 + ) state this ratio still does not reach the exact asymptotic behaviour at the last point where the calculations are converged, which is r n ∼ 12 fm. The abnormally slow convergence to exact asymptotics can be modelled using a two-body potential well formed by two Woods-Saxon potentials: the first of which representing a (normal-ranged) interaction of the nucleon, removed from 12 Be, with the 10 Be core and another (long-ranged) one representing interaction of two valence nucleons when both of them are far away from 10 Be (see Fig. 24b). The long-ranged interaction is represented by a potential with a diffuseness of 1.4-1.5 fm. Including such terms leads to increased cross sections of nucleon knockout reactions from 12 Be [223] and a similar effect would be expected for transfer reactions such as (d, p) and (p, d). Other neutron-rich nuclei away from stability such as 15 B, 16,18,20 C may display a similar behaviour. Whether this phenomenon is manifest in heavier nuclei is currently unknown.

Transfer to the continuum
The final state B populated in a (d, p) reaction can be unbound with respect to one neutron emission and observed as a resonance in the n-A continuum. Resonances occur in light, medium, and heavy nuclear systems. Their energies and widths are of great interest to nuclear astrophysics since many reactions that occur in astrophysical environments proceed through resonances. Furthermore, resonances provide both challenges and stringent tests for nuclear structure models. Predicting their properties involves a proper treatment of the nuclear many-body system including the continuum. In particular, nuclei beyond the neutron and proton driplines exist only as continuum structures, either narrow or broad.
For resonance states, the overlap integral ⟨ψ A |ψ B ⟩ behaves as a n-A scattering wave function, showing oscillating behaviour beyond the nuclear interior and creating important long-ranged contributions that make evaluation of radial integrals in the A(d, p)B amplitude difficult. The tail of ⟨ψ A |ψ B ⟩ is determined by n-A phase shifts in the energy interval associated with the final state B, which are in turn related to the partial width Γ n of the n-A resonance. Thus, on the one hand, the magnitude of (d, p) cross sections populating resonance states is determined by the resonance widths. On the other hand, in some cases these widths could be determined directly from the width of peaks in excitation functions thus providing a unique opportunity to test (d, p) reaction models.
Below we review work that addresses issues related to calculation of (d, p) cross sections for unbound final states. Before that we note that many transfer reaction studies associate the ratio of the experimental (d, p) cross sections to those, calculated with the n-A single-particle potential-model scattering function, with a spectroscopic factor. However, the scattering function has a divergent norm and therefore the definition of the spectroscopic factor must be modified [224].
As in the case of a bound state B the resonance width Γ n is often represented as Γ n = SΓ s.p. n , where Γ s.p. n is the width obtained in the single-particle potential model and S is a number interpreted as a spectroscopic factor. Theoretical models using a bound-state approximation for narrow resonance states can give sensible results for S. For example, in the shell model S values are frequently related to the occupancies of shell model orbits, irrespective of whether they are bound or not. However, determining S experimentally from calculations with some Γ s.p. n will have a significant uncertainty, similar to that discussed in Section 2.4 for bound states, due to the strong dependence of Γ s.p. n on potential model assumptions. We also note that because of the absence of the Coulomb barrier neutron resonances can be narrow either because their energy is much lower than the height of the centrifugal barrier (for ℓ ̸ = 0 resonances) or because the neutron decay requires significant rearrangement of several degrees of freedom in the nucleus. The latter mechanism can also produce s-wave neutron resonances despite the absence of any centrifugal barriers.

Vincent-Fortune contour integration
Early DWBA calculations used a bound state wave function with a very small separation energy to model unbound states. However, given the importance of the contributions to the (d, p) amplitude from the asymptotic region where the bound state approximation completely fails, this approximation is difficult to justify. Vincent and Fortune proposed a different method to deal with this situation [225]. They recognize that populating unbound states is a particular case of the A + d → A + p + n reaction for which an observable is a triple cross section. The general expression for this cross section measured at a proton angle of Ω p and for an energy of E nA in the n-A relative motion is where m d and p d are the mass and the incident momentum of deuteron, E and P are the total energy and momentum of the system and i and f refer to initial and final states of the reaction. The summation is made over all final spins. The reaction amplitude T includes information about the A + n + p continuum state. In most (d, p) experiments the neutron emitted from the unbound state B * is not observed. Therefore, integrating over the energy E nA one obtains Here the integration takes place over the range of energies E nA corresponding to the actual limits accepted by detectors.
The cross section dσ ℓj /dΩ p is determined by the same expression that defines the cross section for a bound final state of B but with the bound state wave function replaced by an n-A continuum wave function in the reaction amplitude.
The continuum state wave function has the quantum numbers ℓj associated with the coordinate r n . Different values of ℓj contribute incoherently to the differential cross sections, which facilitates extraction of the contribution of a resonant state of a given ℓj from the nonresonant background of other ℓj values.
In the DWBA, the amplitude T contains an integrand which is a product of three oscillating functions, two of which χ d and χ p -are the entrance and exit channel distorted waves, respectively, and the third one is the n-A continuum wave function which at large r n behaves as sin(k nA r n + δ n )/r n . The corresponding radial integral converges very slowly.
Vincent and Fortune split this integral into two parts. The integration over the internal nuclear region is carried out using standard integration techniques, while in the region where the nuclear interaction vanishes the integration is done over a contour in the complex plane which involves asymptotically decreasing analytical continuations of the asymptotic wave functions χ p and χ d . Examples of recent application of the Vincent-Fortune method could be find in [226] where continuum states in 10 Li were populated via the 9 Li(d, p) reaction in inverse kinematics. This method was also used in [227] to study intruder states in the neutron-rich isotope 21 O which are above the neutron decay threshold. In all cases the ratio of the experimental to theoretical differential cross sections was associated with the spectroscopic factor, just by analogy with the bound final states without enquiring into the legitimacy of this procedure.

Continuum bins
An alternative description of (d, p) reactions populating continuum is representing the resonance state by the n-A continuum bins discussed in Section 3.5 for the case of n-p continuum. Since these bin states are normalizable and have properties more like a bound state, standard DWBA-like methods can be used to calculate the (d, p) amplitude. In order to account for the sharp variation of the wave function within the energy interval spanned by the resonance width the bins were normalized with a factor w(k) = exp(−δ k ) sin δ k , where δ k is the phase shift in the partial wave ℓ. For narrow resonances a justification of this procedure can be given by considering the known continuum wave function behaviour in the vicinity of the resonance energy, which allows the integration over E nA to be performed in Eq. (133). A comparison of 20 O(d, p) 21 O cross sections calculated using the negative-parity ℓ = 1 and ℓ = 3 continuum bins with those obtained by the Vincent-Fortune method resulted in a similar outcome [227].
Continuum bin methods are not suitable for populating virtual s-wave neutron resonances. In this case, the absence of the centrifugal barrier makes it impossible to construct a resonance continuum bin from the two-body potential model alone. However, populating virtual s-states in neutron-rich light nuclei is of high interest for studies of intruder state evolution beyond the edge of stability. A virtual s-state has been recently populated in the 8 He(d, p) 9 He experiment with the 8 He radioactive beam [228]. In this particular case weakly-bound states were used for the ⟨ 8 He| 9 He(1/2 + )⟩ overlap integral in the analysis of experimental data. Once again, the normalization factor determined from the ratio of experimental to theoretical cross sections was assumed to be the spectroscopic factor. It is unclear, however, how the spectroscopic factor of a virtual state is defined.
The discretized n-A continuum could be used to solve the three-body Schrödinger equation for the final state A + n + p wave function using a CDCC expansion similar to that given by Eq. (62) for the n-p continuum. The CDCC calculations of this type have been performed in [229] to investigate the d 3/2 continuum spectra of the unbound isotope 21 C populated in the 20 C(d, p) reaction. In this work, the n-20 C wave function was expanded onto Gaussian basis pseudo states. Transfer to each continuum pseudo state was calculated and a smoothing technique was applied to obtain the dependence of the (d, p) cross sections on the n-20 C energy. Despite a large range of integration in the variable r n , up to 80 fm, it was still not clear if the cross section had converged.
The above calculations were based on the post-form presentation of the (d, p) reaction amplitude used throughout all the above sections. For populating unbound states a prior-form of this amplitude, T prior dp could be advantageous [230] since it involves radial integrations that are always restricted. For the general derivation of the prior-form representation we refer to [26,27]. We note here that the presentation (134) already assumes that the many-body wave function in the final channel is factorized into a three-body wave function Ψ (−) 3B (R p , r n ) times ψ A and that all the components with excited core states A * are neglected. The difference between Ψ (−) 3B (R p , r n ) and Ψ (−) kp,B (p, B) from Section 3.2 is that the former satisfies the Schrödinger equation that contains all three potentials n-A, p-A and np while the latter omits the n-p interaction. As the result, unlike the three-body amplitude (48) that does not include any remnant term in the transition operator due to its absorption into Ψ (−) kp,B (p, B), the prior-form representation must explicitly include it. This term contains an auxiliary potential U d chosen to describe d-A elastic scattering in the entrance channel. Another potential, U pA is assumed to approximate the p − A interaction ∑ i∈A V ip . It is chosen to be the optical potential describing p-A scattering. When calculating Ψ (−) 3B (R p , r n ) by the CDCC method all potentials are evaluated at a fixed energy. The numerical application of the CDCC to populating continuum states using the prior-form method can be found in [231] where the 9 Li(d, p) 10 Li reaction is studied with several 10 Li structure models. Using this formalism it was found that experiments at two different incident energies of 9 Li radioactive beams can be described with the same 10 Li model. The formalism also helped to identify the contribution from the virtual s-wave state in 10 Li and to demonstrate the persistence of parity inversion beyond the neutron drip line in this region of the nuclear chart.

Surface formulations of the transfer amplitude
Recently, a formulation for transfer reaction amplitude has been proposed in [232,233] that covers transfers to both the bound and resonance states and is general enough to include deuteron breakup. A central point of this formalism is the recasting of the reaction amplitude, which is a volume integral, in terms of a surface integral plus (a presumed small) remnant term that contains contributions from the interior and exterior of the final nucleus. The transition amplitude T for a (d, p) reaction can be split into two parts, T = T int (0, a) + T ext (a, ∞), (135) where a is a specific value of the coordinate r n . The term T int involves integration from r n = 0 to r n = a while T ext involved the integration from r n = a to r = ∞. In the post-form DWBA, considered throughout all this review and also implied in Eq. (135), Green's theorem can be applied to convert T ext (a, ∞) into a surface integral plus a prior-form exterior matrix element T prior ext , T post ext (a, ∞) = T surf (a) + T prior ext (a, ∞).
The latter is defined by Eq. (135). The surface-integral representation (135)- (136) avoids the convergence problems that affect traditional calculations of transfers to resonances because T int (0, a) includes integration over a well-defined internal region while T prior ext (a, ∞) involves a limited range of r n beyond a due to the suppression of the integrand by fast decreasing potentials V nA or U pA − U d . The surface term T surf (a), which is evaluated at a specific radius a, can be parametrized by quantities that are familiar from traditional R-matrix approaches, namely a channel radius (here the separation radius a), logarithmic boundary terms (here logarithmic derivatives of known Hankel functions), and reduced-width amplitudes (here related to the asymptotic normalization of the overlap function) [233]. Thus, a useful link between resonance properties and transfer observables is established because the cross section obtained from the surface integral can be parametrized in terms of quantities that are familiar from R-matrix theory.
The surface-integral formulation has not yet been implemented as readily available computer codes. However, the contribution from the T int (0, a) and T prior ext (a, ∞) terms to the (d, p) DWBA cross section has been investigated for bound final states [234,235] and unbound final states [235] for various isotopes in different mass regions and for a variety of beam energies. It was found that T surf (a) dominated in the region where one expects it to be strong, near a ∼ 5-7 fm for bound states, and at slightly larger radii for resonance states. It does not completely reproduce the exact cross section for any a. Without T int (0, a) and T prior ext (a, ∞) cross section magnitudes could differ by as much as 30%-50% from those obtained in exact calculations. Increasing the surface radius a leads to improved cross sections at the expense of having a contribution T int (0, a) from the nuclear interior that has some model dependence (although less than the model dependence of traditional calculations). Alternatively, decreasing a will require the explicit inclusion of T prior ext (a, ∞) contributions in the DWBA implementation of the method. However, they could be expected to be partially accounted for by the deuteron breakup components in the CDCC extension of the surface-formalism approach [235].
Another formulation of the surface-integral method, which is based on the CDCC, has been presented in [236]. In the post-formulation of this method the amplitude T is expressed in terms of the surface integral over R p at some fixed value of R p rather than over the surface variable of r n as in the previous formulations of [232][233][234][235]. The method is expected to work well for reactions to bound final states. For deuteron stripping to resonance states the prior formulation is used in [236]. There, the amplitude T is expressed in terms of (i) the surface integral over R d , taken at some finite value R d determined by the transition operator and the deuteron bound-state wave function ψ d , and (ii) the volume integral over the variable r, limited by decreasing ψ d (r). In both the post-and prior formulations the amplitude T contains a contribution from an auxiliary matrix element, which determines the contribution from the nuclear interior and has its origin in the inconsistent treatment of V nA in the entrance and exit channels (complex optical potential describing the n-A scattering at half the deuteron incident energy and real potential supporting bound or resonance states in B). Once again, the implementation of this formulation as reaction codes is not yet available.

Conclusions and outlook
We have presented a comprehensive review of modern developments in the theory of deuteron stripping and pickup reactions as used in the analysis of experimental data. The starting point of this theory is a many-body Hamiltonian governed by the NN interaction, but the main object of interest -the cross section -is determined by the projection of an eigenfunction of this Hamiltonian corresponding to particular initial boundary conditions onto the internal wave functions of the fragments in the final reaction channel specified by the experimental observations. Because of the complexity of the direct evaluation of this projection, from the very early stages of the development of reaction theory the neutron transfer reaction amplitude has been expressed in terms of quantities that are linked to other experimental observations such as the elastic cross sections in the entrance and exit channels. The DWBA was formulated along these lines as the matrix element of a potential transition operator sandwiched between distorted waves describing elastic scattering in the initial and final-channel wave functions. It assumes that the probability of direct transitions from the initial to the final state is small enough that transitions via intermediate channels can be neglected. The latter are included implicitly via their contribution to optical potentials in the elastic channels. Because the DWBA is still used today in the analysis of stripping and pick-up data, in Section 2 we recall the origins of the DWBA and discuss the contributions to the reaction amplitude that are left out in this approximation. Section 2 is also designed to link with the many-body structure community who are interested in developing a consistent unified description of nuclear structure and nuclear reactions, especially within the ab-initio approach context. Ab-initio theories can now provide all the input (optical potentials and overlap functions) needed for DWBA and ADWA calculations. However, using this input may not provide a test of the validity of the ab-initio calculations because of contributions not included in the DWBA and ADWA such as 3N forces and proton exchange.
A distinctive feature of the A(d, p)B reaction is the contribution of neutron transfer from deuteron to the target A from A+n+p breakup components of the total wave function, which are strongly coupled to the elastic A+d channel. For almost 50 years the analysis of experimental data has included deuteron breakup within the Johnson-Soper and Johnson-Tandy adiabatic models that go well beyond the DWBA. More recently theoretical research efforts have concentrated on the exact treatment of the A + n + p problem based on a three-body description of the (d, p) reaction within either the CDCC or Faddeev formalism. They confirmed the strong coupling to the A + n + p continuum and quantified the contribution to neutron transfer from these channels. Both the CDCC and the Faddeev method have been extended to take into account the explicit coupling to selected excited states of the target A, usually the lowest state.
Over the last decade significant efforts have been devoted to incorporating the nonlocality of nucleon optical potentials into three-body models of the (d, p) reaction, using the Faddeev and ADWA methods. The approximate treatment of nonlocality has also been developed within the CDCC method. These latest developments are mainly concentrated on methods for solving the equations of three-body formulations and their benchmarking. However, an important question remains concerning the input to the three-body models. Frequently three-body models use a three-body Hamiltonian with nucleon optical potentials fixed at a particular energy, usually half the incident deuteron energy, but this does not represent the physics of the reaction correctly although reasonable (d, p) cross sections angular distributions are often obtained.
Unifying nuclear structure theories with reaction models that include deuteron breakup channels is impossible without a consistent mapping from the physical many-body problem onto a A + n + p three-body model. The Feshbach projection technique applied to this case reveals the existence of a nonlocal three-body operator that cannot be reduced to the sum of two-body n-A and p-A optical potentials plus the n − p interaction V np that binds the incident deuteron. Although this operator does contain terms which include the interactions of, for example, n with the nucleons in A, these terms also depend on the coordinates of p. In other words, the two-body optical potentials that describe N-A elastic scattering are not expected to be the same as the two-body N-A potentials that usually appear in three-body models of a many-body system. In addition, there are new terms that describe physical processes in which the target is excited by the incident neutron and de-excited by the proton and so on in higher orders, thus creating a complicated energy-dependent induced threebody contributions. The necessity of these contributions has been pointed out in two recent Faddeev based calculations. One of these suggests a way of accounting for the energy-dependence of the nucleon optical potentials by modifying the corresponding two-body T -matrices that occur naturally in the Faddeev formalism. The other includes the contribution to the induced three-body force from the first excited state of the target A. It is only within the ADWA framework that an estimate of the contribution to the leading-order nonlocal induced three-body force from all excited target states been reported. However, non-adiabatic contributions are very important for the case of nonlocal potentials implying that further work on many-body treatment of the induced three-body force is needed. Looking forward, we emphasize that future developments should concentrate on understanding the physics missing in the standard DWBA and/or three-body descriptions of deuteron stripping and pick-up reactions. We summarize them as follows.
• A detailed investigation of the remnant ∑ i∈A V ip − U p (R p ) term in the transition operator, based on a modern many-body theory, should be carried out within the DWBA and other reaction models.
• An alternative formalism for the (d, p) reaction in which the remnant term is excluded by an appropriate choice of the exit channel wave function should be extended to include the many-body nature of the nuclei involved.
• To be consistent with modern nuclear Hamiltonians, the contribution to the (d, p) amplitude from realistic 3N forces should be included.
• More research is needed to correctly incorporate induced many-body effects into the three-body dynamical processes that are an important feature of deuteron stripping and pick-up.
• There is also need for further clarification of the role of antisymmetrisation in (d, p) reactions, in particular proton exchange and the Pauli principle, that goes beyond simply multiplying the cross section by a numerical factor to account for neutron identity.
Insights from modern ab-initio approaches to the nuclear reaction and structure theories will play an important role in these developments. They will involve new and unusual matrix elements, that may not be simply related to other observables, and it is there where the help from many-body structure theories, including ab-initio ones, could be sought. The nuclear structure community has learned how to make the many-body problem more manageable by softening NN potentials and constructing effective interactions, for example by using the renormalization group technique. This suppresses the role of high NN momenta but creates important contributions from induced 3N forces. Similar effects will appear in nuclear reaction theory and an understanding of their role is needed. An improved future theory of (d, p) and (p, d) reactions populating isolated nuclear states will also be useful for understanding other reactions such as (d, pγ ) and inclusive reactions, which are of great interest as a surrogate for (n, γ ) cross sections. New developments will also be relevant to proton transfer in (d, n) and (n, d) reactions. Polarization observables and transfer to astrophysically important states with small spectroscopic factors are other areas that will benefit. The developments outlined above are a great challenge for reaction theory and addressing them is important for the correct interpretation of deuteron stripping and pick-up experiments in terms of nuclear structure.