jitr quickstart: fitting an optical potential to \(\alpha\) + \(^{44}\)Ca elastic scattering at 29 MeV

jitr solves the scattering problem for a given nucleus-nucleus interaction. In this first example we

• set up a reaction and its kinematics

• compile a solver for it

• define a parametric optical potential and calculate an elastic differential cross section

• fit the potential parameters to experimental data

The data come from EXFOR F0567 (Oeschler et al., 1972). We will look at the example of of the \(\alpha\)-elastic scattering on \(^{44}\)Ca, fitting an optical potential to differential elastic scattering cross sections (as a ratio to the Rutherford scattering cross section).

A subsequent demo, alpha_ca_calibration.ipynb does a full Bayesian calibration with a simple error model learned from the data.

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

Set up the reaction

A Reaction holds the target, projectile and process. Its kinematics method turns a lab energy into the quantities the solver needs: the center-of-mass energy, the wavenumber \(k\), the reduced mass, \(\mu\), and the Sommerfeld parameter \(\eta\).

from jitr.reactions.reaction import Reaction

reaction = Reaction(target=(44, 20), projectile=(4, 2), process="El")
energy_lab = 29  # MeV
kinematics = reaction.kinematics(energy_lab)
print(reaction)
print(
    f"Ecm: {kinematics.Ecm:1.2f} [MeV]\nk: {kinematics.k:1.2f} [fm^-1]\nmu: {kinematics.mu:1.2f}\neta: {kinematics.eta:1.2f}"
)

44-Ca(alpha,el)
Ecm: 26.58 [MeV]
k: 2.16 [fm^-1]
mu: 3437.68
eta: 2.35

Read in the data

df = pd.read_csv("./alpha_ca_ratio_ruth.csv")
df = df[df["target"] == "Ca44"]
angles_deg = df["angle_cm_deg"].to_numpy()
y_data = df["ratio_to_rutherford"].to_numpy()
print(
    f"{len(y_data)} data points between {angles_deg.min():.0f} and {angles_deg.max():.0f} degrees"
)


def plot_data(ax, y_err=None):
    if y_err is not None:
        ax.errorbar(
            angles_deg,
            y_data,
            y_err,
            marker=".",
            linestyle="none",
            color="k",
            alpha=0.5,
            label="data",
        )
    else:
        ax.plot(angles_deg, y_data, ".", color="k", alpha=0.5, label="data")
    ax.set_yscale("log")
    ax.set_xlabel(r"$\theta_{cm}$ [deg]")
    ax.set_ylabel(
        r"$(\frac{d\sigma}{d\Omega}) / (\frac{d\sigma}{d\Omega})_\text{Rutherford}$"
    )
    ax.set_title(r"$^{44}$Ca($\alpha,\alpha$) at 29 MeV")


fig, ax = plt.subplots(figsize=(6, 4))
plot_data(ax)
plt.tight_layout()
plt.show()

289 data points between 15 and 174 degrees

Compile the solver

jitr uses the \(\mathcal{R}\)-matrix method on a Lagrange-Legendre mesh following Descouvemont, 2016. The solver needs a channel radius (beyond which the nuclear potential is negligible and the wavefunction is matched to Coulomb functions), a number of basis functions \(N\), and a maximum partial wave \(\ell_\text{max}\). The helpers below give rules of thumb for the first two. For the rather sharp-edged \(\alpha\) potentials used here, the default number of nodes is not quite converged at backward angles, so we ask for 8 nodes per half-wavelength (see convergence_channel_radius.ipynb for how to check this yourself).

The DifferentialWorkspace precomputes everything that does not depend on the potential for a fixed set of angles - here the angles of the data. Building it compiles code and takes a few seconds. Once it is compiled, it is fast to pass in an interaction and get back a cross section.

from jitr.rmatrix import Solver
from jitr.utils import utils
from jitr.xs import elastic

