The R-matrix theory in nuclear and atomic physics

From Scholarpedia
Daniel Baye and Pierre Descouvemont (2013), Scholarpedia, 8(1):12360. doi:10.4249/scholarpedia.12360 revision #129959 [link to/cite this article]
Jump to: navigation, search
Post-publication activity

Curator: Daniel Baye

Initially the $R$-matrix theory was aimed at describing resonances in nuclear reactions. At present, the main aim of the $R$-matrix theory is to describe the scattering and reactions resulting from the interaction of particles or systems of particles, which can be nucleons, nuclei, electrons, atoms, or molecules.

Contents

Principle

Cross sections in atomic and nuclear physics may present rapid variations that are known as resonances. The shapes of these resonances often differ strongly and a parametrization of the cross sections in their vicinity as a function of energy and scattering angle is impossible without understanding the underlying physics. An elegant and rather simple solution to this problem was obtained with the $R$-matrix theory. Unexpected consequences of its development were a technique for parametrizing the cross sections of various reactions valid even far from resonances and also, much more surprisingly, an efficient technique for solving coupled-channel Schrödinger equations in the continuum.

The basic idea was introduced by Kapur and Peierls (1938): the configuration space is divided into two regions. The external region allows satisfying scattering properties. The internal region simplifies the treatment of the wave functions at short distances by the use of a square-integrable basis. The basis proposed by Kapur and Peierls however suffered from a complicated energy dependence. The $R$-matrix theory was introduced by Wigner and Eisenbud (1947) according to the same principle, but with a crucial simplification, the use of an energy-independent basis. The advantage of this simplification is that the pivotal quantity in the theory, the $R$-matrix, only involves real energy-independent parameters.

As mentioned before, the $R$-matrix principle relies on a division of the configuration space into the internal and external regions. The boundary between these regions is a parameter known as the channel radius. This radius is chosen large enough so that, in the external region, the different parts of the studied system interact only through known long-range forces and antisymmetrization effects due to the identity of some particles can be neglected. The scattering wave function is approximated in the external region by its asymptotic expression which is known except for coefficients related to the scattering matrix. In the internal region, the system is considered as confined. Its eigenstates thus form a discrete basis. A well chosen square-integrable basis can provide accurate approximations of scattering wave functions over the internal region. Then, the $R$ matrix, which is the inverse of the logarithmic derivative of the wave function at the boundary, can be calculated. A matching with the solution in the external region provides the scattering matrix. This method can also provide the bound states of the system. In this case, the external solution behaves as a decreasing exponential. Since the exponential decrease depends on the unknown binding energy, an iteration is then necessary.

Double evolution of the theory

The $R$-matrix theory was developed into two different directions with little exchange between these variants. Many of its practitioners often ignore the progresses about the other aspect of this double-faced method.

As already mentioned, the original goal was to provide an efficient theory for the treatment of nuclear resonances (Wigner and Eisenbud, 1947; Lane and Thomas, 1958). From information on bound states and low-energy resonances, it soon became clear that the $R$-matrix theory offers an efficient way for accurately parametrizing not only resonances but also the non-resonant part of low-energy cross sections with a small number of parameters (Lane and Thomas, 1958). An important advantage is that most of these parameters have a physical meaning. This first variant of the method is still very important and much employed, in particular to parametrize the low-energy cross sections relevant in nuclear astrophysics. Such parametrizations are useful to provide cross sections at angles and energies where they have not been measured (at least at energies below the maximum energy where measurements took place). In the context of nuclear astrophysics, it is important to reliably extrapolate measured cross sections to the very low energies encountered in stars, at which direct measurements are in general impossible. This version of the $R$-matrix theory will be called hereafter the phenomenological $R$ matrix. Its properties are reviewed in (Lane and Thomas, 1958; Breit, 1959; Descouvemont and Baye 2010).

The other aspect of the $R$-matrix theory is that it can provide a simple and elegant way for solving the Schrödinger equation. It is especially competitive in coupled-channel problems with large numbers of open channels where a direct numerical integration may become unstable. Moreover, the influence of closed channels can be taken into account in a simple way. An additional advantage is that narrow resonances which can escape a purely numerical treatment are easily studied. This other facet of the $R$-matrix theory has been mostly developed in atomic physics although it can also be very useful for nuclear-physics applications. This variant will be called hereafter the calculable or computational $R$ matrix. Its properties are reviewed in (Burke and Robb, 1975; Barrett et al., 1983; Descouvemont and Baye 2010; Burke, 2011).

A very comprehensive review of the phenomenological $R$-matrix method has been given in 1958 by Lane and Thomas. Their article contains most of the important aspects of the phenomenological applications of the $R$ matrix to nuclear physics. Many of their results can also be useful for the calculable $R$ matrix. However, in 1957, just before that review appeared in print, an important improvement of the method was published, which is, therefore, not used in their review. Bloch introduced a singular operator defined on the boundary between the two regions, now known as the Bloch operator, which allows a more elegant and compact presentation of the method (Bloch, 1957). The main interest of the Bloch operator is however that its use led to more general treatments of the resolution of the Schrödinger equation in the internal region and opened the way to accurate methods of resolution in atomic and nuclear physics. Several reviews on the computational $R$ matrix have been published in the context of nuclear physics (Barrett et al., 1973) and of atomic physics (Burke and Robb, 1975; Burke and Berrington (eds), 1993; Burke, 2011). References (Barrett et al., 1983; Descouvemont and Baye 2010) deal with both aspects.

The phenomenological $R$-matrix

Definition

