Class nucleus_rmf (o2scl)¶

class nucleus_rmf¶

Spherical closed-shell nuclei with a relativistic mean-field model in the Hartree approximation.

This code is very experimental.

This class is based on a code developed by C.J. Horowitz and B.D. Serot, and used in [Horowitz81] which was then adapted by P.J. Ellis and used in [Heide94] and [Prakash94ns]. Ellis and A.W. Steiner adapted it for the parameterization in eos_had_rmf for [Steiner05b], and then converted to C++ by Steiner afterwards. This code was also used in [Steiner13cs] to create SFHo and SFHx .

The standard usage is something like:

nucleus_rmf rn;
o2scl_hdf::rmf_load(rn.rmf,"NL4");
rn.run_nucleus(82,208,0,0);
cout << rn.rnrp << endl;

which computes the structure of \( ^{208}\mathrm{Pb} \) and outputs the neutron skin thickness using the model 'NL4'.

This code does not always converge, and this can sometimes be addressed by choosing a different initial guess, stored in nucleus_rmf::ig, which is an object of type nucleus_rmf::initial_guess .

Potential exceptions are

Failed to converge
Failed to solve meson field equations
Energy not finite (usually a problem in the equation of state)
Energy not finite in final calculation
Function iterate() called before init_run()
Not a closed-shell nucleus

The initial level pattern is

S 1/2
// 2 nucleons
P 3/2
P 1/2
// 8 nucleus
D 5/2
D 3/2
S 1/2
// 20 nucleons
F 7/2
// 28 nucleons
F 5/2
P 3/2
P 1/2
// 40 nucleons
G 9/2
// 50 nucleus
G 7/2
D 5/2
H 11/2
D 3/2
S 1/2
// 82 nucleons
H 9/2
F 7/2
I 13/2
F 5/2
P 3/2
P 1/2
// 126 nucleons
G 9/2
I 11/2
J 15/2
D 5/2
S 1/2
G 7/2
D 3/2
// 184 nucleons

Details

Below, \( \alpha \) is a generic index for the isospin, the radial quantum number \( n \) and the angular quantum numbers \( \kappa \) and \( m \). The meson fields are \( \sigma(r), \omega(r) \) and \( \rho(r) \). The baryon density is \( n(r) \), the neutron and proton densities are \( n_n(r) \) and \( n_p(r) \), and the baryon scalar density is \( n_s(r) \). The nucleon field equations are

\[\begin{split}\begin{eqnarray*} F^{\prime}_{\alpha}(r)- \frac{\kappa}{r} F_{\alpha}(r) + \left[ \varepsilon_{\alpha} - g_{\omega} \omega(r) - t_{\alpha} g_{\rho} \rho(r) - \left( t_{\alpha}+\frac{1}{2} \right) e A(r) - M + g_{\sigma} \sigma(r) \right] G_{\alpha}(r) &=& 0 \\ G^{\prime}_{\alpha}(r)+ \frac{\kappa}{r} G_{\alpha}(r) - \left[ \varepsilon_{\alpha} - g_{\omega} \omega(r) - t_{\alpha} g_{\rho} \rho(r) - \left( t_{\alpha}+\frac{1}{2} \right) e A(r) + M - g_{\sigma} \sigma(r) \right] F_{\alpha}(r) &=& 0 \end{eqnarray*}\end{split}\]

where the isospin number, \( t_{\alpha} \) is \( 1/2 \) for protons and \( -1/2 \) for neutrons. The meson field equations are

\[\begin{split}\begin{eqnarray*} \sigma^{\prime \prime}(r) + \frac{2}{r} \sigma^{\prime}(r) - m_{\sigma}^2 \sigma &=& - g_{\sigma} n_s(r) + b M g_{\sigma}^3 \sigma^2 + c g_{\sigma}^4 \sigma^3 - g_{\rho}^2 \rho^2 \frac{\partial f}{\partial \sigma} \\ \omega^{\prime \prime}(r) + \frac{2}{r} \omega^{\prime}(r) - m_{\omega}^2 \omega &=& - g_{\omega} n(r) + \frac{\zeta}{6} g_{\omega}^4 \omega^3 + g_{\rho}^2 \rho^2 \frac{\partial f}{\partial \omega} \\ \rho^{\prime \prime}(r) + \frac{2}{r} \rho^{\prime}(r) - m_{\rho}^2 \rho &=& \frac{g_{\rho}}{2} \left[n_n(r)-n_p(r)\right] + 2 g_{\rho}^2 \rho f + \frac{\xi}{6} g_{\rho}^4 \rho^3 \end{eqnarray*}\end{split}\]

and the Coulomb field equation is

\[ A^{\prime \prime}(r) + \frac{2}{r} A^{\prime}(r) = - e n_p(r) \]

The meson field equations plus a pair of Dirac-like nucleon field equations for each index \( \alpha \) must be solved to compute the structure of a given nucleus.

If the delta meson is included, there is another field equation

\[ \delta^{\prime \prime}(r) + \frac{2}{r} \delta^{\prime}(r) - m_{\delta}^2 \delta = - \frac{g_{\delta}}{2} \left[n_n(r)-n_p(r)\right] + 2 g_{\delta}^2 \delta f + \frac{\xi}{6} g_{\delta}^4 \delta^3 \]

The densities (scalar, baryon, isospin, and charge) are

\[\begin{split}\begin{eqnarray*} n_s(r) &=& \sum_{\alpha} \left\{ \int d^3 r \left[ G(r)^2-F(r)^2 \right] \right\} \\ n(r) &=& \sum_{\alpha} \left\{ \int d^3 r \left[ G(r)^2+F(r)^2 \right] \right\} \\ n_i(r) &=& \sum_{\alpha} \left\{ t_{\alpha} \int d^3 r \left[ G(r)^2-F(r)^2 \right] \right\} \\ n_c(r) &=& \sum_{\alpha} \left\{ \left[t_{\alpha}+\frac{1}{2}\right] \int d^3 r \left[ G(r)^2-F(r)^2 \right] \right\} \end{eqnarray*}\end{split}\]

The total energy is

\[ E = \sum_{\alpha} \varepsilon_{\alpha} (2 j_{\alpha}+1) - \frac{1}{2} \int d^{3} r \left[- g_{\sigma} \sigma(r) \rho_s(r) + g_{\omega} \omega(r) \rho(r) + \frac{1}{2} g_{\rho} \rho(r) + e A(r) n_p(r) \right] \]

The charge density is the proton density folded with the charge density of the proton, i.e.

\[ \rho_{\mathrm{ch}}(r) = \int d^{3} r^{\prime} \rho_{\mathrm{prot}}(r-r^{\prime}) \rho_p(r) \]

where the proton charge density is assumed to be of the form

\[ \rho_{\mathrm{prot}}(r) = \frac{\mu^3}{8 \pi} \exp \left( - \mu |r|\right) \]

and the parameter \(\mu = (0.71)^{1/2}~\mathrm{GeV}\) (see Eq. 20b in [Horowitz81]). The default value of nucleus_rmf::a_proton is the value of \(\mu\) converted into \(\mathrm{fm}^{-1}\).

Generally, the first array index associated with a function is not the value at \( r=0 \), but at \( r=\Delta r \) (stored in step_size) and so the \( r=0 \) part of the algorithm is handled separately.

Todo:

Better documentation

Convert energies() to use EOS and possibly replace sigma_rhs() and related functions by the associated field equation method of eos_had_rmf.

Todo:: I believe currently the center of mass correction for the binding energy per nucleon is not done and has to be added in after the fact

Idea for Future:

Sort energy levels at the end by energy

Improve the numerical methods

Make the neutron and proton orbitals more configurable

Generalize to \( m_n \neq m_p \) .

Allow more freedom in the integrations

Consider converting everything to inverse fermis (AWS 6/12/22: I think this has been at least partially completed)

Convert to zero-indexed arrays (mostly done)

Warn when the level ordering is wrong, and unoccupied levels are lower energy than occupied levels

Connect with o2scl::nucmass ?

Note

This class does not currently have the ability to assign fractional occupation numbers to the various shells.

Note

In this code, the meson fields are internally stored with together with their nucleon couplings, so fields contains \( (g_{\sigma} \sigma, g_{\omega} \omega, g_{\rho} \rho, A ) \) where \( g_{\sigma} \) is given, e.g. by the product of eos_had_rmf::ms times eos_had_rmf::cs .

Numeric configuration

static const int grid_size = 300¶: The grid size.

double a_proton¶

The parameter for the charge density of the proton in \( 1/\mathrm{fm} \) (default is about 4.27073)

This quantity is used in pfold() and gauss()

eos_had_rmf *rmf¶: The base EOS.

std::shared_ptr<table_units<>> profiles¶: The radial profiles.

std::shared_ptr<table_units<>> chden_table¶: The final charge densities.

std::vector<shell> *levp¶: A pointer to the current vector of levels (either levels or unocc_levels)

int verbose¶: Control output (default 1)

shell neutron_shells[n_internal_levels]¶: The starting neutron levels.

shell proton_shells[n_internal_levels]¶: The starting proton levels.

double step_size¶: The grid step size (default 0.04)

double mnuc¶: The nucleon mass (automatically set in init_fun())

bool err_nonconv¶: If true, call the error handler if the routine does not converge or reach the desired tolerance (default true)

int itmax¶: Maximum number of total iterations (default 70)

int meson_itmax¶: Maximum number of iterations for solving the meson field equations (default 10000)

int dirac_itmax¶: Maximum number of iterations for solving the Dirac equations (default 100)

double dirac_tol¶: Tolerance for Dirac equations (default \( 5 \times 10^{-3} \) ).

double dirac_tol2¶: Second tolerance for Dirac equations (default \( 5 \times 10^{-4} \) ).

double meson_tol¶: Tolerance for meson field equations (default \( 10^{-6} \) ).

initial_guess ig¶

Parameters for initial guess.

Default is \( (310/(\hbar c),240/(\hbar c), -6/(\hbar c),25.9/(\hbar c), 6.85, 0.6) \) where the first four are in \( 1/\mathrm{fm} \) and the last two are in fm.

bool generic_ode¶: If true, use the generic ODE solver instead of the internal 4th order Runge-Kutta (default false)

inline void set_step(ode_step<ubvector, ubvector, ubvector, ode_funct, double> &step)¶: Set the stepper for the Dirac differential equation.

int load_nl3(eos_had_rmf &r)¶: Load the default model NL3 into the given eos_had_rmf object.

The meson and photon fields and field equations (protected)

ubmatrix field0¶: Values of the fields from the last iteration.

ubmatrix fields¶: The meson fields (times the nucleon couplings) and photon fields.

ubmatrix gin¶: The Green’s functions inside.

ubmatrix gout¶: The Green’s functions outside.

int surf_index¶: The grid index corresponding to the nuclear surface (computed by init_run())

double sigma_rhs(double sig, double ome, double rho)¶: Scalar density RHS.

double omega_rhs(double sig, double ome, double rho)¶: Vector density RHS.

double rho_rhs(double sig, double ome, double rho)¶: Isovector density RHS.

void meson_init()¶: Calculate meson and photon Green’s functions gin and gout.

void meson_iter(double ic)¶

Calculate meson and photon fields.

The meson and photon field equations are of the Sturm-Liouville form, e.g.

\[ \left[r^2 \sigma^{\prime}(r) \right]^{\prime} - r^2 m_{\sigma}^2 \sigma(r) = f(r) \]

where \( \sigma(0) = \sigma_0 \) and \( \sigma(+\infty) = 0 \). A solution of the homogeneous equation with \( f(r) =0 \) is \( \sigma(r) = \sigma_0 \sinh( m_{\sigma} r)/ (m_{\sigma} r) \). The associated Green’s function is

\[ D(r,r^{\prime},m_{\sigma}) = \frac{-1}{m_{\sigma} r r^{\prime}} \sinh (m_{\sigma} r_{<}) \exp (-m_{\sigma} r_{>}) \]

where \( r_{<} \) is the smaller of \( r \) and \( r^{\prime} \) and \( r_{>} \) is the larger. One can then write the solution of the meson field equation given the density

\[ \sigma(r) = \int_0^{\infty} r^{\prime 2}~d r^{\prime} \left[ - g_{\sigma} n_{s}(r) \right] D\left(r,r^{\prime},m_{\sigma}\right) \]

Since radii are positive, \( \sinh (m r) \approx e^{m r}/2 \) and

\[ D(r,r^{\prime},m_{\sigma}) \approx \left[ \frac{-1}{2 m_{\sigma} r_{<}} \exp (m_{\sigma} r_{<}) \right] \left[ \frac{1}{r_{>}} \exp (-m_{\sigma} r_{>}) \right] \]

The two terms in the Green’s function are separated into