interaction_range_fm = (
    1.2 * (reaction.target.A ** (1 / 3) + reaction.projectile.A ** (1 / 3)) + 2
)
channel_radius_dimensionless = utils.suggested_dimensionless_channel_radius(
    interaction_range_fm, kinematics.k
)
channel_radius_fm = channel_radius_dimensionless / kinematics.k
N = utils.suggested_basis_size(channel_radius_dimensionless, zeros_per_node=8)
lmax = 60
print(
    f"channel radius: {channel_radius_fm:.2f} [fm]\nSystem size per symmetry block {N}x{N}\nMax partial wave: {lmax}"
)

workspace = elastic.DifferentialWorkspace.build_from_system(
    reaction=reaction,
    kinematics=kinematics,
    channel_radius_fm=channel_radius_fm,
    solver=Solver(N),
    lmax=lmax,
    angles=np.deg2rad(angles_deg),
)

channel radius: 11.05 [fm]
System size per symmetry block 64x64
Max partial wave: 60

Define the interaction

The potential is a function radial coordinate – the magnitude of the distance between the projectile and target center of mass – that we evaluate on the workspace’s quadrature grid. We use a real volume Woods-Saxon shaped potential, an imaginary volume and an imaginary surface term, and the Coulomb potential of a uniformly charged sphere. Seven parameters are free: the depths \(V\), \(W_v\), \(W_s\) and the geometry \(r_v, a_v, r_w, a_w\) (radii in units of \(A^{1/3}\) fm; the surface term shares the imaginary geometry). workspace.xs returns the differential cross section; dividing by workspace.rutherford gives the observable in the data.

The parameteric Woods-Saxon shape is

(1)\[\begin{equation} f(r;R,a) = \left[ 1 + \exp{\left( \frac{r-R}{a} \right)} \right]^{-1}, \end{equation}\]

the surface-peaked shape is just the derivative of this with respect to \(r\).

from jitr.optical_potentials.potential_forms import (
    coulomb_charged_sphere as coulomb,
)
from jitr.optical_potentials.potential_forms import (
    woods_saxon_prime_safe as wsp,
)
from jitr.optical_potentials.potential_forms import (
    woods_saxon_safe as ws,
)


def U_central(r, Vv, Wv, Ws, Rv, av, Rw, aw):
    return (
        -Vv * ws(r, Rv, av)
        - 1j * Wv * ws(r, Rw, aw)
        - 1j * Ws * (-4 * aw) * wsp(r, Rw, aw)
    )


rgrid = workspace.radial_grid()
A_factor = reaction.target.A ** (1 / 3)
Zz = reaction.target.Z * reaction.projectile.Z
coulomb_potential = coulomb(rgrid, Zz, 1.3 * A_factor)


def calculate_xs_ratio(theta):
    """Ratio to Rutherford at the data angles for optical model parameters theta."""
    Vv, Wv, Ws, rv, av, rw, aw = theta
    central_potential = U_central(
        rgrid, Vv, Wv, Ws, rv * A_factor, av, rw * A_factor, aw
    )
    xs = workspace.xs(
        central_potential=central_potential, coulomb_potential=coulomb_potential
    )
    return xs.dsdo / workspace.rutherford


names = ["Vv", "Wv", "Ws", "rv", "av", "rw", "aw"]

A typical \(\alpha\)-nucleus potential at these energies has \(V \approx 185\) MeV, \(r_v \approx 1.4\) fm, \(a \approx 0.5\) fm and \(W \approx 25\) MeV.

theta_0 = np.array([185, 25, 25, 1.4, 0.5, 1.4, 0.4])


def plot_potential(theta, ax):
    Vv, Wv, Ws, rv, av, rw, aw = theta
    U = U_central(rgrid, Vv, Wv, Ws, rv * A_factor, av, rw * A_factor, aw)
    ax.plot(rgrid, U.real, "k", label=r"$\mathfrak{Re} U(r)$")
    ax.plot(rgrid, U.imag, "--k", label=r"$\mathfrak{Im} U(r)$")
    ax.legend(ncol=2, framealpha=0)
    ax.set_xlabel(r"$r$ [fm]")
    ax.set_ylabel(r"U(r) [MeV]")
    ax.set_title(r"$^{44}$Ca($\alpha,\alpha$) at 29 MeV")


