Correct Woods-Saxon distribution diffusivity (surface thickness) for finite nucleon size.

corrected-woods-saxon

#!/usr/bin/env python3

"""
Determine how finite nucleon size affects the Woods-Saxon distribution.

Specifically, finite nucleons increase the effective diffusivity $a$ (surface
thickness).  This script computes the change by convolving the radial
Woods-Saxon function with a Gaussian of width $w$, then fitting the resulting
convolution to a Woods-Saxon with corrected diffusivity $a'$.

After repeating the convolution and refitting process for several values of $a$
and nucleon size $w$, it is found that the corrected diffusivity is very nearly

    a'^2 = a^2 + c^2*w^2

where c ~ 0.61 is a universal constant independent of $a$ and $w$.

Hence, given a desired diffusivity $a'$ and nucleon width $w$, one should
sample nucleon positions from a Woods-Saxon distribution with diffusivity

    a = sqrt(a'^2 - c^2*w^2).
"""

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import integrate
from scipy import optimize


def woods_saxon(r, a, R=0):
    return 1/(1 + np.exp((r - R)/a))


def gaussian(x, w, x0=0):
    return 1/(np.sqrt(2*np.pi)*w) * np.exp(-0.5*np.square((x - x0)/w))


def correct_diffusivity(a, w, plot=False):
    r = np.linspace(-6*a, 6*a, 100)

    smeared_woods_saxon = np.array([
        integrate.quad(
            lambda x: woods_saxon(r_ - x, a) * gaussian(x, w),
            -10*a, 10*a
        )[0] for r_ in r
    ])

    a_corrected = optimize.curve_fit(
        woods_saxon, r, smeared_woods_saxon, p0=(a,)
    )[0][0]

    if plot:
        plt.plot(r, woods_saxon(r, a), label=a)
        plt.plot(r, smeared_woods_saxon, label='smeared')
        plt.plot(r, woods_saxon(r, a_corrected), label=a_corrected)
        plt.legend()
        plt.show()

    return a_corrected


def main():
    # correct_diffusivity(0.5, 0.5, plot=True)

    w = np.linspace(0.05, 0.95, 20)

    for a in np.arange(3, 7)/10:
        corrected = np.array([correct_diffusivity(a, w_) for w_ in w])

        def ansatz(x, coeff):
            return np.sqrt(a*a + coeff*coeff*x*x)

        coeff = optimize.curve_fit(ansatz, w, corrected)[0][0]
        X = np.linspace(0, 1, 1000)
        plt.plot(X, ansatz(X, coeff), color='0.5')
        print(a, coeff)

        plt.plot(w, corrected, 'o',
                 label='$a = {:.1f}$, $c = {:.4f}$'.format(a, coeff))

    plt.text(0.05, 0.71, r'$a^{\prime 2} = a^2 + c^2w^2$')
    plt.xlabel('Gaussian nucleon width $w$ [fm]')
    plt.ylabel('Corrected Woods-Saxon diffusivity $a\'$ [fm]')
    plt.xlim(0, 1)
    plt.legend(loc='best')

    plt.tight_layout(pad=0.2)
    plt.savefig('corrected-woods-saxon.pdf')


if __name__ == "__main__":
    main()

corrected-woods-saxon.png