Calculate phase shift of square well

In the book, Introduction to Nuclear Reaction, by C. A. Bertulani, Chapter 2, there is a section that shows an example of phase shift from a square well. I follow the receipt and able to reproduce the result.

---

The actual method I did is little different. I use Mathematica to solve the radial Schrodinger equation in , the differential equation I solve is

where is the Woods-Saxon function, set the diffusiveness parameter to 0.001, the Woods-Saxon shape becomes square well.

In solving the differential equation, I set the right-hand side to be zero when . That solved to problem when , it has to deal with the . In fact, since , the right-hand side is always zero when .

For the solving range, I set

This is because the “first” peak of the Riccati-Bessel function appear after . This will make sure the range has 5 oscillations in the solution.

Next, I table the solution for the last 2 oscillations, find the maximum, normalize the solution, and then fit with

and extract , and .

The phase shift is

.

The negative is due to the is negatively defined in Mathematica.

---

The elastics cross section is

From the calculation, we found that, at very small energy, only s-wave scattering cross section is non-zero, which is agree with theory.

[20220531] update

The Mathematical notebook can be downloaded here.

Sensitivity of Woods-Saxon parameters in Glauber model

Experimentally, people measure the reaction cross section and compare with Glauber model. Since the inputs of the Glauber model are the nuclear density function and the NN cross section. The NN cross section is fixed by experiment. For an assumed density function, one or two parameters have to be fitted to the experimental cross section.

This time, we explore the sensitivity of the Woods-Saxon parameters in Glauber model using the method.

From the above figure, the sensitivity is largest for and is around 10 MeV/u. And the sensitivity drops for higher energy.

Note that match the experimental data even better.

We plot the ratio

And we can see that, all energy range has the same % sensitivity for , but different sensitivity for .

Next, we fixed the energy at 10 MeV/u and study the sensitivity of

We can see that,

So, if the measured cross section is smaller than 40 mb, the uncertainty of is about 0.01 fm at 10 MeV/u for 12C+12C. The last plot at below is the sensitivity for .

By measuring the reaction cross section at two energies with uncertainty less than 50 mb, if the model explain the data, the Woods-Saxon parameters should be extracted with small uncertainty.

I did not consider the uncertainty from the density. The error propagation from the density to the Woods-Saxon parameters need to be studied.

Woods-Saxon density in Glauber model

Following the receipt of the Glauber model. I did the calculation on 12C + 12C using Woods-Saxon density, and compared with the experimental result.

I did it using Mathematica. In order to avoid high dimension integration, the are first tablelized and interpolated. This reduced the complexity and speeded things up.

The NN cross section is taken from the reference S. K. Charagi and S. K. Gupta, PRC 41 (1990) 1610. It gives a parametrized from for the NN cross section. But I did not cross check the value for < 10 MeV.

I also compared the calculation with and without Coulomb correction. Also, finite-range interaction was also calculated.

---

Here is the results. The WS radius is .

The in this case is very close to a Gaussian, and can be approximated as

In order to match the 12C+12C cross section, I have to scale the by 80 times.

The above is the 12C+12C result from Glauber model using zero-range. The calculation is close to the experimental result from the reference. And the low energy is the Coulomb correction. I guess my Coulomb correction is not correct? But for low energy, should be correct, I am more less confident for the high energy. Anyway, the calculation is not that far away. Below is a figure taken from the reference for the 12C+12C reaction

Taken From S. K. Charagi and S. K. Gupta, PRC 41 (1990) 1610

I take the Coulomb KE to be half of the incident KE, i.e. using Center of mass KE. We get a better result at low energy. It seems that the center of mass KE should be used, but why? [update 2021-07-18] it is because the target will be pushed by the Coulomb force from the projectile. Since the net force on the two-particle system is zero, transform the coordinate into the CM frame, and the CM frame KE is used to calculate the closed distance.

Next I compare the calculation using a Gaussian finite-range with finite range parameter of 0.9 fm, below is the :

The finite-range interaction thickness function takes a similar shape but only different amplitude. This is understandable as the NN interaction not only when the nucleon-nucleon at the same position, but also can interaction at 0.9 fm apart. Since the transparency function is insensitive (I will discuss this at below) for the inner part of the thickness function. The increase of the amplitude would not cause much different, but the increase width will make a different.

---

At the end, I like to add some comments on the calculation.

Frist, in the calculation, I used a normalized Woods-Saxon density. Is it necessary? Since the probability of a NN-collision at impact parameter is , and for low energy is 100 ~ 1000 mb or 10 ~ 100 fm². This term is often not fulfilled.

Also, the transparency function, if we take the original form

This requires , but as I discussed above it is usually not the case. However, if we use the exp form

then there seems to be no limitation for the thickness function . So, if the normalization doesn’t matter, then what is the meaning of the thickness function times NN cross section?? not probability anymore? [update July-16, 2021] refer to the end of this post. The density function is normalized to mass number.

Second, I tried to follow the step in the reference S.K. Charagi and S.K. Gupta, but I cannot reproduce the Gaussian parameters as they did. In the paper, they match the Gaussian density to the surface of the Woods-Saxon function. The problem is the normalization of the Woods-Saxon density. I can use the table to interpolate the density and radius, and the Gaussian density is not normalized. Even I use their table values for the density and radius, I still unable to reproduce their result, except the 12C+12C, I can reproduce the trend, but the curve is still not the same as their.

