The aim of this project is to solve the Lippman-Schwinger equation for two interacting nucleons and relate the obtained phase shifts with those extracted from the experimental cross sections.
We are going to solve the Schroedinger equation (SE) for the neutron-proton system in momentum space for positive energies \( E \) in order to obtain the phase shifts. We can rewrite the SE in momentum space as $$ \begin{equation} \frac{k^2}{m}\psi_l(k)+\frac{2}{\pi}\int_0^{\infty}dqq^2V_l(k,q)\psi_l(q)=E\psi_l(k). \label{eq:sem} \end{equation} $$ The derivation is covered by the lecture notes on scattering theory. This derivation is discussed during the regular lectures.
Here we have used units \( \hbar=c=1 \). This means that \( k \) has dimension energy. \( k \) is the relative momentum between the two particles. A partial wave expansion has been used in order to reduce the problem to an integral over the magnitude of momentum only. The subscript \( l \) refers therefore to a partial wave with a given orbital momentum \( l \). Below we will let the interaction to couple different values of \( l \). This leads to what is called a coupled-channel problem. For details, see again the lecture notes on scattering theory.
To obtain the potential in momentum space we used the Fourier-Bessel transform (Hankel transform) $$ \begin{equation} V_l(k,k')= \int j_l(kr)V(r)j_l(k'r)r^2dr, \label{eq:vtrans} \end{equation} $$ where \( j_l \) is the spherical Bessel function. We will just study the case \( l=0 \), which means that \( j_0(kr)=sin(kr)/kr \).
For scattering states, \( E>0 \), the corresponding equation to solve is the so-called Lippman-Schwinger equation. This is an integral equation where we have to deal with the amplitude \( R(k,k') \) (reaction matrix) defined through the integral equation $$ \begin{equation} R_l(k,k') = V_l(k,k') +\frac{2}{\pi}\hat{P} \int_0^{\infty}dqq^2V_l(k,q)\frac{1}{E-q^2/m}R_l(q,k'), \label{eq:ls1} \end{equation} $$ where the total kinetic energy of the two incoming particles in the center-of-mass system is $$ \begin{equation} E=\frac{k_0^2}{m}. \label{_auto1} \end{equation} $$ The symbol \( \hat{P} \) indicates that Cauchy's principal-value prescription is used in order to avoid the singularity arising from the zero of the denominator. We will discuss below how to solve this problem. Eq. \eqref{eq:ls1} represents then the problem you will have to solve numerically.
The matrix \( R_l(k,k') \) relates to the the phase shifts through its diagonal elements as $$ \begin{equation} R_l(k_0,k_0)=-\frac{tan\delta_l}{mk_0}. \label{eq:shifts} \end{equation} $$
From now on we will drop the subscript \( l \) in all equations.
In order to solve the Lippman-Schwinger equation in momentum space, we need first to write a function which sets up the mesh points. We need to do that since we are going to approximate an integral through $$ \begin{equation*} \int_a^bf(x)dx\approx\sum_{i=1}^Nw_if(x_i), \end{equation*} $$ where we have fixed \( N \) lattice points through the corresponding weights \( w_i \) and points \( x_i \).
The transform of a potential on the form \( V(r)=V_0\exp{(-\mu r)}/(\mu r) \) which we rewrite more generally as $$ V_{\eta}\frac{e^{-\eta x}}{x}, $$ with \( x=\mu r \). Hint. The problem for \( l=0 \) is to evaluate the integral $$ I=V_{\eta}\int_{0}^{\infty}dr r^{2}j_{0}(kr)\frac{e^{-\eta x}}{x}j_{0}(k'r) $$ By using the substitution \( x=\mu r \) and converting the product of sines to a difference of cosines through the relation $$ \sin u\sin v=\frac{1}{2}[\cos(u-v)-\cos(u+v)] $$ we obtain $$ I=\frac{V_{\eta}}{2\mu kk'}\int_{0}^{\infty}dx (\cos\alpha x-\cos\beta x) \frac{e^{-\eta x}}{x}. $$ where \( \alpha=\frac{1}{\mu}(k-k') \) and \( \beta=\frac{1}{\mu}(k+k') \). This integral can be converted to a form found in standard integral tables by adding and subtracting 1: $$ I=\frac{V_{\eta}}{2\mu kk'}\left[\int_{0}^{\infty}dx(1-\cos\beta x) \frac{e^{-\eta x}}{x}-\int_{0}^{\infty}dx(1-\cos\alpha x)\frac{e^{-\eta x}} {x}\right]. $$ Both integrals are now of the same form and their values are given as $$ I=\frac{V_{\eta}}{4\mu k k'}\ln\left[\frac{(\mu \eta)^{2}+(k+k')^{2}} {(\mu\eta)^{2}+(k-k')^{2}}\right]. $$
Write a function which calculates the expressions for the potential in momentum space.
We can then use the trick in Eq. \eqref{eq:trick} to rewrite Eq. \eqref{eq:ls1} as $$ \begin{equation} R(k,k') = V(k,k') +\frac{2}{\pi} \int_0^{\infty}dq \frac{q^2V(k,q)R(q,k')-k_0^2V(k,k_0)R(k_0,k') } {(k_0^2-q^2)/m}. \label{eq:ls2} \end{equation} $$ This is the equation you are going to solve numerically in order to calculate the phase shifts of Eq. \eqref{eq:shifts}.We are interested in obtaining \( R(k_0,k_0) \).
How do we proceed in order to solve Eq. \eqref{eq:ls2}?
The first step consists in using the mesh points \( k_j \) and the weights \( \omega_j \). We can rewrite Eq. \eqref{eq:ls2} as $$ \begin{equation} \label{eq:ls3} R(k,k') = V(k,k') +\frac{2}{\pi}\sum_{j=1}^N\frac{\omega_jk_j^2V(k,k_j)R(k_j,k')}{(k_0^2-k_j^2)/m}-\frac{2}{\pi}k_0^2V(k,k_0)R(k_0,k')\sum_{n=1}^N\frac{\omega_n}{(k_0^2-k_n^2)/m}. \end{equation} $$ This equation contains now the unknowns \( R(k_i,k_j) \) (with dimension \( N\times N \)) and \( R(k_0,k_0) \). We can turn Eq. \eqref{eq:ls3} into an equation with dimension \( (N+1)\times (N+1) \) with a mesh which contains the original mesh points \( k_j \) for \( j=1,N \) and the point which corresponds to the energy \( k_0 \). Consider the latter as the 'observable' point. The mesh points become then \( k_j \) for \( j=1,n \) and \( k_{N+1}=k_0 \).
With these new mesh points we define the matrix $$ \begin{equation} \label{eq:aeq} A_{i,j}=\delta_{i,j}-V(k_i,k_j)u_j, \end{equation} $$ where \( \delta \) is the Kronecker \( \delta \) and $$ \begin{equation} u_j=\frac{2}{\pi} \frac{\omega_jk_j^2}{(k_0^2-k_j^2)/m}\hspace{1cm} j=1,N \label{_auto4} \end{equation} $$ and $$ \begin{equation} u_{N+1}=-\frac{2}{\pi} \sum_{j=1}^N\frac{k_0^2\omega_j}{(k_0^2-k_j^2)/m}. \label{_auto5} \end{equation} $$ The first task is then to set up the matrix \( A \) for a given \( k_0 \). This is an \( (N+1)\times (N+1) \) matrix. It can be convenient to have an outer loop which runs over the chosen observable values for the energy \( k_0^2/m \). {\em Note that all mesh points \( k_j \) for \( j=1,N \) must be different from \( k_0 \). Note also that \( V(k_i,k_j) \) is an \( (N+1)\times (N+1) \) matrix}. Write a small function which sets up \( A \).
With the matrix \( A \) we can rewrite Eq. \eqref{eq:ls3} as a matrix problem of dimension \( (N+1)\times (N+1) \). All matrices \( R \), \( A \) and \( V \) have this dimension and we get $$ \begin{equation} A_{i,l}R_{l,j}=V_{i,j}, \label{_auto6} \end{equation} $$ or just $$ \begin{equation} AR=V. \label{eq:final1} \end{equation} $$ Since you already have defined \( A \) and \( V \) (these are stored as \( (N+1)\times (N+1) \) matrices) Eq. \eqref{eq:final1} involves only the unknown \( R \). We obtain it by matrix inversion, i.e., $$ \begin{equation} \label{eq:final2} R=A^{-1}V. \end{equation} $$ Thus, to obtain \( R \), you will need to set up the matrices \( A \) and \( V \) and invert the matrix \( A \). To do that you can use the examples on matrix inversion (this link brings you to the c++ version). With the inverse \( A^{-1} \), performing a matrix multiplication with \( V \) results in \( R \).
With \( R \) you can then evaluate the phase shifts by noting that $$ \begin{equation} R(k_{N+1},k_{N+1})=R(k_0,k_0), \label{_auto7} \end{equation} $$ and you are done.
You can choose to read \( k_0 \) from file or screen, or set up a loop over chosen values of \( k_0 \) and for each \( k_0 \) solve Eq. \eqref{eq:final2}.
When you have \( R(k,k') \) for the given potential, evaluate now the phase-shifts using $$ \begin{equation*} R(k_0,k_0)=-\frac{tan\delta}{mk_0}. \end{equation*} $$
Compare the phase shifts for the potential of Eq. \eqref{eq:realp} with the experimental phase shifts that can be found in the article of the Nijmegen group in Physical Review C 48, 792 (1993). Alternatively look up their website
Here we explore a simple alternative to the momentum-space matrix inversion for the calculation of scattering phase shifts called the Variable Phase Approach (VPA). The VPA (in its simplest formulation) is not as flexible as the matrix inversion method in that it is limited to local potentials (i.e., \( \langle{\bf r'}|V|{\bf r}\rangle = \delta^3({\bf r}-{\bf r'})V(r) \)) without tensor forces. What the VPA lacks in generality, it makes up for in simplicity and the ability to better control the numerical accuracy. The VPA also gives a clean way to infer certain properties of the interaction from the phase shifts (e.g., the existence of a repulsive short-range component, or the presence of bound states) as illustrated below. For simplicity, here we only consider s-waves. Good references here are
Define the truncated potential \( V_\rho(r) \) by $$ V_\rho(r) = V(r) \theta(\rho-r) \;. $$ That is, it is the usual potential for \( r \leq \rho \), but identically zero beyond that. Then we define \( \delta(k,\rho) \) as the phase shift for \( V_\rho \) at momentum \( k \). The phase shift we want is \( \delta(k) = \lim_{\rho\rightarrow\infty} \delta(k,\rho) \). The basis of the variable phase method is a differential equation for \( \delta(k,r) \) at fixed \( k \) (again, this is the s-wave equation): $$ \frac{d\delta(k,r)}{dr} = -\frac{1}{k} 2M V(r) \sin^2[kr + \delta(k,r)], $$ which is a nonlinear first-order differential equation with initial condition \( \delta(k,0) = 0 \). Note that \( M \) here is actually the reduced mass of the NN system, and also ask yourself if there are factors of \( \hbar \) and/or \( c \) that have been set to 1. Think about how you would implement this in your favorite programming language. As an example, the Mathematica notebook SquareWellScattering.nb implements the VPA for a square well. Show that it reproduces the analytically known phase shifts for the square well result.
Changing to a different potential is trivial (see the illustration at the end of the notebook with a combined short-range repulsive square well and a mid-range attractive square well). Implement the toy NN potential that you are using in your momentum space matrix inversion code as a check on the former.
Show from the VPA differential equation that a fully attractive potential gives a positive phase shift and a fully negative potential gives a negative phase shift. This is the cleanest way to see why the s-wave phaseshifts (which change from positive to negative values at \( E_{lab}\approx 270 \) MeV in the \( ^1S_0 \) partial wave) imply a strong short-range repulsion for local NN potentials.
The VPA automatically builds in Levinson's theorem (\( \delta(0) = n\pi \)) about the number of bound states \( n \) and the phase shift at zero energy. How? Hint. What is the condition imposed on the phase shift at large energy for Levinson's theorem? Consider integrating \( d\delta(k,r)/dr \) in \( r \) from zero to infinity. Use \( \sin^2 x \leq 1 \) to put a bound on \( \delta(k) \).
Things to try numerically with the supplied Mathematica or Python notebooks:
Here follows a brief recipe and recommendation on how to write a report for each project.
The preferred format for the report is a PDF file. You can also use DOC or postscript formats or as an ipython notebook file. As programming language we prefer that you choose between C/C++, Fortran2008 or Python. The following prescription should be followed when preparing the report: