phasr is a python package able to calculate scattering phase shifts for arbitrary radial potentials and calculate (elastic) electron nucleus scattering cross sections based on it.
In combination with a fitting routine it can be used to fit elastic electron nucleus scattering data to extract charge density parameterizations. It is also possible to increase the precision sufficiently to quantitatively resolve parity violating electron scattering.
Version 0.4.0 : Add left-right asymmetry.
Version 0.3.2 : Extend docs.
Version 0.3.1 : Add access to reference parameterizations.
Version 0.3.0 : Add overlap integrals.
Version 0.2.2 : Fix some dependencies.
Version 0.2.1 : First PyPI Release.
Version 0.2.0 : Basic features implemented.
Version 0.1.4 : Build dirac_solvers module.
Version 0.1.3 : Build nuclei module.
Version 0.1.2 : Build skeleton.
Version 0.1.0 : First init. Empty.
This package was developed with the goal of extracting charge density parameterizations based on (elastic) electron nucleus scattering cross section data while properly accounting for Coulomb corrections using the so called phase shift model. While this procedure is known since the 60s, the documentation and availability of programs used back then, are very limited and hardly accessible. Furthermore, for a majority of nuclei the extractions of model-independent parameterizations of their charge densities have not been updated since the 80s and do not provide any uncertainty estimates. Since these also applied to the nuclei we were interested in, we build such a framework in python, which calculates elastic electron nucleus scattering cross sections based on input charge densities using the phase shift model, from which this package evolved. With such a more modern implementation and programming language, the calculations can become faster, broader accessible, easier to understand and easier to incorporate into existing programs. The code leading up to this package has already been used successfully in the papers [Noël, Hoferichter - 2024] and [Heinz, Hoferichter, Miyagi, Noël, Schwenk - 2024] to calculate and extract charge distributions for nuclei of interest to \( \mu \to e \) conversion, a beyond the Standard Model process of interest in particle physics. The package has been benchmarked for nuclei reaching from ... to ... .
This program is built with python 3.
It requires the following packages to operate properly:
numpy
scipy
mpmath
To install and use the program please follow the following steps for setup. We assume you already have a working python environment setup on your machine (othewise see REF).
The package is listed on PyPI, which is why it can be easily installed via:
pip install phasr You can import the python package the same way as every other package, i.e.:
import phasr as phr nucleus_Al27 = phr.nucleus(name='Al27',Z=13,A=27,rc=5)
mmu = phr.masses.mmu
boundstates_Al27 = phr.boundstates(nucleus_Al27,kappa=-1,lepton_mass=mmu)
print('E_1s=',boundstates_Al27.energy_levels[0],'MeV') phr throughout this manual.
The code is structured into two main modules.
The nuclei module is home of the nucleus class and manages all the relevant properties of a nucleus/isotope that is considered and calculates other properties based on the information it has available (e.g. the electric field from the charge density).
The dirac_solvers module is home of the boundstates and continuumstates classes and manages solving the Dirac equation for a given nucleus based on its nucleus potential and returns the desired consequential quantities like for example electron wave functions, phase shifts or the crosssection.
In all applicational scenarios one will first construct a nucleus object, which is then be used to construct a boundstates or continuumstates object. For calculating elastic scattering crosssections we supply the function crosssection_lepton_nucleus_scattering,
which can be used to directly calculate the crosssection using the energy, the angle and a nucleus object.
nucleus class To get a nucleus parameterization, use the nucleus class to construct a nucleus object.
A nucleus object requires a few minimum parameters (e.g. a name, a charge \(Z\) and a nucleon number \(A\)) and and has a lot of optional parameters (e.g. some parameters for a specific parameterization).
Based on the parameters supplied the class decides what kind of nucleus it is supposed to construct. You can check what parameterization was choosen based on the attribute nucleus_type
phr.nucleus(name,Z,A,m=None,**args) | nucleus_type | Parameters / **args |
|---|---|
| 'coulomb' | None |
| 'fourier-bessel' | ai (1d-array): Array of parameters \(a_i\) in fm\(^{-3}\). \(N\) given by length of a.
R (float): Cutoff radius R in fm. |
| 'oszillator-basis' | Ci_dict (dict): dictionary of 1d-arrays \(C_i\), with multipole names as keys (see below). |
| 'fermi' | c (float): \(c\) parameter in fm for Fermi parameterization.
z (float): \(z\) parameter in fm for Fermi parameterization. w (float,optional): \(w\) parameter for Fermi parameterization (default: w=0). |
| 'gauss' | b (float): \(b\) parameter in fm for Gauss parameterization. |
| 'uniform' | rc (float): \(r_c\) parameter in fm for Uniform parameterization. |
| 'numerical' | see below |
nucleus object can have the following attributes (if the prerequisites are met), that you might want to access
| Parameters/**args | |
|---|---|
| charge_density (callable,optional): | Function describing the charge density (see below). |
| form_factor (callable,optional): | Function describing the form factor (see below). |
| electric_field (callable,optional): | Function describing the electric field (see below). |
| electric_potential (callable,optional): | Function describing the electric potential (see below). |
| form_factor_dict (dict,optional): | Dictionary containing form factor functions for specific multipoles (keys are 'F'+multipole , see below). |
| density_dict (dict,optional): | Dictionary containing density functions for specific multipoles (keys are 'rho'+multipole , see below). |
| rrange (list,optional): | Sets range parameters for the beginning of the integration and highenergy continuation (in fm, default: [0.,20.,0.02]). |
| qrange (list,optional): | Sets range parameters for the beginning of the integration and highenergy continuation (in MeV, default: [0.,1000.,1.]). |
| renew (bool,optional): | Says if numerical calculations should be done from scratch (aka renewed) or recovered from file if existent (default: False). |
multipoles. Fch(q,L=0), Fw(q,L=0), Fmag(q,L=1), which make it possible to also access the form factors for higher \(L\) if the necessary multipoles were supplied.
dirac_solvers moduleThe dirac_solvers module has two main duties, covered by two classes.
The boundstates class calculates the boundstate solutions of a given nucleus by solving the Dirac equation in the presence of the nucleus potential, which includes finding the boundstate energy first.
The continuumstates class calculates the continuumstate solution for a given energy of a given nucleus by solving the Dirac equation in the presence of the nucleus potential.
Finally, there are some functions that make use of the boundstate and continuumstate solutions by calculating elastic electron nucleus scattering crosssections and overlap integrals for \(\mu \to e\) conversion.
The boundstate and continuumstate wavefunctions are both parameterized according to
\[ \psi = \psi_\kappa^\mu(\vec{r})= \frac{1}{r} \Bigg(\begin{array}{c}
g_\kappa(r) \phi_\kappa^\mu(\hat{r}) \\ i f_\kappa(r) \phi_{-\kappa}^\mu(\hat{r}) \end{array} \Bigg), \]
separating off the angular degrees of freedom, which are the same for all radial symmetric potential. The \(\kappa\) label refers to the angular momentum quantum numbers and spin according to
\[j = |\kappa| - \frac{1}{2}; \quad l = \begin{cases} \kappa, & \quad \kappa > 0 \\ - \kappa - 1 , & \quad \kappa < 0 \end{cases}\]
The numerical calculations then provide descriptions for \(g_\kappa(r)\) and \(f_\kappa(r)\).
The boundstates class:
The boundstates class is used to calculate the boundstates for a given nucleus object, beginning from the ground state for the given partial wave.
First the energy for the boundstate is found by iteratively scanning for flips in the asymptotic behaviour of the wave functions, which indicate a boundstate in between, until a set precision is reached.
For this energy then the solution of the Dirac equation gives the dirac wavefunction for that state.
phr.boundstates(nucleus,kappa,lepton_mass,**args) phr.masses.me or muon mass phr.masses.mmu. needs to be non-zero.
scipy.integrate.solve_IVP). Default values are set by phr.default_boundstate_settings.
boundstates object has the following attributes, which you might want to access:
find_next_energy_level and solve_IVP_at_current_energy consecutively. Can be passed **args to change the solver_setting.
find_next_solution() was called enough times.
Example: Groundstate is wavefunction_g_1s12 (using \(\kappa=-1\))
find_next_solution() was called enough times.
Example: Groundstate is wavefunction_f_1s12 (using \(\kappa=-1\))
continuumstate class: phr.continuumstates(nucleus,kappa,energy,lepton_mass=0,**args) scipy.integrate.solve_IVP). Default values are set by phr.default_continuumstate_settings.
continuumstates object has the following attributes, which you might want to access:
solve_IVP() was called.
solve_IVP() was called.
solve_IVP().
crosssection_lepton_nucleus_scattering function: phr.crosssection_lepton_nucleus_scattering(energy,theta,nucleus,lepton_mass=0,recoil=True,subtractions=3,**args) optimise_crosssection_precision() (see below).
| Parameters/**args | |
|---|---|
| N_partial_waves (int,optional): | angular momenutum of the partial wave at which the series is terminated. (default: 100) |
| phase_difference_limit (float,optional): | limit at which (numerical) value for the phase_difference attribute of a continuumstates all following partial waves are considered indistinguishable to the coulomb solution. (default:0). |
| energy_norm (float,optional): | scaling factor for the units while solving the inital value problem to increase speed (default: \( \hbar c\) fm\(^{-1}\)) |
| method (str,optional): | method keyword passed to scipy.integrate.solve_IVP (default: 'DOP853' / as set in default_continuumstate_settings) |
| atol (float,optional): | atol keyword passed to scipy.integrate.solve_IVP (default: 1e-12 / as set in default_continuumstate_settings) |
| rtol (float,optional): | rtol keyword passed to scipy.integrate.solve_IVP (default: 1e-9 / as set in default_continuumstate_settings) |
crosssection_lepton_nucleus_scattering many times for the same energies and angles theta for similar nuclei (for example for a fit),
it makes sense to assess the necessary precision for the numerical calculation using optimise_crosssection_precision.
This functions can run a few minutes as it is scanning a large **args / parameter space, but can afterwards improve the runtime significantly.
Unfortunately the necessary precision is very dependent on how large the energy and angles are (larger energies and angles require significantly more precision and runtime),
which is why this is the best solution for a reasonable runtime. Of course you can always also try to finetune the parameters yourself, but note that the convergence might break down,
and that you should always check if the calculated crosssection is still reasonable.
phr.optimise_crosssection_precision(energy,theta,nucleus,lepton_mass=0,recoil=True,subtractions=3,crosssection_precision=1e-3,jump_forward_dist=1) crosssection_lepton_nucleus_scattering which are still consistent with crosssection calculations with the maximum precision up to a relative precision of the value of crosssection_precision.
To calculate the scattering cross section base on a analytic charge distribution parametrisation, you first construct the nucleus object, which you then use to calculate the cross section. As an example, we use a Fourier Bessel parameterisation for \(^{48}\)Ti with \( R = 9.25\) fm and \[a_{1} = 0.03392 ~ \text{fm}^{-3}, ~ a_{2}=0.05913 ~ \text{fm}^{-3}, ~ a_{3}=0.01547 ~ \text{fm}^{-3}, ~ a_{4}=-0.02550 ~ \text{fm}^{-3}, \] \[a_{5}=-0.0152 ~ \text{fm}^{-3}, ~ a_{6}=0.0029 ~ \text{fm}^{-3}, ~ a_{7}=0.0037 ~ \text{fm}^{-3}, \] The nucleus is then constructed by
ai_48Ti=np.array([0.03392,0.05913,0.01547,-0.02550,-0.0152,0.0029,0.0037])
R_48Ti=9.25
nucleus_48Ti = phr.nucleus('48Ti_FB',Z=22,A=48,ai=ai_48Ti,R=R_48Ti)import matplotlib.pyplot as plt
energy=249.5
theta=np.arange(30,120,1e-1)
plt.plot(theta,phr.crosssection_lepton_nucleus_scattering(energy,theta*pi/180,nucleus_48Ti),label=r'$E=249.5~$MeV')
plt.yscale('log')
plt.title(r"Cross section $^{48}$Ti")
plt.xlabel(r"$\theta$ in deg")
plt.ylabel(r"$\frac{\operatorname{d}\sigma}{\operatorname{d}\Omega}$ in MeV$^{-2}$")
plt.legend()
plt.xlim(30,120)
args=phr.optimise_crosssection_precision(energy,theta*pi/180,nucleus_48Ti,crosssection_precision=1e-2)import matplotlib.pyplot as plt
energy=249.5
theta=np.arange(30,120,1e-1)
plt.plot(theta,phr.crosssection_lepton_nucleus_scattering(energy,theta*pi/180,nucleus_48Ti,**args),label=r'$E=249.5~$MeV')
plt.yscale('log')
plt.title(r"Cross section $^{48}$Ti")
plt.xlabel(r"$\theta$ in deg")
plt.ylabel(r"$\frac{\operatorname{d}\sigma}{\operatorname{d}\Omega}$ in MeV$^{-2}$")
plt.legend()
plt.xlim(30,120)Calculating the scattering crosssection for a numerical charge density is for the most part identical to the analytical part. The only difference is that the nucleus object requires slighlty more preparation, as for example the nucleus potential needs to be calculated numerically. Let's construct the same nucleus as in the analytical case as a numerical nucleus. We do
ai_48Ti=np.array([0.03392,0.05913,0.01547,-0.02550,-0.0152,0.0029,0.0037])
R_48Ti=9.25
N_48Ti=len(ai_48Ti)
qi_48Ti = np.arange(1,N_48Ti+1)*np.pi/R_48Ti
def rho_48Ti(r): return phr.nuclei.parameterizations.fourier_bessel.charge_density_FB(r,ai_48Ti,R_48Ti,qi_48Ti)
nucleus_48Ti_num = phr.nucleus('48Ti_num',Z=22,A=48,charge_density=rho_48Ti)nucleus_48Ti_num.set_electric_field_from_charge_density()
nucleus_48Ti_num.set_electric_potential_from_electric_field()nucleus_48Ti_num.fill_gaps()energy=249.5
theta=np.arange(30,120,1e-1)
plt.plot(theta,phr.crosssection_lepton_nucleus_scattering(energy,theta*pi/180,nucleus_48Ti_num),label=r'$E=249.5~$MeV')
plt.yscale('log')
plt.title(r"Cross section $^{48}$Ti")
plt.xlabel(r"$\theta$ in deg")
plt.ylabel(r"$\frac{\operatorname{d}\sigma}{\operatorname{d}\Omega}$ in MeV$^{-2}$")
plt.legend()
plt.xlim(30,120)
plt.show()TBD
The package is currently in an development stage, and partly incomplete.
Use the package at your own risk and without any garanties.
This package was written by Frederic Noël.
I would like to thank my collaborators and office mates for meaningfull inputs. This website was build based on a template from github user charlyllo.