Third, in the Coulomb correction, the wave number is non-relativistic.

Fourth, at energy below the Coulomb barrier, it becomes mostly elastic scattering, and the effective impact parameter, is very large, the reaction cross section should becomes zero. The Coulomb corrected Glauber model should reflect this.

Fifth, the interplay between the NN cross section, finite-range, the thickness function, and the transparency function. The transparency function is

We use a Gaussian approximation for the for illustration. And the transparency function becomes

This is the exactly the form from using Gaussian density.

In below, ,

We can see that, when , and . We can see that, the increase of the NN cross section, the increase of the WS density, and the finite-range effect, has similar effect and all contribute to increase the width of . When the width increased, the position for also increased, and the integration or the reaction cross section increased.

Here is the calculation for the 40Ca + 40Ca, the WS density also 80 times more.

Sixth, the looks like a Woods-Saxon function, can we find a simple relation between the density parameters to the ?

---

[update July-16, 2021]

In the above, we have to multiple the by 80 to reproduce the experimental result. At the calculation, the is calculated using Woods-Saxon function that is normalized by 1. Thus, this 80 is the multiplication of the mass number. But the multiplication of the mass number is . And 80 is only 55%. This 55% could be a correction to the density in the compensation for the radius and diffusiveness parameters.. For Woods-Saxon potential that normalized to the mass number,

We can see that the density has 20% difference, and we know that the density should be almost constant as the short-range nuclear force saturated at the center of nucleus. And the experimental nuclear density is 0.16 fm^-3.

I redo the calculation using , and the scaling is 144.

The result is closer to the experimental result.

Gaussian density for Glauber model

We are going to calculate the reaction cross section using Glauber model with Gaussian density. This already done in reference S. K. Charagi and S. K. Gupta, PRC 41 (1990) 1610, But I will provide more detail of the calculation. And also, in the reference, there is a factor 10, that is the conversion

---

The density function is

The finite-range function takes the form

when , the function becomes 2-D delta function .

The goal is to compute

The is a 6-dimension integral, but thanks to the Gaussian density and Gaussian finite-range function, all 6-dimension can be integrated individually in Cartesian coordinate.

---

Let integrate the beam direction.

Thus,

In Cartesian coordinate,

Lets integrate

Similarly for the , and we have

We can see a similar pattern in here, rearrange and simplify

Thus,

When , it reduces to zero-range.

We simplify :

where

We substitute at the last step. We can check that

which is expected as is the probability of a NN-collision at impact parameter .

---

The reaction cross section take the form

This is exponential of exponential, and integrate. I never deal with this before. Lets get a feeling on the . When is large, the inside exp approaches to zero, so, the outside gives 1. When is small, the inside exp gives a finite number, and the outside . In fact, if we plot using

We can see that, the is the probability of a NN-collision, and the is the transparency function.

Mathematica gives the integral as:

where is the Euler Gamma constant and is the exponential integral function.

The expression agree with the reference. We plot the below,

For a fixed , i.e. the size of the target and projectile, the reaction cross section change with non-linearly. But for a fixed the reaction cross section saturated to .

In fact, mathematically,

---

For some I don’t know why reason, the radial distance with

is called the strong absorption radius. In here, the density normalization is to the mass number, not 1.

requires or

When the energy is more than 100 MeV, the , this impiles , which is satisfied for most nuclei.

Then, a classical geometrical cross section is

Compare with

Radial density from Woods-Saxon potential

In this post, we estimated the radial density from a harmonic oscillator. Now, we do a similar job from Woods-Saxon potential.

The radial function (for neutron) is a solution of the Schrodinger equation :

I will use the python code to calculate the eigen wave functions.

The density will be

where is the orbital id, and is the occupation number of the orbital. For magic nuclei,

---

At the edge of the density function, it is similar to the Woods-Saxon, but at the center, depends on how much the s-orbital contributes, it can be a dip, an anti-dip, or flat.

Interestingly, the density saturates around 1 unit. There is no constraint, requirement, or normalization for that in the calculation. Is it a mathematical coincidence that it is a property of Woods-Saxon potential? or it is a general property for the Schrodinger equation?

---

[update 2022-03-05]

In text book, the mean of nuclear density is the same for all nuclei because the balance between the strong nuclear force and electromagnetic repulsion for proton, or the nuclear force is short-range. A situation similar to a liquid drop. And yes, this is the fundamental reason.

And the fact that the WS calculation reproduces a constant density for all nuclei is that the radius is

.

That makes sure the number density is

---

[update 2023-12-29]

The concepts of the wavefunction and the nuclear density are closely related. The square of the wave function in a given volume is the probability of finding a nucleon, while the density is the number of nucleons in a given volume. Mathematically, we demand the nuclear density to be

i.e. the total number of nucleons in the whole nucleus is recovered by integrating the density. From the wave function, we can construct

integrate out the angular dependence,

if we treat to be the radial density distribution, thus,

which is larger than it should be. Thus, a proper construction would be

Since the experimental nuclear mass density is nucleon/fm³. Thus, the neutron (or proton) mass density is approximately half and equal to 0.08 nucleon/fm³. From the above plots, if we divided them with , the density is roughly 0.08 nucleon/fm³. For light nuclei, since the Coulomb force is not strong and the proton and neutron wave functions are similar, we can simply multiply the neutron mass density by to approximate the mass density.