fig, ax = plt.subplots(figsize=(6, 4))
plot_potential(theta_0, ax)

Let’s see how that does.

fig, ax = plt.subplots(figsize=(6, 4))
plot_data(ax)
ax.plot(angles_deg, calculate_xs_ratio(theta_0), color="#3366CC", label=r"$\theta_0$")
ax.legend(framealpha=0)
plt.tight_layout()
plt.show()

Not a bad start.

Fit the potential to the data

The paper gives no uncertainties, so we have to assume something. The simplest reasonable choice is an unknown relative error on every point, \(\sigma_i = \beta\, y_i^{\text{pred}}\), with the same fraction \(\beta\) everywhere.

We will use scipy.optimize.least_squares for fitting, to find the parameters \(\alpha\) that minimize the \(\chi^2\):

(2)\[\begin{equation} \chi^2(\alpha) \equiv \sum_i \frac{ (y_i - y_m(x_i;\alpha) )^2}{(\beta y_i)^2}. \end{equation}\]

The minimum of this function is \(\beta\)-independent, so this does not require making any assumptions other than that the error is relative.

%%time
from scipy.optimize import least_squares

theta_bounds = np.array(
    [[120, 250], [0, 60], [0, 60], [1.15, 1.6], [0.4, 0.8], [1.0, 1.8], [0.2, 0.9]]
)


def relative_residuals(theta):
    y_pred = calculate_xs_ratio(theta)
    return (y_data - y_pred) / y_pred


fit = least_squares(relative_residuals, theta_0, bounds=theta_bounds.T, x_scale="jac")
theta_mle = fit.x
error_fraction = np.sqrt(np.mean(fit.fun**2))
print(f"converged: {fit.success} after {fit.nfev} evaluations")
print(f"inferred relative error: {100 * error_fraction:.0f}%")

converged: True after 32 evaluations
inferred relative error: 53%
CPU times: user 24 s, sys: 11.7 ms, total: 24 s
Wall time: 3.1 s

pd.DataFrame({"start": theta_0, "fit": theta_mle}, index=names).T.round(3)

	Vv	Wv	Ws	rv	av	rw	aw
start	185.000	25.000	25.000	1.400	0.500	1.400	0.400
fit	215.744	21.822	0.104	1.315	0.577	1.694	0.332

fig, ax = plt.subplots(figsize=(6, 4))
plot_data(ax, y_err=0.1 * y_data)
ax.plot(
    angles_deg,
    calculate_xs_ratio(theta_0),
    color="#3366CC",
    alpha=0.4,
    label=r"$\theta_0$",
)
ax.plot(
    angles_deg,
    calculate_xs_ratio(theta_mle),
    color="#3366CC",
    linewidth=2,
    label="maximum likelihood",
)
ax.legend(framealpha=0)
plt.tight_layout()
plt.show()

What did we learn?

The fit roughly reproduces the angular spacing of diffraction pattern as well as the general rate at which the ratio decreases with angle (which tells us something about the amplitude for non-elastic processes), but the depth of each of the dips is poorly reproduced after about 20-30\(^\circ\).

To make matters worse, light ion-nucleus potentials are famously ambiguous - several distinct families of parameters, differing in \(V\) by tens of MeV, give nearly the same cross section. Where the optimizer lands depends on where it starts.

Finally, we have said nothing about the uncertainty of our predictions. Surely, with a best fit line that does not pass through all of our data, we would expect the corresponding uncertainty in the inferred parameters to actually be large.

The subsequent notebook, alpha_ca_calibration.ipynb does a full Bayesian calibration in which the scales of the errors are learned from the data alongside the potential, and checks the result against the data.