On Colab the first import exfor_tools downloads and unpacks the EXFOR database, which takes a few minutes.
1# Google Colab setup: installs jitr and this notebook's extra dependencies.
2# Does nothing when the notebook is run anywhere else.
3import sys
4
5if "google.colab" in sys.modules:
6 import subprocess
7
8 subprocess.run(
9 [sys.executable, "-m", "pip", "install", "-q", "jitr", "exfor-tools"],
10 check=True,
11 )
Analyzing power of the quasi-elastic \((\vec{p},n)\) reaction¶
In 1976 Gosset, Mayer and Escudié measured the differential analyzing power \(A_y(\theta)\) of the quasi-elastic \((p,n)\) reaction to the isobaric analog state (IAS) on nine targets, with a 22.8 MeV polarized proton beam at Saclay (Phys. Rev. C 14, 878). Those data are EXFOR entry O0669.
Their point was that \(A_y\) is a far sharper probe of the Lane potential than the cross section is. In the macroscopic Lane model the nucleon-nucleus optical potential splits into isoscalar and isovector parts,
which for the proton-target, neutron-analog and neutron-target channels gives
with \(T_0 = (N - Z)/2\). The off-diagonal part of \(U_1\) drives the charge-exchange transition between \(|T_0 T_0\rangle\) and its analog \(|T_0, T_0 - 1\rangle\). Because the transition carries \(l = s = j = 0\), only two amplitudes survive — one conserving and one inverting the projectile spin projection — and the analyzing power is
The spin-flip amplitude \(Y\) is built from the difference between the \(j = l + \tfrac{1}{2}\) and \(j = l - \tfrac{1}{2}\) radial matrix elements.
This notebook solves the two coupled Lane channels exactly with
jitr.xs.lane_pn
rather than in DWBA, and propagates the posteriors of three uncertainty-quantified global optical
potentials — KDUQ, WLH and CHUQ — plus a deterministic JLMB calculation, through to \(A_y\) for all
nine targets.
Notw that, while \(^{49}\)Ti, \(^{117}\)Sn and \(^{165}\)Ho all have non-zero ground-state spin, the measured transitions transfer no angular momentum, so the target spin is a spectator and the spin-1/2 geometry used here is still reasonable.
1import os
2
3import matplotlib.pyplot as plt
4import numpy as np
5import periodictable
6
7import jitr
8from jitr.folding import ILDAFolder
9from jitr.folding.jlm import (
10 lambda_v0,
11 lambda_v1,
12 lambda_vso,
13 lambda_w0,
14 lambda_w1,
15 lambda_wso,
16 potential_JLMB,
17 spin_orbit_jlmb,
18)
19from jitr.optical_potentials import chuq, kduq, wlh
20from jitr.utils.density import density_table
21from jitr.xs import lane_pn
22
23# set JITR_QUICK=1 to thin the posterior draws; used by the CI notebook run
24QUICK = os.environ.get("JITR_QUICK", "0") == "1"
Pulling the measurements from EXFOR¶
Each of the nine targets is a subentry of O0669, tabulated as an analyzing power against the center-of-mass angle, together with the excitation energy of the analog state in the residual nucleus. That excitation energy is what fixes the outgoing neutron energy, so we read it from the entry to set the kinematics in the exit channels.
1import exfor_tools
2
3print(exfor_tools.__version__)
Using database version X4-2024-12-31 located in: /home/kyle/db/exfor/unpack_exfor-2024/X4-2024-12-31
1.3
1from exfor_tools import ExforEntry
2from exfor_tools.reaction import Reaction as ExforReaction
3
4# the nine targets of Table I, in the order the paper presents them
5TARGETS = [
6 (49, 22),
7 (56, 26),
8 (64, 28),
9 (70, 30),
10 (90, 40),
11 (96, 40),
12 (117, 50),
13 (165, 67),
14 (208, 82),
15]
16PROTON, NEUTRON = (1, 1), (1, 0)
17
18
19def nuclide(A, Z):
20 return rf"$^{{{A}}}${periodictable.elements[Z].symbol}"
21
22
23datasets = []
24for A, Z in TARGETS:
25 entry = ExforEntry(
26 "O0669",
27 ExforReaction(
28 target=(A, Z), projectile=PROTON, product=NEUTRON, residual=(A, Z + 1)
29 ),
30 quantity="Ay",
31 )
32 (measurement,) = entry.measurements
33 datasets.append(
34 {
35 "A": A,
36 "Z": Z,
37 "name": f"{A}{periodictable.elements[Z].symbol}",
38 "residual_name": f"{A}{periodictable.elements[Z + 1].symbol}",
39 "label": nuclide(A, Z),
40 "subentry": measurement.subentry,
41 "E_lab": measurement.Einc,
42 "E_IAS": measurement.Ex,
43 "theta": measurement.x,
44 "Ay": measurement.y,
45 "Ay_err": measurement.statistical_err,
46 }
47 )
48
49print(entry.bibtex())
@article{JGosset1976_O0669001,
author = {J.Gosset and B.Mayer and J.L.Escudie},
title = {Quasielastic (P,N) Reactions Induced By Polarized Protons.},
journal = {Physical Review, Part C, Nuclear Physics},
year = {1976},
volume = {14},
pages = {878},
note = {EXFOR reference: Physical Review, Part C, Nuclear Physics 14, 878 (1976); EXFOR subent: O0669001}
}
1header = f"{'target':>8} {'residual':>10} {'subentry':>10} {'E_IAS [MeV]':>12} {'points':>7} {'angular range [deg]':>22}"
2print(header)
3print("-" * len(header))
4for d in datasets:
5 angles = f"{d['theta'][0]:.1f} - {d['theta'][-1]:.1f}"
6 print(
7 f"{d['name']:>8} {d['residual_name']:>10} {d['subentry']:>10} "
8 f"{d['E_IAS']:>12.2f} {len(d['theta']):>7} {angles:>22}"
9 )
10print(
11 f"\n{sum(len(d['theta']) for d in datasets)} points, all at E_p = {datasets[0]['E_lab']} MeV"
12)
target residual subentry E_IAS [MeV] points angular range [deg]
--------------------------------------------------------------------------
49Ti 49V O0669002 5.36 18 16.6 - 103.2
56Fe 56Co O0669003 3.51 23 15.7 - 125.1
64Ni 64Cu O0669004 6.70 19 16.7 - 106.4
70Zn 70Ga O0669005 8.12 20 16.0 - 114.3
90Zr 90Nb O0669006 5.03 22 15.4 - 124.3
96Zr 96Nb O0669007 11.07 18 14.4 - 98.7
117Sn 117Sb O0669008 11.18 17 15.0 - 94.9
165Ho 165Er O0669009 15.49 8 16.9 - 61.0
208Pb 208Bi O0669010 15.33 4 17.3 - 31.4
149 points, all at E_p = 22.8 MeV
Compiling one coupled-channels workspace per target¶
jitr.xs.lane_pn.Workspace solves the two-channel problem
for every \((l, j)\), with an incoming wave in the proton channel only.
1E_LAB = 22.8
2LMAX = 25
3ANGLES = np.radians(np.linspace(1.0, 180.0, 180))
4
5solvers = {}
6for d in datasets:
7 A, Z = d["A"], d["Z"]
8 reaction = jitr.reactions.Reaction(
9 target=(A, Z), projectile=PROTON, product=NEUTRON, residual=(A, Z + 1)
10 )
11 # "El" so that the neutron potential is evaluated for n + residual
12 exit_reaction = jitr.reactions.Reaction(
13 target=(A, Z + 1), projectile=NEUTRON, process="El"
14 )
15 kinematics_entrance = reaction.kinematics(E_LAB)
16 kinematics_exit = reaction.kinematics_exit(
17 kinematics_entrance, residual_excitation_energy=d["E_IAS"]
18 )
19
20 a = jitr.utils.interaction_range(A) * kinematics_entrance.k + 2 * np.pi
21 nbasis = jitr.utils.suggested_basis_size(a)
22 solver = solvers.setdefault(nbasis, jitr.rmatrix.Solver(nbasis))
23
24 d.update(
25 reaction=reaction,
26 exit_reaction=exit_reaction,
27 kinematics_entrance=kinematics_entrance,
28 kinematics_exit=kinematics_exit,
29 nbasis=nbasis,
30 channel_radius_fm=a / kinematics_entrance.k,
31 workspace=lane_pn.Workspace(
32 reaction,
33 kinematics_entrance,
34 kinematics_exit,
35 solver,
36 ANGLES,
37 LMAX,
38 a / kinematics_entrance.k,
39 ),
40 )
41
42for d in datasets:
43 print(
44 f"{d['name']:>8} E_n = {d['kinematics_exit'].Elab:5.2f} MeV "
45 f"a = {d['channel_radius_fm']:5.2f} fm nbasis = {d['nbasis']:2d}"
46 )
49Ti E_n = 15.92 MeV a = 10.47 fm nbasis = 20
56Fe E_n = 13.78 MeV a = 10.66 fm nbasis = 20
64Ni E_n = 13.50 MeV a = 10.85 fm nbasis = 20
70Zn E_n = 13.11 MeV a = 10.99 fm nbasis = 20
90Zr E_n = 10.74 MeV a = 11.40 fm nbasis = 20
96Zr E_n = 10.99 MeV a = 11.52 fm nbasis = 20
117Sn E_n = 8.96 MeV a = 11.88 fm nbasis = 20
165Ho E_n = 6.05 MeV a = 12.58 fm nbasis = 25
208Pb E_n = 3.72 MeV a = 13.10 fm nbasis = 25
Propagating the optical model posteriors¶
Let’s propagate the KDUQ, WLH and CHUQ posteriors through the \((p,n)\) \(A_y\).
Nothing here is a fit to these data: the potentials were calibrated to elastic scattering, and the Lane transition potential is taken to be the isovector difference
evaluated separately for the central and spin-orbit terms — so the isovector spin-orbit form factor that the analyzing power is sensitive to comes entirely from the difference between the proton and neutron spin-orbit potentials of each global parameterization. Note that the neutron and proton potentials going into \(U_1\) are evaluated at the energies of their respective channels.
CHUQ deserves a word: CH89 is Lane-consistent by construction, so it is the one built-in whose proton and neutron potentials really are the \(U_{pT}\) and \(U_{nA}\) of the equations above, and whose difference therefore really is \(U_1\). For KDUQ and WLH the same construction is an approximation, and the diagonal potentials may already absorb part of the coupling.
1from tqdm import tqdm
2
3OMPS = {
4 "KDUQ": (
5 kduq.KDUQ(PROTON),
6 kduq.KDUQ(NEUTRON),
7 kduq.get_samples(PROTON),
8 kduq.get_samples(NEUTRON),
9 ),
10 "WLH": (
11 wlh.WLH(PROTON),
12 wlh.WLH(NEUTRON),
13 wlh.get_samples(PROTON),
14 wlh.get_samples(NEUTRON),
15 ),
16 # CH89 is Lane consistent: one parameter set serves both projectiles
17 "CHUQ": (chuq.CHUQ(), chuq.CHUQ(), chuq.get_samples(), chuq.get_samples()),
18}
19
20
21def propagate(omp_p, omp_n, samples_p, samples_n):
22 n_draws = min(len(samples_p), len(samples_n))
23 if QUICK:
24 n_draws = min(n_draws, 50)
25
26 bands = []
27 for d in tqdm(datasets):
28 workspace = d["workspace"]
29 rgrid = workspace.radial_grid()
30 Ay = np.zeros((len(ANGLES), n_draws))
31 dsdo = np.zeros((len(ANGLES), n_draws))
32 for j in range(n_draws):
33 U_p, U_p_so, U_coulomb = omp_p(
34 rgrid, d["reaction"], d["kinematics_entrance"], *samples_p[j]
35 )
36 U_n, U_n_so, _ = omp_n(
37 rgrid, d["exit_reaction"], d["kinematics_exit"], *samples_n[j]
38 )
39 _, S = workspace.rsmatrix(U_coulomb, U_p, U_p_so, U_n, U_n_so)
40 observables = workspace.observables_from_smatrix(S)
41 Ay[:, j] = observables.Ay
42 dsdo[:, j] = observables.dsdo
43 bands.append(
44 {
45 "Ay": np.percentile(Ay, [16, 84], axis=1),
46 "Ay_median": np.median(Ay, axis=1),
47 "dsdo": np.percentile(dsdo, [16, 84], axis=1),
48 "dsdo_median": np.median(dsdo, axis=1),
49 }
50 )
51 return n_draws, bands
1predictions = {}
2for name, args in OMPS.items():
3 n_draws, bands = propagate(*args)
4 predictions[name] = bands
5 print(f"{name}: {n_draws} posterior draws x {len(datasets)} targets")
0%| | 0/9 [00:00<?, ?it/s]/home/kyle/umich/jitr/src/jitr/optical_potentials/kduq.py:470: RuntimeWarning: overflow encountered in exp
d2 = d2_0 + d2_A / (1 + np.exp((A - d2_A3) / d2_A2))
100%|█████████████████████████████████████████████████████████████████| 9/9 [00:26<00:00, 2.90s/it]
KDUQ: 416 posterior draws x 9 targets
100%|█████████████████████████████████████████████████████████████████| 9/9 [00:45<00:00, 5.01s/it]
WLH: 1000 posterior draws x 9 targets
100%|█████████████████████████████████████████████████████████████████| 9/9 [00:09<00:00, 1.04s/it]
CHUQ: 208 posterior draws x 9 targets
A semi-microscopic comparison: JLMB¶
JLMB is a different kind of potential: a local-density folding of the Jeukenne-Lejeune-Mahaux nuclear-matter g-matrix over a Gogny D1M density, with energy-dependent normalizations \(\lambda_V, \lambda_W, \lambda_{V_1}, \lambda_{W_1}\) and a Scheerbaum spin-orbit form factor. It is the one potential here with a genuine isovector form factor rather than one obtained by differencing.
It is not uncertainty quantified, so it contributes a single curve, not a band.
1def jlmb_potentials(d):
2 workspace = d["workspace"]
3 rgrid = workspace.radial_grid()
4 folder = ILDAFolder(r_max=15.0, n_quad=200)
5
6 def channel(A, Z, projectile, energy, coulomb):
7 table = density_table(A, Z, model="d1m")
8 rho_n = folder.interp_to_quad(table.radial_grid, table.neutron_density_grid)
9 rho_p = folder.interp_to_quad(table.radial_grid, table.proton_density_grid)
10 V_C_quad = folder.V_coulomb(rho_p) if coulomb else None
11 V, W = potential_JLMB(
12 folder,
13 rho_n,
14 rho_p,
15 projectile,
16 (A, Z),
17 energy,
18 V_C=V_C_quad,
19 parameterization="talys",
20 lambda_V=float(lambda_v0(energy)),
21 lambda_W=float(lambda_w0(energy, mode=0)),
22 lambda_V1=float(lambda_v1(energy)),
23 lambda_W1=float(lambda_w1(energy, mode=0)),
24 t_r=1.25,
25 t_i=1.35,
26 r_out=rgrid,
27 )
28 spin_orbit = (
29 float(lambda_vso(energy)) + 1j * float(lambda_wso(energy))
30 ) * spin_orbit_jlmb(folder.r_q, rho_n, rho_p, projectile, r_out=rgrid)
31 coulomb_out = folder.V_coulomb(rho_p, r_out=rgrid) if coulomb else None
32 return (
33 np.asarray(V + 1j * W, dtype=np.complex128),
34 np.asarray(spin_orbit, dtype=np.complex128),
35 (
36 None
37 if coulomb_out is None
38 else np.asarray(coulomb_out, dtype=np.complex128)
39 ),
40 )
41
42 U_p, U_p_so, U_coulomb = channel(
43 d["A"], d["Z"], PROTON, d["kinematics_entrance"].Elab, coulomb=True
44 )
45 U_n, U_n_so, _ = channel(
46 d["A"], d["Z"] + 1, NEUTRON, d["kinematics_exit"].Elab, coulomb=False
47 )
48 return U_coulomb, U_p, U_p_so, U_n, U_n_so
49
50
51jlmb = []
52for d in datasets:
53 jlmb.append(d["workspace"].observables(*jlmb_potentials(d)))
54print(f"JLMB evaluated for {len(jlmb)} targets")
JLMB evaluated for 9 targets
Plot the measurements against the calculations¶
1BANDS = {"KDUQ": "#0072B2", "WLH": "#D55E00", "CHUQ": "#009E73"}
2JLMB_INK = "#333333"
3angles_deg = np.degrees(ANGLES)
4
5
6def grid_figure(ylabel):
7 fig, axes = plt.subplots(3, 3, figsize=(12, 10), sharex=True)
8 fig.subplots_adjust(hspace=0.12, wspace=0.22)
9 for ax in axes[-1]:
10 ax.set_xlabel(r"$\theta_{\rm c.m.}$ [deg]")
11 for ax in axes[:, 0]:
12 ax.set_ylabel(ylabel)
13 for ax in axes.flat:
14 ax.set_xlim(0, 140)
15 ax.tick_params(direction="in", top=True, right=True)
16 for side in ax.spines.values():
17 side.set_linewidth(0.8)
18 return fig, axes
19
20
21fig, axes = grid_figure(r"$A_y$")
22for i, (ax, d, curve) in enumerate(zip(axes.flat, datasets, jlmb, strict=True)):
23 for name, color in BANDS.items():
24 lo, hi = predictions[name][i]["Ay"]
25 ax.fill_between(angles_deg, lo, hi, color=color, alpha=0.35, lw=0, label=name)
26 # ax.plot(angles_deg, predictions[name][i]["Ay_median"], color=color, lw=1.2)
27 ax.plot(angles_deg, curve.Ay, color=JLMB_INK, lw=1.2, ls="--", label="JLMB")
28 ax.errorbar(
29 d["theta"],
30 d["Ay"],
31 yerr=d["Ay_err"],
32 fmt="o",
33 ms=3.5,
34 color="k",
35 lw=0.9,
36 capsize=1.5,
37 label="Gosset (1976)",
38 )
39 ax.axhline(0.0, color="0.75", lw=0.7, zorder=0)
40 ax.set_ylim(-1.05, 1.05)
41 ax.text(0.04, 0.06, d["label"], transform=ax.transAxes, fontsize=13)
42
43axes[0, 0].legend(loc="upper right", fontsize=8, frameon=False, ncols=2)
44fig.suptitle(r"$(\vec{p},n)$ analyzing power to the IAS at $E_p = 22.8$ MeV", y=0.92)
45plt.show()