RMS radius of various nuclear density functions

In the Glauber model, one of the input is the density functions. People will adjust the parameters of the density function to fit the experimental cross section.

There are many choice of nuclear density functions, for example, Woods-Saxon, Harmonic Oscillator, Square well, and Gaussian etc. Among all of them, the root-mean-squared (RMS) radius of these functions usually represent the nuclear RMS radius. So, there are many RMS radius depends on the density function. In here, I will make some comparisons.

All density functions are normalized

and the mean-squared radius is

Using the above definition, I have

Name	Density function	RMS radius
Square Well
Woods-Saxon
Gaussian

The rms radius of Square well and Gaussian is simple ratio to the radius parameter.

As we shown before, the Woods-Saxon will approach to Square well, as it should be. But the Gaussian will be different. How to translate the rms radius between each of them? For example, If I say the rms radius is 4 fm +- 0.1 fm with Gaussian density function, what is the rms radius for Woods-Saxon ?

Here is other comparisons

---

I guess, to translate the rms radius between different density functions in Glauber model, we have to compare the . If the is similar. To define the similarity of functions, we using

When , the integral is equal to 1.

I compare the Woods-Saxon and Gaussian,

The 3rd order approximation for the peak (red) line is

Here is the plot for the Woods-Saxon with latex R=4 , a = 0.67$, and the Gaussian with ,

Glauber Model – 2

I compared the calculations from hard sphere and Woods-Saxon type density functions, to study the sensitivity of the Glauber model. The Glauber model calculate, for spherical nuclei

which is a 4-dimension integration. Once the is calculated, all other things follow. The reaction cross section is

When the is small and is large, make analog,

an approximation:

and call the exponential be transparency function

Because the form of the reaction cross section is also like

And the meaning of is also the transparency. Also, in semi-classical, the average angular momentum is , Thus, people would link

.

---

The hard sphere is relatively less time consuming, and the density function is easy to play with. It is easy to know that

,

the corresponding is

But it is hard to get the analytical form of the for Woods-Saxon, so I have to integrate the density numerically using 4-dimensional 25-points Gauss quadrature. Because it is 4-dimensional, the number of sum is , comparing to 2-dimensional 25-point GQ, which only need .

The normalized Woods-Saxon function is

Here is the comparison for radius 2 and radius 1, and diffusiveness is 0.3.

The simple hard-sphere means that I use 2-dimension GQ, “hard-sphere” means the 4-D GQ was used. And that verified the calculations for 4-dimension GQ works well.

The Woods-Saxon gives a smooth curve because the function itself is smooth. In Hard-sphere, the function is not smooth and that could create problem when using GQ method.

If we interpolate the two functions and integrate, i.e.

The hard-sphere gives 0.160056 and WS gives 0.158891. So, I think I can conclude that the Glauber model is not sensitive in these two cases. i.e., If someone try to extract the nucleus radius by fitting the experimental reaction cross section using Glauber model, using hard-sphere and WS is basically the same.

---

[Update July-16, 2021]

The “approximation” of mentioned above is not correct.

For ,

And we can absorb the into and this lead to re-define of the density function.

That means, the density function is no longer normalized to 1, but it normalized to equal the number of nucleons.

If we assume the density function is normalized to 1 and takes a Woods-Saxon shape with parameters , the maximum . And this limited the NN cross section be . And that translate to . This estimation is only based on one nucleus. , We calculated the energy limit versus mass number using

We can see that for most of the nuclei, the lab energy can be as low as 4 MeV/u (green line). And that is close to the Coulomb barrier.

Solving Woods-Saxon potential using Python

For seek of nice drawing, I code a python code for that. Although the scipy package can solve differential equation, the way is breaking a 2nd order equation into two 1st order equations. So, I simply move the RK4 method into python. Also, the plotly package can draw a nice plot.

import numpy
import math
import plotly.graph_objects as go
import scipy.integrate
import scipy.special
import scipy.interpolate

class WoodsSaxon:

  #grid setting
  rStart = 0.0001
  dr = 0.1
  nStep = 300

  #constant
  mn = 939.56539 #MeV/c2
  hbarc = 197.326979 #MeV.fm

  #RK4 constants
  parC = [0, 1./2, 1./2, 1.]
  parD = [1./6, 2./6, 2./6, 1./6]

  mu = 0
  R0 = 0
  RSO = 0

  E = 0 #energy of the state
  L = 0
  J = 0

  #inital condition
  solu0 = 0.0
  dsolu0 = 1.0

  trancation = 100 #wavefunction will trancate if value < max/ 100

  rpos = numpy.arange(rStart, rStart+nStep*dr, dr)
  kpos = numpy.arange(0, 1000, 1) #momentum in MeV/c
  kzpos = numpy.arange(0, 500, 1 ) #para-momenutem in MeV/c

  SolU = []
  maxSolU = 0.0
  WF = numpy.empty([nStep], dtype=float) # radial wave function
  maxWF = 0.0
  WFk = numpy.empty([len(kpos)], dtype=float) #momentum wave function
  WFkz = numpy.empty([len(kzpos)], dtype=float) #para-momentum  dist function

  minEigenE = -100 # the lowest eigen energy
  eigenEnergy = numpy.empty([0,6], dtype=float) #energy, N, L, J, dEnergy, rms(raidus)
  eigenWFList = numpy.empty([0,nStep], dtype=float) # list of eigen WF. 
  eigenWFkList = numpy.empty([0, len(kpos)], dtype=float) #list of eigen momentum wave function
  eigenWFkzList = numpy.empty([0, len(kzpos)], dtype=float) # list of para-momentum dist function