The $R$ matrix allows parametrizing various physical processes and its determination provides collision matrices and cross sections. For each set of good quantum numbers, i.e. total angular momentum and parity, the dimension of the phenomenological $R$ matrix is equal to the number of channels relevant to the physical properties. When a single channel is considered, the $R$ matrix for a partial wave with orbital momentum $l$ (for simplicity, we consider zero spins) is a function of the energy $E$ parametrized by the formula \[\tag{1} R_{l}(E) = \sum_{n=1}^N \frac{\gamma_{nl}^2} {E_{nl} - E} \] where the $\gamma_{nl}^2$ are the reduced widths (Lane and Thomas, 1958). In principle, the $R$ matrix possesses an infinity of poles at the real energies $E_{nl}$ but only a limited number $N$ of such poles affect the low-energy cross sections. The lowest poles are closely related to bound states at negative energies or to narrow resonances at positive energies. Nevertheless, the poles and the energies of physical states are slightly different (see Eq. (14) below for an explanation). Because of this shift, the determination of these parameters from data requires some skill. The real parameters $\gamma_{nl}$ are known as the reduced width amplitudes because their square is a crucial factor of the width of non-overlapping resonances (see Eq. (13) below). Contrary to the reduced width, the width is also very sensitive to the effects of transmission through the Coulomb barrier. The parameters $E_{nl}$ and $\gamma_{nl}$, and thus the $R$ matrix, depend on the chosen channel radius $a$.

The $R$ matrix allows calculating the phase shift $\delta_l$ of the $l$th partial wave at an energy $E$ corresponding to the wavenumber $k$ and velocity $v$, \[\tag{2} \tan \delta_{l} = - \frac{F_l(\eta,ka) - ka R_{l}(E) F_l'(\eta,ka)} {G_l(\eta,ka) - ka R_{l}(E) G_l'(\eta,ka)} \] where $F_l$ and $G_l$ are the regular and irregular Coulomb functions and $\eta = Z_1Z_2e^2/\hbar v$ is the dimensionless Sommerfeld parameter. Low-energy elastic cross sections require fitting all phase shits up to some value $l_{\rm max}$ depending on the maximum energy considered. Expression (2) is often written with an additional parameter, the boundary parameter $B$. As shown by Barker, 1972, the results are independent of the choice of $B$ and the simple value $B=0$ can be chosen without loss of generality.

When studying a reaction, it is necessary to consider more than one channel. In this case the good quantum numbers are the total angular momentum $J$ and parity $\pi$. The spins of the particles in the different channels play a role. The number $N_c$ of channels on which the $R$ matrix depends is in general larger than the number of open physical channels. At energy $E$, the multichannel $R$ matrix is defined as \[\tag{3} R_{c c'}(E) = \sum_{n=1}^N \frac{\gamma_{nc}\gamma_{nc'}}{E_n-E} \] where $c$, $c'$ are channel indices varying from $1$ to $N_c$. Typically, notation $c$ combines the definition of the reaction channels (nature and spins of the nuclei) and the relative orbital momentum. For given $J^\pi$, the $R$ matrix depends on $N$ real poles $E_n$ and $N N_c$ reduced-width amplitudes $\gamma_{nc}$ of pole $E_n$ in channel $c$. The total number of parameters is thus $N (N_c +1)$. It is more difficult to adjust such a number of parameters to available data, specially when several partial waves contribute. For this reason, the number of channels considered in practice is often restricted to $N_c = 2$.