\[ \mathrm{gin}(r) = \frac{e^{m_{\sigma} r}}{2 m_{\sigma} r} \]

and

\[ \mathrm{gout}(r) = \frac{e^{-m_{\sigma} r}}{r} \, . \]

These functions are computed in meson_init() . Then the field is given by

\[ \sigma(r)= \left(\frac{e^{-m_{\sigma} r}}{r}\right) \int_0^{r} r^{\prime 2} g_{\sigma} n_{s} \left(\frac{e^{m_{\sigma} r^{\prime}}}{2 m_{\sigma} r^{\prime}} \right)~d r^{\prime} + \left(\frac{e^{m_{\sigma} r}}{2 m_{\sigma} r} \right) \int_r^{\infty} r^{\prime 2} g_{\sigma} n_{s} \left(\frac{e^{-m_{\sigma} r^{\prime}}} {r^{\prime}}\right)~d r^{\prime} \]

where the first integral is stored in xi2 and the second is in xi1 in the function meson_iter() . The part of xi2 at \( r^{\prime}=0 \) is stored in xi20.

int meson_solve()¶: Solve for the meson profiles.

Density information (protected)

ubmatrix xrho¶: The densities times radius squared.

ubvector xrhosp¶: The proton scalar density times radius squared.

ubvector xrhos¶: The scalar field RHS.

ubvector xrhov¶: The vector field RHS.

ubvector xrhor¶: The isubvector field RHS.

ubvector chden1¶: Charge density.

ubvector chdenc¶: Charge density.

ubvector arho¶: Baryon density.

ubvector energy¶: Energy integrand.

void init_meson_density()¶: Initialize the meson and photon fields, the densities, etc.

int energy_radii(double xpro, double xnu, double e)¶: Calculate the energy profile.

void center_mass_corr(double atot)¶: Compute the center of mass correction.

Calculating the form factor, etc. (protected)

interp_vec<ubvector> *gi¶: Interpolation object.

void pfold(double x, double &xrhof)¶: Fold in proton form factor.

double xpform(double x, double xp, double a)¶: Function representing proton form factor.

void gauss(double xmin, double xmax, double x, double &xi)¶: Perform integrations for form factor.

double xrhop(double x1)¶: Desc.

Solving the Dirac equations (protected)

ode_rkck_gsl<ubvector, ubvector, ubvector, ode_funct, double> def_step¶: The default stepper.

ode_step<ubvector, ubvector, ubvector, ode_funct, double> *ostep¶: The ODE stepper.

bool init_called¶: True if init() has been called.

ubvector ode_y¶: ODE functions.

ubvector ode_dydx¶: ODE derivatives.

ubvector ode_yerr¶: ODE errors.

int dirac(int ilevel)¶

Solve the Dirac equations.

Solves the Dirac equation in from 12 fm to the match point and then out from .04 fm and adjusts eigenvalue with

\[ \Delta \varepsilon = -g(r=\mathrm{match\_point}) \times (f^{+}-f^{-}) \]

void dirac_step(double &x, double h, double eigen, double kappa, ubvector &varr)¶: Take a step in the Dirac equations.

int odefun(double x, size_t nv, const ubvector &y, ubvector &dydx, odparms &op)¶: The form of the Dirac equations for the ODE stepper.

void field(double x, double &s, double &v, ubvector &varr)¶: Compute the fields for the Dirac equations.

double dirac_rk4(double x, double g1, double f1, double &funt, double eigen, double kappa, ubvector &varr)¶: Integrate the Dirac equations using a simple inline 4th order Runge-Kutta.

Gauss-Legendre integration points and weights

double x12[6]¶

double w12[6]¶

double x100[50]¶

double w100[50]¶

Results

int nlevels¶: The number of levels.

std::vector<shell> levels¶

The levels (protons first, then neutrons)

An array of size nlevels

int nuolevels¶: The number of unoccupied levels (equal to unocc_Z + unocc_N)

std::vector<shell> unocc_levels¶

The unoccupied levels (protons first, then neutrons)

An array of size nuolevels

double stens¶

Surface tension (in \( \mathrm{fm}^{-3} \) )

Computed in post_converge() or automatically in run_nucleus()

double rnrp¶

Skin thickness (in fm)

Computed every iteration in iterate() or automatically in run_nucleus()

double rnrms¶