#====================== Setting

  # must declare at the initialization of the class
  def __init__(self, A, V0, r0, a0, VSO, rSO, aSO):
    self.A = A
    self.V0 = V0
    self.r0 = r0
    self.a0 = a0
    self.VSO = VSO
    self.rSO = rSO
    self.aSO = aSO
    self.mu = self.mn * A/(A+1)
    self.R0 = r0 * math.pow(A,1./3)
    self.RSO = rSO * math.pow(A, 1./3)
    #self.PrintParameters()
    self.minEigenE = V0 + 5  #assume the 1st s-state is 5 MeV above the V0

  # set energy, can be not eigen status
  def SetE(self, E):
    self.E = E

  # set the mometum L and J
  def SetLJ(self,L,J):
    self.L = L
    self.J = J

  # set the range in fm
  def SetRange(self, rStart, dr, nStep):
    self.rStart = rStart
    self.dr  = dr
    self.nStep = nStep
    self.WF=numpy.empty([nStep], dtype=float)
    rpos = numpy.arange(self.rStart, self.rStart+self.nStep*dr, self.dr)

  # Woods-Saxon potential
  def __WS(self, x):
    return self.V0/(1 + math.exp((x-self.R0)/self.a0))

  # spin-obital potential
  def __SO(self, x):
    return self.LS(self.L, self.J)*(self.VSO * math.exp((x-self.RSO)/self.aSO))/(self.aSO*math.pow(1+math.exp((x-self.RSO)/self.aSO),2))/x

  # LS coupling factor                                  
  def LS(self, L, J) :
    return (J*(J+1)-L*(L+1)-0.5*(1.5))/2.

  # The G-function, u''[r] = G[r, u[r], u'[r]]
  def __G(self, x, y, dy):
    #return  -2*x*dy -2*y  # solution of gaussian
    return 2*self.mu/math.pow(self.hbarc,2)*(self.__WS(x) - self.__SO(x) - self.E)*y + self.L*(1+self.L)/x/x*y

