One-Pion Exchange NN Scattering

Set up the OPE potential

In this notebook, we compute the phase shifts of the one-pion exchange (OPE) potential using the Variable Phase Approach (VPA) generalized to l >0. The reduced mass is μ = Mn/2.

gA = 1.26; mpi= 139.; Mn = 939.; fpi=92.4;  hbarc = 197.33;

mu = Mn/2;

Set Snn, Lnn, Jnn, Tnn [remember that -1^(Tnn+Snn+Lnn) = -1]

Snn = 1; Lnn = 6; Jnn = 6; Tnn = 0;    (* 3I6, seems to work *)

(* Snn = 0; Lnn = 6; Jnn = 6; Tnn = 1; *)    (* 1I6, seems to work *)

(* Snn = 1; Lnn = 6; Jnn = 5; Tnn = 0; *)      (* 3I5, seems to work *)

(* Snn = 1; Lnn = 6; Jnn = 7; Tnn = 0; *)    (* 3I7, seems to work *)

Set the spin and isospin factors

sigmaFactor := If [Snn==0, -3, 1];
tauFactor := If [Tnn==0, -3, 1];

We multiply S12Factor by Snn to make it zero for singlet states.

S12Factor := Snn*If[Lnn==Jnn,2/3,
                                 If[Lnn==Jnn-1,-2/3 (Jnn-1)/(2Jnn+1),
                                 If[Lnn == Jnn+1, -2/3(Jnn+2)/(2Jnn+1), 0]
                                    ]
                                ]

sigmaFactor

tauFactor

S12Factor

Define the OPE:

Yukawa[z_]:= Exp[-z]/z

YukawaTensor[z_] := Yukawa[z]*(1 + 3/z + 3/z^2)

Vope[r_] := 1/(12 Pi) gA^2/(4fpi^2) mpi^3 tauFactor *(3*YukawaTensor[mpi*r]*S12Factor+Yukawa[mpi*r]*sigmaFactor)

Plot[Vope[r/hbarc],{r,0.01,6}]  (* Uncomment to check potential *)

Other definitions

convert = 180/Pi;  (* to convert from radians to degrees for plotting *)

Ek[k_] := k^2/(2μ)   (* Convenient definition of kinetic energy *)

Define the “hat” functions we need using the current value of Lnn.

jhat[z_] := z*SphericalBesselJ[Lnn,z]

nhat[z_] := z*SphericalBesselY[Lnn,z]

Phase shift from the Variable Phase Approach

Use the formula: deltarho'(r)==-((2 μ Vope(r) (cos(deltarho(r))*jhat(k r)-sin(deltarho(r))nhat(k r)))^2)/k) 
with the initial condition that deltarho(0) = 0.

We’ll need a numerical differential equation solver.  In Mathematica this is NDSolve:

?NDSolve

In principle we integrate out to infinity.  In practice, we choose a value of k and integrate out to Rmax, chosen to be well beyond the range of the potential (at which point the right side of the equation for deltarho’(r) is zero), and then evaluate deltarho(Rmax), which is the phase shift for momentum k.  Here is an implementation:

deltaVPA[Elab_] :=(
krel = Sqrt[Mn Elab/2]; (* krel from Elab in MeV *)
Rmin= 0.3; (* start the integration away from zero *)
Rmax = 30; (* integrate out to Rmax in fm *)ans = NDSolve[{deltarho'[r]==-(1/krel)Mn Vope[r/hbarc]/hbarc (Cos[deltarho[r]]*jhat[krel  r/hbarc] - Sin[deltarho[r]]*nhat[krel r/hbarc])^2,deltarho[Rmin]==0},deltarho,{r,Rmin,Rmax},AccuracyGoal->6,PrecisionGoal->6];
(deltarho[r] /. ans)[[1]] /. r->Rmax  (* evaluate at r=Rmax *)
)

Do a quick check against the analytic result with sample values for k and V0 to see the accuracy we are getting (we use NumberForm to get 10 digits displayed):

NumberForm[ deltaVPA[10],10]

To get more digits correct, increase AccuracyGoal and PrecisionGoal (at the cost of slower evaluation of the function).

Check the cutoff phase shift out to Rmax.  Does this plot make sense?

Plot[Evaluate[convert*deltarho[r]/.ans],{r,Rmin,Rmax}]

Now that it seems to be working, try plotting δ against lab energy to test against NN-online.

Plot[convert*deltaVPA[Elab],{Elab,0.001,150},PlotRange->Full]

Make a table to compare directly to table output from NN-online.

Table[{Elab,convert*deltaVPA[Elab]},{Elab,1,150,1}]//TableForm