Neutron RMS radius (in fm)

Computed every iteration in iterate() or automatically in run_nucleus()

double rprms¶

Proton RMS radius (in fm)

Computed every iteration in iterate() or automatically in run_nucleus()

double etot¶

Total energy (in MeV)

Computed every iteration in iterate() or automatically in run_nucleus()

double r_charge¶

Charge radius (in fm)

Computed in post_converge() or automatically in run_nucleus()

double r_charge_cm¶

Charge radius corrected by the center of mass (in fm)

Computed in post_converge() or automatically in run_nucleus()

inline std::shared_ptr<table_units<>> get_profiles()¶

Get the radial profiles.

The profiles are calculated each iteration by iterate().

inline std::shared_ptr<table_units<>> get_chden()¶: The final charge densities.

Equation of state

eos_had_rmf def_rmf¶

The default equation of state (default NL3)

This is set in the constructor to be the default model, NL3, using the function load_nl3().

thermo hb¶

thermo object for the EOS

This is just used as temporary storage.

fermion n¶

The neutron.

The mass of the neutron is ignored and set by init_run() to be eos_had_rmf::mnuc from rmf.

fermion p¶

The proton.

The mass of the proton is ignored and set by init_run() to be eos_had_rmf::mnuc from rmf.

inline int set_eos(eos_had_rmf &r)¶

Set the base EOS to be used.

The equation of state must be set before run_nucleus() or init_fun() are called, including the value of eos_had_rmf::mnuc.

Basic operation

int run_nucleus(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)¶

Computes the structure of a nucleus with the specified number of levels.

Note that rmf must be set before run_nucleus() is called.

This calls init_run(), and then iterate() until iconverged is 1, and then post_converge().

inline void set_verbose(int v)¶: Set output level.

Lower-level interface

void init_run(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)¶

Initialize a run.

Note that rmf must be set before init_run() is called.

int iterate(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N, int &iconverged, int &dirac_converged, int &meson_converged)¶: Perform an iteration.

int post_converge(int nucleus_Z, int nucleus_N, int unocc_Z, int unocc_N)¶: After convergence, make CM corrections, etc.

Public Functions

nucleus_rmf()¶

~nucleus_rmf()¶

Protected Attributes

bool include_delta¶: If true, include the scalar isovector meson.

Protected Static Attributes

static const int n_internal_levels = 29¶: The total number of shells stored internally.

Private Types

typedef boost::numeric::ublas::vector<double> ubvector¶

typedef boost::numeric::ublas::matrix<double> ubmatrix¶

struct initial_guess¶

Initial guess structure.

The initial guess for the meson field profiles is a set of Fermi-Dirac functions, i.e.

\[ \sigma(r)=\mathrm{sigma0}/ [1+\exp((r-\mathrm{fermi\_radius})/\mathrm{fermi\_width})] \]

Public Members

double sigma0¶: Scalar field times the nucleon coupling at \( r=0 \) in \( 1/\mathrm{fm} \).

double omega0¶: Vector field times the nucleon coupling at \( r=0 \) in \( 1/\mathrm{fm} \).

double rho0¶: Isovector field times the nucleon coupling at \( r=0 \) in \( 1/\mathrm{fm} \).

double A0¶: Coulomb field at \( r=0 \) in \( 1/\mathrm{fm} \).

double fermi_radius¶: The radius (in fm) for which the fields are half their central value

double fermi_width¶: The “width” of the Fermi-Dirac function (in fm)

struct odparms¶

A convenient struct for the solution of the Dirac equations.

Public Members

double eigen¶: Eigenvalue.

double kappa¶: Quantum number \( \kappa \) .

ubmatrix *fields¶: The meson fields times the nucleon couplings.

ubvector *varr¶: Desc.

struct shell¶

A shell of nucleons for nucleus_rmf.

Public Members

int twojp1¶: Degeneracy \( 2 j+1 \) .

int kappa¶: \( \kappa \)

double energy¶: Energy eigenvalue.

double isospin¶: Isospin ( \( +1/2 \) or \( -1/2 \) ).

std::string state¶: Angular momentum-spin state \( ^{2s+1} \ell_{j} \).

double match_point¶: Matching radius (in fm)

double eigen¶: Energy eigenvalue.

double eigenc¶: Correction to the energy eigenvalue.

int nodes¶: Number of nodes in the wave function.