#============================== Radial wave function
  # Plot potential and wavefunction for give energy E
  def PlotAll(self, E, maxR = 20):
    self.E = E
    ypos = []
    rSpinObi = []
    SpinObi = []
    rL = []
    yL = []
    rtot = []
    ytot = []
    yE = []

    maxWS = abs(self.V0)

    for i in range(self.nStep):
      x = self.rStart + self.dr*i
      if( x < maxR ):
        central = self.__WS(x)
        ypos.append(central)
        yE.append(self.E)
        so = self.__SO(x)
        centrapetal = math.pow(self.hbarc,2)/2/self.mu*(self.L*(1+self.L))/x/x
        tot = central-so + centrapetal
        if abs(so) < maxWS  :
          rSpinObi.append(x)
          SpinObi.append(so)
        if abs(centrapetal) < maxWS :
          rL.append(x)
          yL.append(centrapetal)
        if abs(tot) < maxWS :
          rtot.append(x)
          ytot.append(tot)
    fig = go.Figure()
    fig.add_trace(go.Scatter(x=self.rpos, y=ypos, name="Woods-Saxon"))
    fig.add_trace(go.Scatter(x=rSpinObi, y=SpinObi, name="Spin-Orbit"))
    fig.add_trace(go.Scatter(x=rL, y=yL, name="Centrapetal"))
    fig.add_trace(go.Scatter(x=rtot, y=ytot, name="Total"))

    if E != 0.0 :
      fig.add_trace(go.Scatter(x=self.rpos, y=yE, name="Energy", marker=dict(color = 'black')))
      self.SolveByRK4()
      self.CovertSolUtoWF()
      fig.add_trace(go.Scatter(x=self.rpos, y=(self.WF/self.maxWF*maxWS + self.E), name="WaveFunction"))

    fig.update_xaxes(range=[0, maxR])
    fig.update_layout(xaxis_title="radius [fm]")
    fig.update_layout(yaxis_title="Energy [MeV]")
    fig.show()

  # plot the Eigen wave function of ID, will calculate Wave function
  def PlotEigenWF(self, ID, maxR = 20):
    if -1 < ID and ID < len(ws.eigenEnergy) : 
      self.SetLJ(self.eigenEnergy[ID, 2],self.eigenEnergy[ID, 3])
      self.PlotAll(self.eigenEnergy[ID,0], maxR)

  # plot all eigen wave function
  def PlotAllEigenWFs(self, isSquare = False, maxR = 20):
    fig = go.Figure()
    for i in range(len(self.eigenEnergy)) :
      haha = self.eigenEnergy[i]
      wf = self.eigenWFList[i]
      if isSquare == True :
        fig.add_trace(go.Scatter(x=self.rpos, y=(wf*wf), name=self.FormatOrbital(haha)))
      else:
        fig.add_trace(go.Scatter(x=self.rpos, y=wf, name=self.FormatOrbital(haha)))
    fig.update_xaxes(range=[0, maxR])
    fig.update_layout(xaxis_title="radius [fm]")
    fig.update_layout(yaxis_title=r"$\phi(r)$")
    fig.show()

  # Using Rungu-Kutta 4th method to solve u''[r] = G[r, u[r], u'[r]]
  def SolveByRK4(self):
    #initial condition
    self.SolU = [self.solu0]
    dSolU = [self.dsolu0]

    dyy = numpy.array([1., 0., 0., 0., 0.])
    dzz = numpy.array([1., 0., 0., 0., 0.])

    self.maxSolU = 0.0

    for i in range(self.nStep-1):
      r = self.rStart + self.dr * i
      y = self.SolU[i]
      z = dSolU[i]

      for j in range(4):
        j = j+1
        dyy[j] = self.dr * (z + self.parC[j-1]*dzz[j-1])
        dzz[j] = self.dr * self.__G(r + self.parC[j-1]*self.dr, y + self.parC[j-1]*dyy[j-1], z + self.parC[j-1]*dzz[j-1])
        #if( i < 10 ):
        #  print("%d, %f, (%f, %f ), %f, %f, %f" % (j, parC[j-1], dyy[j-1], dzz[j-1], dyy[j], dzz[j], G(r, y,z)))
      
      #if( i < 10 ): 
      #  print(dyy, dzz)

      dy = 0
      dz = 0
      for j in range(4):
        j = j+1
        dy = dy + self.parD[j-1]*dyy[j]
        dz = dz + self.parD[j-1]*dzz[j]

      #if( i < 10 ):
      #  print("%7.5f, %8.4f, %8.4f, %8.4f" % (r, y, G(r,y,z),  y +dy))

      self.SolU.append(y + dy)
      dSolU.append(z + dz)
      
      if abs(y+dy) > self.maxSolU and r < self.R0*6:
        self.maxSolU = abs(y+dy)

    return self.SolU

  # normalize SolU as int( SolU^2 r^2 ) = 
  def NormalizeSolU(self):
    SolU2 = []
    for i in range(self.nStep):
      if( self.rpos[i] > 5* self.R0 and self.SolU[i] < self.maxSolU / self.trancation):
        SolU2.append(0.0)
      else:
        SolU2.append(math.pow(self.SolU[i],2))
    norm = math.sqrt(scipy.integrate.simps(SolU2, self.rpos))
    for i in range(self.nStep):
      self.SolU[i] = self.SolU[i]/norm

  # convert the SolU = u[r] to radial wave function as R[r] = u[r]/r
  def CovertSolUtoWF(self):
    self.NormalizeSolU()
    #convert to wave function
    self.maxWF = 0.0
    for i in range(self.nStep):
      r = self.rStart + self.dr * i
      self.WF[i] = self.SolU[i]/r
      if( self.maxWF < abs(self.SolU[i]/r)):
        self.maxWF = abs(self.SolU[i]/r)
    if self.L == 0 :
      self.WF[0] = self.WF[1]

  # Find Eigen energies, should save the wave function one by one?
  def FindEigenEnergis(self, isPrint=True):
    Estep = 0.2 # this limit the weakly bound state
    self.eigenEnergy = numpy.empty([0,6], dtype=float)
    if self.A < 40 : 
      lmax = 3
    else :
      lmax = 6
    for l in range(lmax+1):
      if l == 0 :
        jList = [1./2]
      else :
        jList = [l + 1./2, l - 1./2]
      for j in jList:
        self.SetLJ(l, j)
        count = 0
        eigenN = 0
        EList = numpy.arange(self.minEigenE, 0, Estep) # assume there is no state with only 5 MeV above the V0
        oldhaha = []
        for E in EList:
          self.E = E
          Enow = E
          haha = self.SolveByRK4()
          if count > 0 :
            # if the tail change sign, use bisection method to find the best E    
            repeat = 0
            while oldhaha[-1] * haha[-1] < 0:
              EMiddle = (Enow + Epast)/2.
              if  (Enow - Epast < 1e-4 and abs(oldhaha[-1] - haha[-1]) < self.maxSolU * 0.01 ) or repeat > 200:
                if eigenN == 0 :
                  self.minEigenE = EMiddle #this set the starting energy for next j-value 

                self.eigenEnergy = numpy.append(self.eigenEnergy, [[EMiddle, eigenN, l, j, abs(Enow-Epast), 0]], axis=0)
                eigenN = eigenN + 1
                break
              else:
                self.E = EMiddle
                hahaMiddle = self.SolveByRK4()
                if oldhaha[-1] * hahaMiddle[-1] < 0  :
                  haha = hahaMiddle
                  Enow = EMiddle
                else:
                  oldhaha = hahaMiddle
                  Epast = EMiddle
              repeat = repeat + 1
          oldhaha = haha
          Epast = E
          count = count + 1

    #sorting energy energy with N, L, J
    self.eigenEnergy = self.eigenEnergy[self.eigenEnergy[:,0].argsort()]

    self.eigenWFList = numpy.empty([0,self.nStep], dtype=float) # list of eigen WF. 
    for i in range(len(self.eigenEnergy)):
      haha = self.eigenEnergy[i]
      self.E = haha[0]
      self.SetLJ(haha[2], haha[3])
      self.SolveByRK4()
      self.NormalizeSolU()
      self.CovertSolUtoWF()
      self.eigenWFList = numpy.append(self.eigenWFList, [self.WF], axis=0)
      self.eigenEnergy[i][5] = self.CalRMSRadius()

    if isPrint :
      self.PrintEigenEnergies()

  # Get Eigen radial wave function
  def GetEigenWF(self, ID):
    return self.eigenWFList[ID]

  # Calculate RMS radius for the WF stored in the momeory
  def CalRMSRadius(self):
    return math.sqrt(scipy.integrate.simps(self.WF*self.WF*numpy.power(self.rpos,4), self.rpos))

  # Print Woods-Saxon parameters
  def PrintParameters(self):
    print("======================================================")
    print("A = %d,  mu = %f MeV/c2" % (self.A, self.mu))
    print("V0  = %7.2f MeV, r0  = %6.3f(%6.3f) fm, a0  = %6.3f fm" % (self.V0, self.r0, self.R0, self.a0))
    print("VSO = %7.2f MeV, rSO = %6.3f(%6.3f) fm, aSO = %6.3f fm" % (self.VSO, self.rSO, self.RSO, self.aSO))
    print("Range = Start : %.6f fm,  End : %.6f | stepSize : %.2f fm, num. Step : %d" % (self.rStart, self.rStart + self.dr*self.nStep, self.dr, self.nStep))
    print("======================================================")

  # Format te orbital symbol from N, L, J using the eigen energy array
  def FormatOrbital(self, eigen):
    if eigen[2] == 0.0 :  Lsym = "s"
    if eigen[2] == 1.0 :  Lsym = "p"
    if eigen[2] == 2.0 :  Lsym = "d"
    if eigen[2] == 3.0 :  Lsym = "f"
    if eigen[2] == 4.0 :  Lsym = "g"
    if eigen[2] == 5.0 :  Lsym = "h"
    if eigen[2] == 6.0 :  Lsym = "i"
    if eigen[2] == 7.0 :  Lsym = "j"
    output = "%.0f%s%.0f/2" % (eigen[1], Lsym, 2*eigen[3])
    return output

  # print eigen energies
  def PrintEigenEnergies(self):
    print("======================================================")
    count = 0
    for haha in self.eigenEnergy:
      print("%2d | %4s   %14.10f MeV +- %14.10e | rms(Radius): %10.6f fm" % (count, self.FormatOrbital(haha), haha[0], haha[4], haha[5]))
      count = count + 1
    print("======================================================")

