The aim here is to develop a program which solves the Lippman-Schwinger equation for a simple parametrization for the \( ^1S_0 \) partial wave. This partial wave is given by a central force only and is parametrized in coordinate space as $$ V(r)=V_a \frac{e^{-ax}}{x}+V_b \frac{e^{-bx}}{x}+V_c \frac{e^{-cx}}{x} $$ with \( x=\mu r \), \( \mu=0.7 \) fm (the inverse of the pion mass), \( V_a=-10.463 \) MeV and \( a=1 \), \( V_b=-1650.6 \) MeV and \( b=4 \) and \( V_c=6484.3 \) MeV and \( c=7 \).
a) Find an analytical expression for the Fourier-Bessel transform (Hankel transform) to momentum space for \( l=0 \) using $$ \left\langle k \right | V_{l} \left | k' \right\rangle = \int j_l(kr)V(r)j_l(k'r)r^2dr. $$
b) Write a small program which calculates the latter expression and use this potential to compute the \( T \)-matrix at positive energies for \( l=0 \). Compare your results to those obtained with a box potential given by $$ V(r)=\left\{ \begin{array}{cc} -V_0& r < R_0 \\ 0 & r > R_0 \end{array} \right. $$ Make a plot of the two \( T \)-matrices for energies up to 300 MeV in the lab frame and comment your results.
Finally, a warning, the above central potential is fitted to data from approximately 20 MeV to some 300 MeV. This means that results outside the data set should be taken seriously.
The following Fortran 95 program solves the above Lippmann-Schwinger equation. Python and C++ codes will be added later. Discussions of the results will also be added.
C *******************************************************
C Example program used to evaluate the
C T-matrix following Kowalski's method (eqs V88 & V89
C in Brown and Jackson)
C for positive energies only
C The program is set up for S-waves only
C Coded by : Morten Hjorth-Jensen
C Language : FORTRAN 90
C *******************************************************
C ******************************
C Def of global variables
C ******************************
MODULE constants
DOUBLE PRECISION , PUBLIC :: p_mass, hbarc
PARAMETER (p_mass =938.926D0, hbarc = 197.327D0)
END MODULE constants
MODULE mesh_variables
INTEGER, PUBLIC :: n_rel
PARAMETER(n_rel=48)
DOUBLE PRECISION, ALLOCATABLE, PUBLIC :: ra(:), wra(:)
END MODULE mesh_variables
C ******************************
C Begin of main program
C ******************************
PROGRAM t_matrix
USE mesh_variables
IMPLICIT NONE
INTEGER istat
ALLOCATE( ra (n_rel), wra (n_rel), STAT=istat )
CALL rel_mesh ! rel mesh & weights
CALL t_channel ! calculate the T-matrix
DEALLOCATE( ra,wra, STAT=istat )
END PROGRAM t_matrix
C *********************************************************
C obtain the t-mtx
C vkk is the box potential
C f_mtx is equation V88 og Brown & Jackson
C *********************************************************
SUBROUTINE t_channel
USE mesh_variables
IMPLICIT NONE
INTEGER istat, i,j
DIMENSION vkk(:,:),f_mtx(:),t_mtx(:)
DOUBLE PRECISION, ALLOCATABLE :: vkk,t_mtx,f_mtx
DOUBLE PRECISION t_shell
ALLOCATE(vkk (n_rel,n_rel), STAT=istat)
CALL v_pot_yukawa(vkk) ! set up the box potential in routine vpot
ALLOCATE(t_mtx (n_rel), STAT=istat) ! allocate space in heap for T
ALLOCATE(f_mtx (n_rel), STAT=istat) ! allocate space for f
DO i=1,n_rel ! loop over positive energies e=k^2
CALL f_mtx_eq(f_mtx,vkk,i) ! solve eq. V88
CALL principal_value(vkk,f_mtx,i,t_shell) ! solve Eq. V89
DO j=1,n_rel ! the t-matrix
t_mtx(j)=f_mtx(j)*t_shell
IF(j == i) WRITE(6,*) ra(i) ,t_mtx(i)
c & DATAN(-ra(i)*t_mtx(i))
ENDDO
ENDDO
DEALLOCATE(vkk , STAT=istat)
DEALLOCATE(t_mtx, f_mtx, STAT=istat)
1000 FORMAT( I3, 2F12.6)
END SUBROUTINE t_channel
C ***********************************************************
C The analytical expression for the box potential
C of exercise 1 and 12
C vkk is in units of fm^-2 (14 MeV/41.47Mevfm^2, where
C 41.47= \hbarc^2/mass_nucleon),
C ra are mesh points in rel coordinates, units of fm^-1
C ***********************************************************
SUBROUTINE v_pot_box(vkk)
USE mesh_variables
USE constants
IMPLICIT NONE
INTEGER i,j
DOUBLE PRECISION vkk, r_0, v_0, a, b, fac
PARAMETER(r_0=2.7d0,v_0=0.33759d0) !r_0 in fm, v_0 in fm^-2
DIMENSION vkk(n_rel,n_rel)
DO i=1,n_rel ! set up the free potential
DO j=1,i-1
a=ra(i)+ra(j)
b=ra(i)-ra(j)
fac=v_0/(2.d0*ra(i)*ra(j))
vkk(j,i)=fac*(DSIN(a*r_0)/a-DSIN(b*r_0)/b)
vkk(i,j)=vkk(j,i)
ENDDO
fac=v_0/(2.d0*(ra(i)**2))
vkk(i,i)=fac*(DSIN(2.d0*ra(i)*r_0)/(2.d0*ra(i))-r_0)
ENDDO
END SUBROUTINE v_pot_box
C ***********************************************************
C The analytical expression for a Yukawa potential
C in the l=0 channel
C vkk is in units of fm^-2,
C ra are mesh points in rel coordinates, units of fm^-1
C The parameters here are those of the Reid-Soft core
C potential, see Brown and Jackson eq. A(4)
C ***********************************************************
SUBROUTINE v_pot_yukawa(vkk)
USE mesh_variables
USE constants
IMPLICIT NONE
INTEGER i,j
DOUBLE PRECISION vkk, mu1, mu2, mu3, v_1, v_2, v_3, a, b, fac
PARAMETER(mu1=0.49d0,v_1=-0.252d0)
PARAMETER(mu2=7.84d0,v_2=-39.802d0)
PARAMETER(mu3=24.01d0,v_3=156.359d0)
DIMENSION vkk(n_rel,n_rel)
DO i=1,n_rel ! set up the free potential
DO j=1,i
a=(ra(j)+ra(i))**2
b=(ra(j)-ra(i))**2
fac=1./(4.d0*ra(i)*ra(j))
vkk(j,i)=v_1*fac*DLOG((a+mu1)/(b+mu1))+
& v_2*fac*DLOG((a+mu2)/(b+mu2))+
& v_3*fac*DLOG((a+mu3)/(b+mu3))
vkk(i,j)=vkk(j,i)
ENDDO
ENDDO
END SUBROUTINE v_pot_yukawa
C **************************************************
C Solves eq. V88
C and returns < p | f_mtx | n_pole =k>
C **************************************************
SUBROUTINE f_mtx_eq(f_mtx,vkk,n_pole)
USE mesh_variables
USE constants
IMPLICIT NONE
INTEGER i, j, int, istat, n_pole
DOUBLE PRECISION f_mtx,vkk,dp,deriv,pih,xsum
DIMENSION dp(1),deriv(1)
DIMENSION f_mtx(n_rel),vkk(n_rel,n_rel),a(:,:),fu(:),q(:),au(:)
DOUBLE PRECISION, ALLOCATABLE :: fu, q, au, a
pih=2.D0/ACOS(-1.D0)
ALLOCATE( a (n_rel,n_rel), STAT=istat)
DO i=1,n_rel
ALLOCATE(fu(n_rel), q(n_rel), au(n_rel), STAT=istat)
DO j=1,n_rel
fu(j)=vkk(i,j)-vkk(i,n_pole)*vkk(n_pole,j)/
& vkk(n_pole,n_pole)
ENDDO
DO j=1,n_rel
IF(j /= n_pole ) THEN ! regular part
a(j,i)=pih*fu(j)*wra(j)*(ra(j)**2)/
& (ra(j)**2-ra(n_pole)**2)
ELSEIF(j == n_pole) THEN ! use l'Hopitals rule to get pole term
dp(1)=ra(j)
CALL spls3(ra,fu,n_rel,dp,deriv(1),1,q,au,2,0)
a(j,i)=pih*wra(j)*ra(j)/2.d0*deriv(1)
ENDIF
ENDDO
DEALLOCATE(fu, q, au, STAT=istat) ! free space in heap
a(i,i)=a(i,i)+1.D0
ENDDO
CALL matinv(a, n_rel) ! Invert the matrix a
DO j=1,n_rel ! multiply inverted matrix a with dim less pot
xsum=0.D0
DO i=1,n_rel
xsum=xsum+a(i,j)*vkk(i,n_pole)/vkk(n_pole,n_pole) ! gives f-matrix in V88
ENDDO
f_mtx(j)=xsum
ENDDO
DEALLOCATE (a, STAT=istat)
END SUBROUTINE f_mtx_eq
C **************************************************
C Solves the principal value integral of V89
C returns the t-matrix for k=k, t_shell
C **************************************************
SUBROUTINE principal_value(vkk,f_mtx,n_pole,t_shell)
USE mesh_variables
IMPLICIT NONE
DOUBLE PRECISION vkk, f_mtx, t_shell, sum, pih, deriv, term
DIMENSION deriv(1)
DIMENSION vkk(n_rel, n_rel), f_mtx(n_rel),fu(:), q(:), au(:)
DOUBLE PRECISION, ALLOCATABLE :: fu, q, au
INTEGER n_pole, i, istat
ALLOCATE(fu(n_rel), q(n_rel), au(n_rel), STAT=istat)
sum=0.D0
pih=2.D0/ACOS(-1.D0)
DO i=1,n_rel
fu(i)=vkk(n_pole,i)*f_mtx(i)
ENDDO
DO i=1,n_pole-1 ! integrate up to the pole - 1 mesh
term=fu(i)*(ra(i)**2)-fu(n_pole)*(ra(n_pole)**2)
sum=sum+pih*wra(i)*term/(ra(i)**2-ra(n_pole)**2)
ENDDO ! here comes the pole part
CALL spls3(ra,fu,n_rel,ra(n_pole),deriv,1,au,q,2,0)
sum=sum+pih*wra(n_pole)*(fu(n_pole)+ra(n_pole)*deriv(1)/2.d0)
DO i=n_pole+1,n_rel ! integrate from pole + 1mesh pt to infty
term=fu(i)*(ra(i)**2)-fu(n_pole)*(ra(n_pole)**2)
sum=sum+pih*wra(i)*term/(ra(i)**2-ra(n_pole)**2)
ENDDO
t_shell=vkk(n_pole,n_pole)/(1.d0+sum)
DEALLOCATE (fu, q, au, STAT=istat)
END SUBROUTINE principal_value
C ***********************************************
C Set up of relative mesh and weights
C ***********************************************
SUBROUTINE rel_mesh
USE mesh_variables
IMPLICIT NONE
INTEGER i
DOUBLE PRECISION pih,u,s,xx,c,h_max
PARAMETER (c=0.75)
DIMENSION u(n_rel), s(n_rel)
pih=ACOS(-1.D0)/2.D0
CALL gausslegendret (0.D0,1.d0,n_rel,u,s)
DO i=1,n_rel
xx=pih*u(i)
ra(i)=DTAN(xx)*c
wra(i)=pih*c/DCOS(xx)**2*s(i)
ENDDO
END SUBROUTINE rel_mesh
C *********************************************************
C Routines to do mtx inversion, from Numerical
C Recepies, Teukolsky et al. Routines included
C below are MATINV, LUDCMP and LUBKSB. See chap 2
C of Numerical Recepies for further details
C Recoded in FORTRAN 90 by M. Hjorth-Jensen
C *********************************************************
SUBROUTINE matinv(a,n)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION a(n,n)
INTEGER istat
DOUBLE PRECISION, ALLOCATABLE :: y(:,:)
INTEGER, ALLOCATABLE :: indx(:)
ALLOCATE (y( n, n), STAT =istat)
ALLOCATE ( indx (n), STAT =istat)
DO i=1,n
DO j=1,n
y(i,j)=0.
ENDDO
y(i,i)=1.
ENDDO
CALL ludcmp(a,n,indx,d)
DO j=1,n
call lubksb(a,n,indx,y(1,j))
ENDDO
DO i=1,n
DO j=1,n
a(i,j)=y(i,j)
ENDDO
ENDDO
DEALLOCATE ( y, STAT=istat)
DEALLOCATE ( indx, STAT=istat)
END SUBROUTINE matinv
SUBROUTINE LUDCMP(A,N,INDX,D)
IMPLICIT REAL*8(A-H,O-Z)
PARAMETER (TINY=1.0E-20)
DIMENSION A(N,N),INDX(N)
INTEGER istat
DOUBLE PRECISION, ALLOCATABLE :: vv(:)
ALLOCATE ( vv(n), STAT = istat)
D=1.
DO I=1,N
AAMAX=0.
DO J=1,N
IF (ABS(A(I,J)) > AAMAX) AAMAX=ABS(A(I,J))
ENDDO
IF (AAMAX == 0.) PAUSE 'Singular matrix.'
VV(I)=1./AAMAX
ENDDO
DO J=1,N
IF (J > 1) THEN
DO I=1,J-1
SUM=A(I,J)
IF (I > 1)THEN
DO K=1,I-1
SUM=SUM-A(I,K)*A(K,J)
ENDDO
A(I,J)=SUM
ENDIF
ENDDO
ENDIF
AAMAX=0.
DO I=J,N
SUM=A(I,J)
IF (J > 1)THEN
DO K=1,J-1
SUM=SUM-A(I,K)*A(K,J)
ENDDO
A(I,J)=SUM
ENDIF
DUM=VV(I)*ABS(SUM)
IF (DUM >= AAMAX) THEN
IMAX=I
AAMAX=DUM
ENDIF
ENDDO
IF (J /= IMAX)THEN
DO K=1,N
DUM=A(IMAX,K)
A(IMAX,K)=A(J,K)
A(J,K)=DUM
ENDDO
D=-D
VV(IMAX)=VV(J)
ENDIF
INDX(J)=IMAX
IF(J /= N)THEN
IF(A(J,J) == 0.) A(J,J)=TINY
DUM=1./A(J,J)
DO I=J+1,N
A(I,J)=A(I,J)*DUM
ENDDO
ENDIF
ENDDO
IF(A(N,N) == 0.) A(N,N)=TINY
DEALLOCATE ( vv, STAT = istat)
END SUBROUTINE LUDCMP
SUBROUTINE LUBKSB(A,N,INDX,B)
implicit real*8(a-h,o-z)
DIMENSION A(N,N),INDX(N),B(N)
II=0
DO I=1,N
LL=INDX(I)
SUM=B(LL)
B(LL)=B(I)
IF (II /= 0)THEN
DO J=II,I-1
SUM=SUM-A(I,J)*B(J)
ENDDO
ELSE IF (SUM /= 0.) THEN
II=I
ENDIF
B(I)=SUM
ENDDO
DO I=N,1,-1
SUM=B(I)
IF (I < N)THEN
DO J=I+1,N
SUM=SUM-A(I,J)*B(J)
ENDDO
ENDIF
B(I)=SUM/A(I,I)
ENDDO
END SUBROUTINE lubksb
c) The parameters of the box potential are chosen to fit a potential with a bound state at zero energy. What does this mean for your \( T \)-matrix with this potential when \( k\rightarrow 0 \)?