The unitary collision matrix is obtained with \[\tag{4}\newcommand{\ve}[1]{\boldsymbol{#1}} \ve{U}=\ve{Z}^{-1} \ve{Z}^*. \] An element of matrix $\ve{Z}$ reads \[\tag{5} Z_{c c'} = (k_{c'} a)^{-1/2} \Bigl[ O_{c} (\eta_c, k_{c} a)\delta_{c c '} - k_{c'} a R_{c c'} O'_{c'} (\eta_c, k_{c'} a) \Bigr], \] where $k_c = \sqrt{2\mu_c (E-E_c)/\hbar^2}$ is the wavenumber in open channel $c$ with reduced mass $\mu_c$ and threshold energy $E_c < E$, $\eta_c$ is the corresponding Sommerfeld parameter and $O_c (\eta_c,x) = G_{l_c} (\eta_c,x) + i F_{l_c} (\eta_c,x)$ is an outgoing Coulomb wave with orbital momentum $l_c$. The size of the collision matrix depends on the number of open channels. It may be smaller than the size of the $R$ matrix, which may take into account the influence of closed channels.

For complex potentials as encountered in the optical model, expression (4) must be modified into \[\tag{6} \ve{U}=\ve{Z}_O^{-1} \ve{Z}_I, \] where $\ve{Z}_O = \ve{Z}$ and $\ve{Z}_I \neq \ve{Z}_O^*$ is given by a similar expression with outgoing functions replaced by incoming ones $I_c (\eta_c,x) = G_{l_c} (\eta_c,x) - i F_{l_c} (\eta_c,x)$ (Descouvemont and Baye 2010; Druet et al., 2010). The collision matrix $\ve{U}$ is then not unitary.

Resonances

Resonances arise at energies close to a pole $E_{nl}$ of the $R$ matrix. The $R$ matrix is approximated by a single term \[\tag{7} R_l(E) \approx \frac{\gamma_{nl}^2}{E_{nl} - E}. \] The phase shift becomes \[\tag{8} \delta_l(E) \approx \phi_l(E) + \arctan \frac{\gamma_{nl}^2 P_l (E)} {E_{nl} - \gamma_{nl}^2 S_l (E) - E} \] where \[\tag{9} \phi_l(E) = -\arctan [F_l(\eta,ka)/G_l (\eta,ka)] \] is the hard-sphere phase shift, \[\tag{10} P_l (E) = \frac{ka}{F_l(\eta,ka)^2 + G_l (\eta,ka)^2} \] is the penetration factor and \[\tag{11} S_l (E) = P_l (E) [F_l(\eta,ka) F_l'(\eta,ka) + G_l (\eta,ka) G_l'(\eta,ka)] \] is the shift factor.

Expression (8) is close to the Breit-Wigner form of a resonant phase shift \[\tag{12}\newcommand{\dem}{\frac{1}{2}} \delta_l^{\rm BW}(E) \approx \phi_l(E) + \arctan \frac{\dem \Gamma(E)}{E_R - E} \] where the so-called formal width is defined by \[\tag{13} \Gamma(E) = 2 \gamma_{nl}^2 P_l (E). \] While the reduced width and penetration factor depend on $a$, the width is a physical quantity and may thus not depend on the channel radius. It is an energy-dependent quantity whose asymmetric shape depends on the behaviour of $P_l$. Roughly, one can interpret $\gamma_{nl}^2$ as the internal (or nuclear) component of the width and $P_l$ as its Coulomb component. The resonance energy $E_R$ is given by \[\tag{14} E_R = E_{nl} - \gamma_{nl}^2 S_l (E_R) \] and is thus shifted with respect to the pole energy. The factor $S_l (E_R)$ in (14) is calculated at the resonance energy $E_R$. Hence, $E_R$ is defined by an implicit equation, which can be solved by iteration.

Because of the energy dependence of $S_l$, the value of the formal width (13) cannot be directly compared with experiment. Since the shift factor weakly depends on energy, one can make the Thomas approximation (Lane and Thomas, 1958), i.e. use a Taylor expansion as a function of $E - E_R$ limited at first order, \[\tag{15} S_l (E) \approx S_l (E_R) + (E - E_R) \left( \frac{dS_l}{dE} \right)_{E_R}. \] Equation (8) becomes \[\tag{16} \delta_l(E) \approx \phi_l(E) + \arctan \frac{\gamma_{nl}^2 P_l (E)} {(E_R - E) \left[1 + \gamma_{nl}^2 \left( \frac{dS_l}{dE} \right)_{E_R} \right]}. \] One recovers the standard Breit-Wigner expression (12) with an energy-dependent width called observed width given by \[\tag{17} \Gamma_{\rm obs}(E) = \frac{2 \gamma_{nl}^2 P_l (E)} {1 + \gamma_{nl}^2 \left( \frac{dS_l}{dE} \right)_{E_R}} \] which is directly comparable with experiment. In many cases, the reduced width is small and the difference between formal and observed width is negligible with respect to the experimental error bar (see section XII.$3$ of (Lane and Thomas, 1958) for a discussion). A counterexample can be found in (Delbar et al., 1992).

An example: $^{12}$C + p elastic scattering

As an example, let us consider the elastic scattering of protons by $^{12}$C. Three resonances are known in the energy range covered by the data: $1/2^+$ at $0.424$ MeV, $3/2^-$ at $1.558$ MeV and $5/2^+$ at $1.604$ MeV. Data sets are available at the c.m. angles $\theta=89.1^{\circ}$ and $146.9^{\circ}$. They are fitted simultaneously by using $E_R$ and $\Gamma_R$ of the resonant partial waves as adjustable parameters in the single-pole approximation ($N = 1$) of (1). For other partial waves, the hard-sphere phase shift $\phi_l$ is used.

Table 1: $R$-matrix parameters from a simultaneous fit of $^{12}$C+p scattering data (Meyer et al., 1976) at $\theta=89.1^{\circ}$ and $146.9^{\circ}$. Resonance energies $E_R$ are expressed in MeV and widths $\Gamma_R$ in keV
$J^{\pi}=1/2^+$ $J^{\pi}=3/2^-$ $J^{\pi}=5/2^+$
$E_R$ $\Gamma_R$ $E_R$ $\Gamma_R$ $E_R$ $\Gamma_R$
$a=4$ fm $0.427$ $33.8$ $1.560$ $51.4$ $1.603$ $48.1$
$a=5$ fm $0.427$ $32.9$ $1.559$ $51.4$ $1.604$ $48.1$
$a=6$ fm $0.427$ $30.9$ $1.558$ $51.3$ $1.606$ $47.8$
Exp. (Meyer et al., 1976) $0.424$ $33$ $1.558$ $55$ $1.604$ $50$

The fitted resonance properties are given in Table 1 for different channel radii, $a = 4$, $5$ and $6$ fm (Descouvemont and Baye, 2010). The results are almost independent of $a$. The corresponding cross sections are shown in Figure 1. The $R$-matrix parametrization reproduces the data very well, not only in the vicinity of the resonances, but also between them, where the process is mostly non-resonant. It provides a good approximation of the cross section at any angle for $E$ smaller than, or close to, $2$ MeV. The three channel radii provide fits which are indistinguishable at the scale of the figure. The width of the $1/2^+$ resonance depends a little on the channel radius: it decreases from $33.8$ keV at $a = 4$ fm to $30.9$ keV at $a = 6$ fm. This indicates that experimental widths must sometimes be considered with caution. They may depend on the way they are derived from the data.

Figure 1: $R$-matrix fits of $^{12}$C+p experimental excitation functions at two c.m. angles (Meyer et al., 1976) with the parameters of Table 1 (Descouvemont and Baye, 2010).

Criticisms

A drawback of the phenomenological $R$ matrix is that the poles and widths depend on the choice of channel radius, i.e. on a rather arbitrary value. This aspect of the $R$ matrix has been criticized by a number of authors and has led to further developments of the competing phenomenological $K$ matrix (Humblet, 1990). The $K$ matrix, which provides an alternative formulation of the collision matrix, is also expanded in a series involving an infinity of poles. This approach is based on a delicate treatment of Coulomb functions. In spite of the fact that the $K$ matrix does not contain an arbitrary parameter such as the channel radius, its parametrization is more difficult because its parameters are not necessarily real and may have a less direct physical interpretation.

Applications

The phenomenological $R$ matrix and most of its applications were already exhaustively described fifty years ago in (Lane and Thomas, 1958). Among these applications, let us mention a detailed study of resonances and an extension of the method to the description of electromagnetic processes. All these applications have been made in nuclear physics, i.e. for the scattering of neutrons on nuclei or of nuclei on nuclei with the presence of a repulsive Coulomb barrier. Nevertheless, using the method still revealed a number of difficulties. In a series of papers, Barker and collaborators provided practical solutions to the determination of the $R$-matrix parameters from experimental data (Barker et al., 1968; Barker, 1972; Barker, 1988) and applied this framework to the spectroscopy and reactions of light nuclei (Barker, 1987b; Barker and Ferdous, 1980). They also explained non-intuitive effects such as the Thomas-Ehrman shift (Barker and Ferdous, 1980), ghosts of resonances (Barker et al., 1976) and extended the method to further processes such as radiative-capture reactions (Barker, 1980; Barker, 1987a; Barker and Kajino, 1991) and delayed $\beta$ decay (Barker, 1989). The approach developed by Barker and collaborators has become a standard tool for the analysis of low-energy radiative-capture reactions useful in astrophysics. Progresses in the adjustment of $R$-matrix parameters in the single-channel and multichannel cases have been performed in (Angulo and Descouvemont, 2000; Brune, 2002).

More recently, the phenomenological $R$-matrix method was used to analyze scattering experiments involving radioactive beams. These analyses essentially provide information on resonance properties (spin, energy, width). In (Pellegriti et al., 2008), the $^{18}$Ne+p elastic cross section and the $^{18}$Ne(p,p') inelastic cross section were measured simultaneously, and interpreted within the $R$-matrix formalism. This analysis provided evidence for $^{19}$Na states where core excitations are important. Other recent applications on elastic scattering can be found for example in (deBoer et al., 2012; Jung et al., 2012).

In nuclear astrophysics, the main goal of the $R$-matrix formalism (Azuma et al., 2010) is to fit available data in order to extrapolate them to stellar energies, which are much lower than the Coulomb barrier. At those energies, the cross sections are too small to be measured in the laboratory. The important $^{12}$C($\alpha,\gamma)^{16}$O reaction fixes the $^{12}$C/$^{16}$O ratio in the galaxy. Simultaneous fits of the capture and elastic cross sections, and of the $^{16}$N $\beta$-decay spectrum to $^{12}$C$+\alpha$ (Azuma et al., 1994) provide strong constraints on the astrophysical $S$ factor at 300 keV, the relevant energy for this capture reaction in stellar models. Another important reaction in nuclear astrophysics is $^{18}$F(p,$\alpha)^{15}$O (see (Mountford et al., 2012) and references therein). The $\beta$-delayed spectra of $^{12}$B and $^{12}$N are analyzed in (Hyldegaard et al., 2010) with the $R$-matrix formalism to derive information on unbound states of $^{12}$C.

The calculable $R$-matrix

Origin

The aim of the calculable $R$ matrix is to provide an efficient way of solving the Schrödinger equation both at positive and negative energies. It was proposed in 1965 by Haglund and Robson and applied to a two-channel problem involving square-well potentials (Haglund and Robson, 1965). An expansion over a finite basis was introduced by Buttle (Buttle, 1967). He performed the first realistic application on $^{12}$C + n scattering (Buttle, 1967). He also proposed a correction to the truncation of the $R$ matrix to a finite number of poles, that is now named after him. A more serious problem is a discontinuity of the derivative of the wave function at the boundary between the regions that occurs with the traditional choice of basis states inspired by the original ideas in (Wigner and Eisenbud, 1947). Various solutions to the lack of matching at the boundary have been suggested (see Barrett et al., 1983 for a review). This apparent problem has attracted a lot of attention even long after an efficient technique where it does not occur was introduced (Lane and Robson, 1966; Lane and Robson, 1969). By dropping an unnecessary condition, the $R$-matrix method can be very accurate without matching problems and without need for a Buttle correction (see Choice of basis and misconceptions and Descouvemont and Baye, 2010 for details).

Some users of the phenomenological $R$ matrix consider the channel radius as a parameter, which must be optimized when fitting the data. Even if this dependence on a parameter without strong physical meaning is weak, this is a drawback that would not be acceptable when aiming at accurately solve the Schrödinger equation. Hence a crucial test of the results of the calculable $R$ matrix is an almost perfect independence with respect to the choice of channel radius. This test provides a measure of the accuracy of the calculations.

Calculation of the $R$-matrix for scattering by a potential

After separation of the angular part, the radial Schrödinger equation in partial wave $l$ can be written as \[\tag{18} (H_l - E) u_l = 0. \] In this expression, the radial Hamiltonian $H_l$ is defined as \[\tag{19} H_l = T_l + V = -\frac{\hbar^2}{2\mu} \left( \frac{d^2}{dr^2} - \frac{l(l+1)}{r^2} \right) + V(r), \] where $\mu$ is the reduced mass and potential $V$ is well approximated by the Coulomb potential beyond the channel radius $a$. The bounded radial solutions $u_l$ vanish at the origin and have at positive energies the asymptotic behaviour \[\tag{20}\newcommand{\arrow}[2]{\ _\overrightarrow{#1 \rightarrow #2}\ } u_{l} (r) \arrow{r}{\infty} \cos \delta_l\, F_l(\eta,kr) + \sin \delta_l\,G_l (\eta,kr). \]

An important progress of the calculable $R$-matrix method was the introduction of the Bloch operator (Bloch, 1957) \[\tag{21}\newcommand{\cL}{\cal L} {\cL} = \frac{\hbar^2}{2\mu}\, \delta (r-a) \left( \frac{d}{dr} - \frac{B}{r} \right). \] The boundary parameter $B$ can be chosen arbitrarily and the physical results are independent of its value (Descouvemont and Baye, 2010). The wave function $\newcommand{\uint}{u_{l}^{\rm int}}\uint$ in the internal region is approximated by solutions of the inhomogeneous Bloch-Schrödinger equation \[\tag{22}\newcommand{\uext}{u_{l}^{\rm ext}} (H_l + {\cL} - E) \uint = {\cL} \uext \] where the external solution $\uext$ approximated by (20) is used in the right-hand member. This equation is complemented by the continuity condition \[\tag{23} \uint (a) = \uext (a) \] at the boundary. The Bloch operator has two important advantages: (i) it corrects the non-Hermiticity of the Hamiltonian $H_l$ over the internal region and (ii) it matches the logarithmic derivatives of $\uint$ and $\uext$ at $r=a$, thus providing a smooth connection between the two parts of the radial function.

In the internal region, the wave function $\uint$ is expanded over some finite basis involving $N$ linearly independent square integrable functions $\varphi_i$ as \[\tag{24} \uint (r) = \sum_{i=1}^N c_i \varphi_i (r). \] The functions $\varphi_i$ vanish at the origin but are not necessarily orthogonal. At $r = a$, they are not assumed to satisfy a specific boundary condition. (Many papers impose such a condition, but it has unfavourable effects on the convergence; see Choice of basis and misconceptions and Descouvemont and Baye (2010) for a discussion.)

For $B = 0$, the $R$ matrix at energy $E$ is defined by \[\tag{25} u_l (a) = R_l(E) a u_l' (a). \] The inverse of the $R$ matrix is thus the dimensionless logarithmic derivative of the radial wave function at the boundary between both regions. This "matrix" has dimension $1$ in a single-channel case and is just a function of energy. It also depends on the channel radius. Expansion (24) is introduced in (22) and the resulting equation is projected on $\varphi_i (r)$. The calculable $R$ matrix is then obtained as \[\tag{26} R_l (E) = \frac{\hbar^2}{2\mu a} \sum_{i,j=1}^N \varphi_i(a) (\ve{C}^{-1})_{ij} \varphi_j(a), \] where the elements of the symmetric matrix $\ve{C}$ are defined as \[\tag{27}\newcommand{\la}{\langle}\newcommand{\ra}{\rangle} C_{ij}(E) = \la \varphi_i | T_l + {\cL} + V - E | \varphi_j \ra. \] Dirac brackets correspond here to one-dimensional integrals over the variable $r$ from $0$ to $a$.

When the basis functions $\varphi_i (r)$ are orthonormal, the eigenvalues $E_{nl}$ and the corresponding normalized eigenvectors $\ve{v}_{nl}$ of matrix $\ve{C}(0)$ lead to the spectral decomposition \[\tag{28} [\ve{C}(E)]^{-1} = \sum_{n=1}^{N} \frac{\ve{v}_{nl} \ve{v}_{nl}^{\rm T}}{E_{nl} - E} \] where T means transposition. The $R$ matrix (26) can be written as in Eq. (1), \[\tag{29} R_l(E) = \sum_{n=1}^N \frac{\gamma_{nl}^2}{E_{nl} - E}. \] In this expression, the reduced width amplitudes are given by \[\tag{30} \gamma_{nl} = \left( \frac{\hbar^2}{2\mu a} \right)^{1/2} \phi_{nl} (a) \] with \[\tag{31} \phi_{nl} (r) = \sum_{i=1}^N v_{nl,i} \varphi_i(r), \] where $v_{nl,i}$ is the $i$th component of $\ve{v}_{nl}$. The reduced width amplitudes are proportional to the value at the channel radius of variational approximations $\phi_{nl}$ of the eigenfunctions of the Hermitian operator $H_l + {\cL}$. The functions $\phi_{nl}$ corresponding to the lowest energies $E_{nl}$ thus represent approximate eigenfunctions of the physical problem confined over the interval $(0,a)$ with a fixed logarithmic derivative at $r = a$. The other functions, though unphysical, are important to ensure a smooth matching between the internal and external wave functions (Baye et al., 1998).

The traditional expression for the theoretical $R$ matrix is obtained when $N$ tends towards infinity in a complete basis as \[\tag{32} R_l(E) = \sum_{n=1}^\infty \frac{\gamma_{nl}^2}{E_{nl} - E}. \] The energies $E_{nl}$ are now the exact eigenvalues of the operator $H_l + {\cL}$ and the reduced width amplitudes $\gamma_{nl}$ are related to the values at $r=a$ of its exact eigenfunctions. The $R$ matrix is a real function when $V$ is real. It has an infinity of real simple poles, bounded from below. Its derivative is always positive at regular points.

The $R$-matrix theory can be extended to several channels and to non-local potentials. A partial wave $JM\pi$ of the total wave function of the collision at energy $E$ can be written as \[\tag{33} \Psi_{(c_0)}^{JM\pi} = \sum_{c} |c \ra r^{-1} u_{c(c_0)} (r), \] where $|c \ra$ is a channel state describing the internal properties of the colliding nuclei and the angular part of their relative motion (Descouvemont and Baye, 2010), and index $c_0$ refers to the entrance channel. Introducing (33) in the Schrödinger equation and projecting over $|c\ra$ leads to the coupled equations \[\tag{34} \left(T_c + E_{c} - E \right) u_{c(c_0)}(r) + \sum_{c'} V_{cc'}(r) u_{c'(c_0)}(r) = 0 \] for $c = 1$ to $N_c$. In this expression, $E_c$ is the threshold energy of channel $c$, \[\tag{35} T_c = -\frac{\hbar^2}{2\mu_c} \left( \frac{d^2}{dr^2} - \frac{l_c(l_c+1)}{r^2} \right) \] is the kinetic energy operator for a channel with reduced mass $\mu_c$ and orbital momentum $l_c$, and $V_{cc'}$ is an element of a potential matrix with the property $V_{cc'} (r) \rightarrow Z_{1c}Z_{2c}e^2 r^{-1} \delta_{cc'}$ when $r$ tends to infinity. The potential matrix could be non local without much additional complication. The asymptotic form of the radial wave functions $u_{c(c_0)}(r)$ and the radial wave functions in the external region at energy $E$ are defined as \[\tag{36} u^{\rm ext}_{c(c_0)}(r)= \left\{\begin{array}{ll} v_{c}^{-1/2} \Bigl( I_{c} (k_{c} r)\delta_{c c_0} - U_{c c_0} O_{c} (k_{c} r) \Bigr) & {\rm \ for\ } E > E_c \\ A_{cc_0} W_c(2\kappa_{c}r) & {\rm \ for\ } E < E_c, \end{array} \right . \] where $c_0$ is the entrance channel ($E_{c_0} < E$) and $W_c(x) \equiv W_{-\eta_{c},l_c+\frac{1}{2}}(x)$ is a Whittaker function. Here closed channels must be taken into account. In these cases, the real parameters $\kappa_c$ and $\eta_c$ are calculated with the positive energy difference $E_c - E$. The Bloch operator is defined in the multichannel formalism as \[\tag{37} {\cL} = \sum_{c} |c\ra {\cL}_{c} \la c |,\ \ \ \ \ \ \ \ \ \ \ {\cL}_{c}= \frac{\hbar^2}{2\mu_{c}} \, \delta (r-a) \left( \frac{d}{dr} - \frac{B_{c}}{r} \right), \] where coefficients $B_{c}$ are chosen as zero for open channels and as \[\tag{38} B_c = 2\kappa_c a \frac{W_c'(2\kappa_c a)}{W_c(2\kappa_c a)} \] for closed channels because this choice suppresses the right-hand side of the Bloch-Schrödinger system of equations (39) in channel $c$. Notice that these coefficients then depend on energy for closed channels. The Bloch-Schrödinger system of equations is given by \[\tag{39} \sum_{c'} \Bigl[(T_{c}+{\cL}_{c}+E_c-E)\delta_{c c'} + V_{c c'} \Bigr] u^{\rm int}_{c'}(r) = {\cL}_{c} u^{\rm ext}_{c}. \] The internal parts of the radial wave functions are expanded over a basis $\varphi_i(r)$, \[\tag{40} u^{\rm int}_c (r) = \sum_{i=1}^N c_{c i} \varphi_i(r). \] The $R$ matrix is symmetric, with elements given by \[\tag{41} R_{c c'}(E) = \frac{\hbar^2}{2\sqrt{\mu_{c}\mu_{c'}} a} \sum_{i,i'=1}^N \varphi_i(a) (\ve{C}^{-1})_{c i, c'i'} \varphi_{i'}(a), \] where \[\tag{42} C_{c i, c'i'} = \la \varphi_{i}| T_{c}+{\cL}_{c}+E_c-E |\varphi_{i'} \ra \delta_{c c'} + \la \varphi_{i}| V_{c c'} |\varphi_{i'} \ra. \] The spectral decomposition of the symmetric matrix $\ve{C}$ for an orthonormal basis provides the canonical form (3) of the multichannel $R$ matrix, where the real poles $E_n$ are the eigenvalues of $\ve{C}$ and the reduced-width amplitude of pole $E_n$ in channel $c$ is expressed as a function of the components of the corresponding normed eigenvector $\ve{v}_n$ as \[\tag{43} \gamma_{nc} = \left( \frac{\hbar^2}{2\mu_{c} a} \right)^{1/2} \sum_{i=1}^N v_{n,ci} \varphi_i(a). \] The size of calculable $R$ matrices is often much larger than the size of phenomenological $R$ matrices. Indeed, the size of the calculable $R$ matrix is determined by a physical choice, the number of configurations included in the calculation. It is only limited by the available computation times. The size of the phenomenological $R$ matrix is limited by the availability of experimental data. Such data are rarely available for more than two channels.

Choice of basis and misconceptions

Considerable confusion exists in the literature about the properties that basis states $\varphi_i$ should have. Improper choices have led to the introduction of corrections and to attempts to use the boundary parameter $B$ to correct drawbacks of the basis. However single-channel results are independent of the choice of $B$. This confusion has sometimes led to an undeserved reputation of poor convergence for the calculable $R$-matrix method.

In their seminal paper, Wigner and Eisenbud wanted to provide a phenomenological description of resonances (Wigner and Eisenbud, 1947). They did not intend to propose a technique of resolution of the Schrödinger equation. To reach their goal they assume that the basis functions all satisfy the boundary conditions $\varphi_i(0)=0$ and \[\tag{44} a\varphi'_i(a) -B \varphi_i(a) = 0. \] This procedure leads to $R$ matrix (32). When used as a technique of resolution, the finite-basis $R$ matrix (26) or (29) obtained with this procedure does not converge uniformly. The reason is simple. The first derivative of the wave function suffers from a discontinuity at $r=a$ (Szmytkowski and Hinze, 1996). The limit of ${\uint}'$ when $r$ tends towards $a$ to the left is equal to ${\uext}'(a)$ but not to ${\uint}'(a)$, \[\tag{45} \lim_{r \rightarrow a^-} {\uint}'(r) = {\uext}'(a) \ne {\uint}'(a). \] For example, if $B = 0$, $\varphi'_i (a)$ vanishes for all $i$ values and one readily sees that ${\uint}'(a) = 0$ at all energies. This property has unfavourable consequences on the convergence of numerical methods when the basis is truncated since the logarithmic derivative of the external solution depends on the phase shift (and thus on energy) and can not be matched with the internal solution (see Figure 3). Buttle (Buttle, 1967) has proposed a correction to the $R$-matrix truncation. Although this correction improves the phase shifts, it does not really solve the problem because it does not improve the wave functions.

This problem received a solution with the works of Lane and Robson (Lane and Robson, 1966; Lane and Robson, 1969; Philpott, 1975). Their method was successfully applied in nuclear physics where traditional basis functions do not satisfy (44) and, on the contrary, display a variety of behaviours at the channel radius. With oscillator basis functions, accurate results for neutron-nucleus scattering could be obtained (Philpott, 1975; Philpott and George, 1974). At the same time, a microscopic extension of the $R$ matrix using a fully antisymmetrized two-centre harmonic-oscillator model provided accurate phase shifts for collisions between light nuclei with few basis states (Baye and Heenen, 1974; Baye et al., 1977). The success of these calculations relies on the fact that the Bloch operator makes condition (44) unnecessary. Since the results do not depend on $B$, the choice $B=0$ can be used. A general though economical method for solving coupled-channel problems is described in (Hesse et al., 1998).

The negative role of condition (44) remained long unnoticed in atomic physics where in many cases basis states were imposed to satisfy (44). To illustrate this negative role, let us compare two families of basis functions, satisfying or not this condition. Since many types of basis functions exist and are used in practical calculations, the present choice is specifically made to emphasize the drawbacks of condition (44). Sine basis functions are given for $i = 1$ to $N$ by \[\tag{46} \varphi_i(r)=\sin \left[ \frac{\pi r}{a}(i-1/2) \right]. \] Their derivatives vanish at $r=a$ and thus satisfy the unnecessary condition (44) with $B = 0$. The Legendre basis is defined for $i = 1$ to $N$ by \[\tag{47} \varphi_i(r)=r P_{i-1}(2r/a-1), \] where $P_n$ is the Legendre polynomial of degree $n$. The factor $r$ ensures that the wave function vanishes at the origin. The logarithmic derivatives of these basis functions present a variety of behaviours at $r=a$. In fact, it is convenient to replace these Legendre polynomials by the equivalent Lagrange-Legendre basis (Hesse et al., 2002) \[\tag{48} \varphi_i (r)=(-1)^{N+i} \sqrt{\frac{1-x_i}{a x_i}}\, \frac{r P_N(2r/a-1)}{r-a x_i}, \] where the $x_i$ are the $N$ zeros of \[\tag{49} P_N(2x_i-1)=0. \] This basis is equivalent to the Legendre basis since it also involves $N$ linearly independent polynomials of degree at most $N$ vanishing at the origin.

The Lagrange basis is particularly useful when an additional approximation, which does essentially not cost any loss of accuracy, is made (Baye et al., 2002) (see Appendix). The $R$-matrix method on a Lagrange mesh (Malegat, 1994; Baye et al., 1998; Hesse et al., 1998; Hesse et al., 2002) is accurate and economical. It is very convenient in calculations where the matrix elements of the potential, local or non local, are long to compute such as in ab initio calculations (Quaglioni and Navrátil, 2008).

An example: $^{12}$C + p scattering

As an example, let us again consider the elastic scattering of protons by $^{12}$C. As seen in the preceding example, the $^{12}$C+p system presents a narrow $1/2^+$ resonance at low energies (see Table 1). The nuclear and Coulomb potentials \[\tag{50} V_N(r)=-73.8 \exp[-(r/2.70)^2],\\ V_C(r)=6e^2/r, \] reproduce the resonant behaviour of the $s$ phase shifts near 0.42 MeV. The units in the potential are fm and MeV for lengths and energies, respectively ($e^2 = 1.44$ MeV fm).

Figure 2: $^{12}$C+p $R$-matrix phase shifts (in degrees) as a function of the energy $E$ for different bases and conditions ($l=0$): (a) Lagrange functions at $a=8$ fm (the exact results are superimposed to the $N=10$ curve), (b) sine functions at $a=8$ fm, (c) Lagrange functions for $N=15$ (the $a=7$ fm curve corresponds to the exact results) (Descouvemont and Baye, 2010).

Exact phase shifts from numerical resolutions of the Schrödinger equation are compared with $R$-matrix calculations using Lagrange (a) and sine (b) functions in Figure 2. The channel radius is chosen as $a=8$ fm, where the nuclear interaction is negligible. Figure 2 (a) shows that with the Lagrange functions (48), the convergence is reached with $N \geq 10$. The value $N=7$ is an example of a number of points not large enough. Figure 2 (b) illustrates the sine functions (46), poorly adapted to a good matching since the left derivative of the wave function at $r=a$ is zero. Even $N=20$ is far from the exact calculation. Figure 2 (c) presents the convergence as a function of the channel radius with the Lagrange functions ($N=15$ is fixed). As could be expected, $a=5$ fm is too small ($\left|V_N(a)/V_C(a)\right|=1.4$). To obtain a good stability, channel radii larger than 6 fm should be used.

The matching problem is illustrated in Figure 3, where is shown the wave function at $E=2$ MeV with $a=8$ fm and $N=15$. With the Lagrange functions, the matching between the internal and external wave functions is quite smooth. For sine functions, the matching is poor, which has a direct impact on the phase shift.

Figure 3: $^{12}$C+p $l=0$ wave functions with Lagrange and sine functions at $E=2$ MeV ($a=8$ fm, $N=15$). Solid lines represent the internal wave function, and dotted lines the external ones. The Lagrange wave function is superimposed to the exact result (Descouvemont and Baye, 2010).

Applications

After its introduction for nuclear-physics problems (Haglund and Robson, 1965; Buttle, 1967), the $R$-matrix method was first extensively developed to study electron (or positron) collisions on atoms and molecules (Burke and Robb, 1975; Burke and Berrington (eds), 1993; Burke, 2011). Important and difficult aspects of these atomic-physics problems are the non-locality of the interaction due to electron exchanges and the long-range nature of the interactions, due to polarization effects. The non-locality is well treated in the $R$-matrix approach. The long range of the force implies that the asymptotic behaviour of the solution is only reached for very large values of the interparticle distance. To avoid using a very large channel radius, propagation methods have been introduced (Light and Walker, 1976). They involve an intermediate region where the interaction can be simplified, for example with an asymptotic expansion. The applications of the $R$ matrix in atomic physics are reviewed in (Burke, 2011).

In nuclear physics, the computational $R$ matrix is less used. The microscopic $R$-matrix method is much applied in microscopic cluster calculations in which the $A$-body Schrödinger equation is solved in an approximate way by assuming that clusters exist. These clusters are substructures such as an $\alpha$ particle for which a frozen internal wave function is employed in the model. The difficult antisymmetrization is taken into account exactly, but in the internal region only (Baye and Heenen, 1974; Baye et al., 1977). Many microscopic calculations of elastic and inelastic scattering can be performed. The versatility of the $R$ matrix finds interesting applications in processes where bound and scattering states are mixed, such as radiative capture or delayed $\beta$ decay (Baye and Descouvemont, 1983; Baye and Descouvemont, 1987b).

The $R$ matrix can be very useful in large coupled-channel calculations. A significant simplification occurs when the Lagrange functions (48) are chosen as basis functions and the consistent Gauss quadrature is used as explained in the Appendix (Baye et al., 1998; Hesse et al., 1998; Descouvemont and Baye, 2010). Coupled-channel calculations are simplified because calculations of matrix elements of the potentials are avoided (Malegat, 1994; Baye et al., 1998; Hesse et al., 1998). This approach has been extended to non-local interactions (Hesse et al., 2002). It is applied in recent ab initio calculations (Quaglioni and Navrátil, 2008).

The $R$ matrix has recently been applied to the Continuum Discretized Coupled Channel (CDCC) method. The CDCC method was suggested by Rawitscher (Rawitscher, 1974) and developed and used by several groups (Austern et al., 1987; Nunes and Thompson, 1999). The purpose of the CDCC method is to determine, as accurately as possible, the scattering and dissociation cross sections of a nucleus which can be easily broken up in the nuclear and Coulomb fields of a target. The final states may thus involve three or more particles: the target and the fragments of the projectile. The principle of the CDCC method is to replace the continuum states describing the relative motion of the fragments of the projectile by square integrable approximations at discrete energies. The CDCC equations then take the form of standard coupled-channel equations. The $R$ matrix can be useful at two different levels (Descouvemont and Baye, 2010; Druet et al., 2010). On one hand, it can be used to describe the discretized continuum states in a simple way by averaging internal wave functions over energy in the form of bins. On the other hand, as mentioned above, it can be used to solve the coupled-channel equations in an efficient way. Here propagation methods (Light and Walker, 1976; Descouvemont and Baye, 2010) can be useful.

Another recent application concerns three-body scattering. When three unbound particles interact, such as in the final state of a breakup process, the scattering is described by three-body collision matrices. These matrices are infinite because of the infinity of ways the three particles can share the total angular momentum. A truncation is thus necessary in practice. The $R$ matrix is used as a tool to derive the three-body phase shifts in a practical way and to provide manageable wave functions (Descouvemont et al., 2006; Descouvemont and Baye, 2010). Propagation methods are necessary, given the large value of the channel radius. These wave functions are useful to describe breakup reactions (Baye et al., 2009; Pinilla et al., 2012).

Conclusion and outlook

In this overview, we have presented the two variants of the $R$-matrix method. To date, the phenomenological $R$ matrix remains close to its original spirit and is very much used in nuclear physics to parametrize low-energy cross sections. Its main merits are that all parameters are real and have a physical meaning. Although resonances often play a crucial role in these parametrizations, non-resonant cross sections are accurately described as well. The calculable $R$ matrix is an efficient technique to solve the Schrödinger equation in various situations as well as its relativistic extensions. It underwent many of its developments in atomic physics but becomes more and more useful in nuclear physics.

Most progresses can be expected in the calculable $R$ matrix. A challenge in scattering problems is the occurrence of large coupled-channel systems. This situation is often met in CDCC problems (see for example Druet and Descouvemont (2012)), and in particular when the projectile presents a three-body structure (see for example Matsumoto et al. (2004)). Large coupled-channel systems also show up in three-body continuum calculations (see for example Descouvemont et al. (2006)). The problem of a large number of coupled equations has been addressed in finite-difference methods, where the resolution of the coupled-channel system is replaced by a single inhomogeneous equation (Raynal, 1972; Rhoades-Brown et al., 1980). An iterative procedure provides specific elements of the collision matrix (Thompson, 1988). Although it may not converge for strong coupling potentials, its main advantage is to deal with a single equation. For weak coupling potentials, the convergence can be accelerated by using Padé approximants (Rhoades-Brown et al., 1980). These techniques might be introduced in the $R$-matrix formalism and provide an efficient way to deal with large coupled-channel systems.

Optimizing the propagation technique is another challenge. The propagation of the $R$ matrix is necessary in problems where the channel radius is too large. One or several intermediate regions are introduced between the internal and external regions. The $R$ matrix must be accurately propagated from boundary to boundary. Reducing the computational cost of propagation while keeping the accuracy is an important issue.

Appendix: The Lagrange-mesh method

The Lagrange-mesh method combines the basis (48) with the Gauss quadrature associated with the mesh points defined from (49). The Lagrange-Legendre basis functions (48) satisfy the Lagrange conditions \[\tag{51} \varphi_i(a x_j)=(a \lambda_i)^{-1/2}\delta_{ij}, \] where $\lambda_i$ is the weight of the Gauss-Legendre quadrature corresponding to the $[0,1]$ interval. If the matrix elements with basis functions (48) are computed at the Gauss approximation of order $N$, consistent with the $N$ mesh points $ax_i$, their calculation is strongly simplified. At this approximation, the overlap is given by \[\tag{52} \langle \varphi_i|\varphi_j\rangle =\int_0^a \varphi_i(r)\varphi_j(r) dr\approx \delta_{ij}. \] The local-potential matrix is diagonal with elements simply given by the value of the potential at the mesh points, \[\tag{53} \langle \varphi_i|V|\varphi_j\rangle =\int_0^a \varphi_i(r) V(r) \varphi_j(r) dr \approx V(a x_i)\delta_{ij}. \] This introduces an important simplification since no integrals involving the potential need be calculated. This method is easily extended to non-local potentials (Hesse et al., 2002).

References

Internal references

Further reading

External links

See also

Nuclear astrophysics, Bloch operator, K-matrix, Hamiltonian

Sponsored by: Dr. Nicolas Alamanos, CEA, IRFU; F-91191 Gif-sur-Yvette, France, France
Reviewed by: Prof. Filomena Nunes, National Superconducting Cycltron Laboratory and Department of Physicis and Astronomy Michigan State University
Reviewed by: Prof. Carl R. Brune, Institute for Nuclear and Particle Physics Department of Physics and Astronomy Ohio, University Athens, Ohio, Ohio 45701, USA
Accepted on: 2013-01-21 16:56:40 GMT
Retrieved from "http://www.scholarpedia.org/w/index.php?title=The_R-matrix_theory_in_nuclear_and_atomic_physics&oldid=129959"
Personal tools
Namespaces

Variants
Actions
Navigation
Focal areas
Activity
Tools