#================================= Momentum
  # Calculate momentum wave function. This use the WF in the memory.
  def __CalMomentumWF(self, L):
    for i in range(len(self.kpos)):
      kfm = self.kpos[i] / self.hbarc
      function = numpy.empty([self.nStep], dtype=float) # this is SpheicalBeseel(kr)*WF(r)*r^2
      for j in range(self.nStep):
        if( self.rpos[j] > 5 * self.R0 and self.WF[j] < self.maxWF / self.trancation) :
          function[j] = 0.
        else:
          sb = scipy.special.spherical_jn(int(L) , kfm * self.rpos[j])
          function[j] = sb * self.WF[j] * self.rpos[j] * self.rpos[j]
      self.WFk[i] = scipy.integrate.simps(function, self.rpos)
    maxWFk = numpy.sign(self.WFk[1])*max(abs(numpy.amax(self.WFk)), abs(numpy.amin(self.WFk)))
    self.WFk = self.WFk/maxWFk

  # Calculate eigen momentum wave function and store it
  def CalEigenMomentumWF(self):
    self.eigenWFkList = numpy.empty([0, len(self.kpos)], dtype=float)
    for i in range(len(self.eigenWFList)) :
      self.WF = self.eigenWFList[i]
      self.__CalMomentumWF(self.eigenEnergy[i][2])
      self.eigenWFkList = numpy.append(self.eigenWFkList, [self.WFk], axis=0)

  # [Obsolete?] normaalize momentum Wave Function as int(wfk^2 k^2) = 1
  def NormalizedMomentumWF(self):
    if len(self.WFk) == len(self.kpos) :
      norm = math.sqrt(scipy.integrate.simps((self.WFk*self.WFk*self.kpos*self.kpos), self.kpos))
      self.WFk/norm

  # return the eigen momentum wave function
  def GetEigenMomentumWF(self, ID):
    return self.eigenWFkList[ID]

  # Plot all Eigne momentum wave functions, this will recalculate from the eigen energy
  def PlotAllEigenMomentumWF(self, isSquare = False, maxk=500):
    fig = go.Figure()
    if len(self.eigenEnergy) != len(self.eigenWFkList) :
      self.CalEigenMomentumWF()

    for i in range(len(self.eigenWFkList)) :
      haha = self.eigenEnergy[i]
      if isSquare == False :
        wfk = self.eigenWFkList[i]
      else:
        wfk = numpy.power(self.eigenWFkList[i],2)
      fig.add_trace(go.Scatter(x=self.kpos, y=wfk, name=self.FormatOrbital(haha)))
    fig.update_xaxes(range=[0, maxk])
    fig.update_layout(xaxis_title="Momemtum [MeV/c]")
    if isSquare == True :
      fig.update_layout(yaxis_title=r"$\phi^2(k)$")
    else:
      fig.update_layout(yaxis_title=r"$\phi(k)$")
      
    fig.show()

  # calculate the parallel momentum distribution using the WFk in the memory
  # the para-momentum distribution is f(kz) = int( WFk(sqrt(kz^2 + kr^2 ))^2 kr )
  def __CalParaMomentumDist(self):
    wfk = scipy.interpolate.interp1d(self.kpos, self.WFk)
    krList = self.kzpos
    for i in range(len(self.kzpos)):
      kz = self.kzpos[i]
      kList = numpy.sqrt(kz*kz + krList*krList)
      func = wfk(kList)
      self.WFkz[i] = scipy.integrate.simps(func * func * krList, krList )
    maxWFkz = max(abs(numpy.amax(self.WFkz)), abs(numpy.amin(self.WFkz)))
    self.WFkz= self.WFkz/maxWFkz

  # Calculate eigen para-momentum distribution
  def CalEigenParaMomentumDist(self):
    self.eigenWFkzList = numpy.empty([0, len(self.kzpos)], dtype=float)
    if len(self.eigenEnergy) != len(self.eigenWFkList) :
      self.CalEigenMomentumWF()

    for haha in self.eigenWFkList :
      self.WFk = haha
      self.__CalParaMomentumDist()
      self.eigenWFkzList = numpy.append(self.eigenWFkzList, [self.WFkz], axis=0)

  # Plot all parallel momentum distribution 
  def PlotAllParaMomentumDists(self, maxk = 500):
    if len(self.eigenEnergy) != len(self.eigenWFkzList) :
      self.CalEigenParaMomentumDist()

    fig = go.Figure()
    for i in range(len(self.eigenWFkzList)) :
      haha = self.eigenEnergy[i]
      wfkz = self.eigenWFkzList[i]
      fig.add_trace(go.Scatter(x=self.kpos, y=wfkz, name=self.FormatOrbital(haha)))
    fig.update_xaxes(range=[0, maxk])
    fig.update_layout(xaxis_title="para-Momemtum [MeV/c]")
    fig.show()

The usage is very straight forward.

[A, V0, r0, a0, Vso, rso, aso ] = [15, -48.5, 1.24, 0.67, 16.5, 1.24, 0.67]
ws = WoodsSaxon(A, V0, r0, a0, Vso, rso, aso)
ws.SetLJ(1, 1.5)
ws.PlotAll(-13.191423125, 30)

The result is

It can also find all eigen energies

ws = WoodsSaxon(12, -58, 1.22, 0.69, 35, 0.75, 0.56)
ws.FindEigenEnergis()
ws.PlotAllEigenWFs(False)
ws.PlotAllEigenMomentumWF(False, 500)
ws.PlotAllParaMomentumDists() #parallel momentum distribution

The result is

======================================================
A = 12
V0  =  -58.00 MeV, r0  =  1.220( 2.793) fm, a0  =  0.690 fm
VSO =   35.00 MeV, rSO =  0.750( 1.717) fm, aSO =  0.560 fm
Range = Start : 0.000100 fm,  End : 30.000100 | stepSize : 0.10 fm, num. Step : 300
======================================================
found E : -31.1438878912, NLJ: 00(1/2) | E-limit: (-3.114389e+01, -3.114389e+01),  U-limit: (3.275552e-05, -1.037775e-02) | maxSolU: 0.975954
found E :  -1.6901489258, NLJ: 10(1/2) | E-limit: (-1.690161e+00, -1.690137e+00),  U-limit: (-3.547450e-03, 3.666512e-03) | maxSolU: 0.872269
found E : -16.4853724474, NLJ: 01(3/2) | E-limit: (-1.648537e+01, -1.648537e+01),  U-limit: (1.056252e-01, -1.962114e-01) | maxSolU: 40.553369
found E :  -9.5338185416, NLJ: 01(1/2) | E-limit: (-9.533819e+00, -9.533819e+00),  U-limit: (1.387732e-01, -1.342873e-01) | maxSolU: 53.446166
found E :  -1.3501820670, NLJ: 02(5/2) | E-limit: (-1.350188e+00, -1.350176e+00),  U-limit: (2.848967e-01, -1.905022e+01) | maxSolU: 2095.640491
======================================================
 0 | 0s1/2   -31.1438878912 MeV +- 3.5527136788e-15 | rms(Radius):   1.976819 fm
 1 | 0p3/2   -16.4853724474 MeV +- 1.1642242725e-11 | rms(Radius):   2.579059 fm
 2 | 0p1/2    -9.5338185416 MeV +- 1.4901164747e-09 | rms(Radius):   2.976584 fm
 3 | 1s1/2    -1.6901489258 MeV +- 2.4414062500e-05 | rms(Radius):   5.055909 fm
 4 | 0d5/2    -1.3501820670 MeV +- 1.2207031250e-05 | rms(Radius):   3.667586 fm
======================================================

Momentum distribution From Woods-Saxon potential

In this post, we derived the transformation between radial distribution and momentum distribution. This is a simple “Fourier Transform” in 3D with spherical Bessel function. In that post we also plot some simple functions.

The effective potential for most nuclei can be approximated (or well described) using the Woods-Saxon potential. Here we study the 15C, find the radial wave functions and the corresponding momentum distributions.

There are many way to calculate the radial function for a give set of Woods-Saxon parameters, using C++ code, MS excel, or Mathematica. In this post, I will use Mathematica. I found that MS Excel has the Bessel function, but it only accept integer n, but spherical Bessel needs half integer n, so it is not easy for calculating the momentum distribution.

In Mathematica, solving the Schrodinger equation for the radial part is extremely easy, for give energy , angular momentum , total spin

sol = NDSolve[{u''[
r] == (2 m)/hbar^2 (V0/(1 + Exp[(r - R0)/a0]) - 
1/r (Exp[(r - RSO)/aSO] VSO)/(a0 (1 + Exp[(r - RSO)/aSO)]^2)
LS - T) u[r] + (L (L + 1))/r^2 u[r], u[0] == 0, u'[0] == 1}, 
u, { r, 0, 30}];

where is the potential depth in MeV, is the hald-radius in fm, is the diffusiveness parameter in fm, and the subscript for central and spin-orbit. is the LS coupling factor. is the reduced radial function.

The normalization of the reduced radial function is

---

Usually, we don’t know the energy for a give state, we have to scan through all possible energy from . When the energy is not the eigen energy, the tail of the wave function will diverge and oscillate with respect to the energy. Thus, we can use this property to find the eigen energy that the tail of the wave function is flat.

In here, I am requiring the sum to be minimum.

For 15C, the single particle energies for the 1s1/2 and 0d5/2 is known from the 14C(d,p) reaction, which are -1.2181 MeV and -0.4741 MeV, respectively.

After searching for best fit parameters, I found

are best fit. The energies of different states are:

Orbital	Energy (WS) [MeV]	Energy (exp) [MeV]
0s1/2	-25.634
0p3/2	-12.522
0p1/2	-10.351
1s1/2	-1.210	-1.2181
0d5/2	-0.463	-0.474

And, I have the radial function.

We can see that, the 1s1/2 orbital has largest radial distribution and the 0s1/2 has the smallest.

Next, we transform the radial distribution to the momentum distribution. The unit conversion is that

In the above figure, we scaled the distribution so that they have similar magnitude. We can see, the momentum distribution is smallest for the 1s1/2 as expected.

It is interesting that the momentum spread of the 0s1/2 and 1s1/2 are very different, the spread of the 0s1/2 is almost 3 times more than that of the 1s1/2, and that of the p-wave also about double.

And the d-wave momentum distribution has 2 bumps, this looks something wrong. (??? I am no solid answer for that. I suspect this is related to the tail of the wave function, when the tail is not zero, the spherical Bessel J times r^2 will magnify the tail. )

---

The radial solution of the Woods-Saxon potential (or the Schrodinger equation?) for the bound state is always take this form:

where is the short range behaviour, is the long range behaviour, and is a modifier for the middle part of the function.

Lets assume the for simplicity.

The momentum distribution using this simplified radial function is

Normalized with

Here are the plot with for

the is similar to orbital.

The momentum distribution is good for distinguish and , but not so good for higher .

And it is counter intuitive that more spread of the radial function, the momentum distribution also spread out more. I guess the effect of the angular momentum has to take into account. For same , more spread of momentum still means the radial spread is smaller.

Coulomb Energy for a s-orbital under Woods-Saxon potential

Pick up from the Coulomb Energy calculation, For beta-decay, a neutron change to a proton and emit an electron. I assume the proton-orbital is as same as that of the neutron. In particular, I am working on 11Be beta-decay.

The outer most neutron is in the 1s1/2 orbital with separation energy of 0.502 MeV.

First, I have to find the 1s1/2 neutron wave function using Woods-Saxon potential. Since it is the s-orbital, the spin-orbital part can be ignored. And only 3 parameters left: the potential depth, the half-depth radius, and the diffuseness parameter. The diffuseness parameter is fixed at 0.67 fm. The half-depth radius is for . ( May be I should use ? ). I adjusted the depth to give the -0.502 MeV neutron separation energy. The depth is -51.56 MeV.

The proton/charge distribution is also same as the Woods-Saxon potential, but with a normalization that

In general, a normalized spherical Woods-Saxon potential is

,where is the polylogarithm function.

The charge distribution and the squared s-wave function is

Now, we have to integrate the wave function and the charge distribution

The potential is follow

The Coulomb potential of the Woods-Saxon shape is “smaller” than that of the Uniform sphere is expected.

The Coulomb energy with Z = 4 is then

The Coulomb energy using the potential of a uniform sphere with Z  = 4 is

---

The experimental Coulomb displacement energy for the 11Be (T = 3/2 ) state can be extracted using the mass and the excited state of 11B(T=3/2). The excitation energy of 11B(T=3/2) = 12.554 MeV. The experimental Coulomb energy is

The calculation is ~ 0.7 MeV smaller. The difference should be due to the shape of the wave function, and also co-related to the shape of the potential, as the Woods-Saxon shape does not accurately descript light nuclei. (or, take the electron mass into account ?? ) 

I changed the WS parameter fm. So that the wave function is more concentrated toward the center that the 1st node is around 1.9 fm (the potential depth becomes -59.2 MeV). The Coulomb energy with uniform sphere is 1.47 MeV. 

a [fm]	Coulomb Energy with Uniformly charge sphere [MeV] (R=1.25(11)^3 = 2.78 fm)
0.67	1.24
0.3	1.47
0.1	1.62
0.02	1.65

R [fm]	Coulomb Energy with Uniformly charge sphere [MeV] ( a= 0.67 fm)
2.78	1.24
2.5	1.28
2.1	1.34

After a few trials, I found that when R=2.2 fm, a = 0.1 fm, The Coulomb energy of the wavefunction of the Woods-Saxon potential and Uniformly charge sphere is 1.85 MeV. The 1st node of the 1s1/2 wavefunction is around 1.5 fm. But the depth of the potential is -98 MeV, which is too deep.