Class nucleus_rmf (o2scl)¶
-
class o2scl::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 in eos_had_rmf for [Steiner05b], and then converted to C++ by Steiner afterwards.
The standard usage is something like:
which computes the structure of \( ^{208}\mathrm{Pb} \) and outputs the neutron skin thickness using the modelnucleus_rmf rn; o2scl_hdf::rmf_load(rn.rmf,"NL4"); rn.run_nucleus(82,208,0,0); cout << rn.rnrp << endl;
'NL4'
.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
1 S 1/2 // 2 nucleons 1 P 3/2 1 P 1/2 // 8 nucleus 1 D 5/2 1 D 3/2 2 S 1/2 // 20 nucleons 1 F 7/2 // 28 nucleons 1 F 5/2 2 P 3/2 2 P 1/2 // 40 nucleons 1 G 9/2 // 50 nucleus 1 G 7/2 2 D 5/2 1 H 11/2 2 D 3/2 3 S 1/2 // 82 nucleons 1 H 9/2 2 F 7/2 1 I 13/2 2 F 5/2 3 P 3/2 3 P 1/2 // 126 nucleons 2 G 9/2 1 I 11/2 1 J 15/2 3 D 5/2 4 S 1/2 2 G 7/2 3 D 3/2 // 184 nucleons
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.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 ref 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.
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 ?
Numeric configuration
-
static const int grid_size = 300¶
The grid size.
-
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,240,-6,25.9,6.85,0.6}
-
bool generic_ode¶
If true, use the generic ODE solver instead of the internal 4th order Runge-Kutta.
-
double a_proton¶
The parameter for the charge density of the proton (default is about 4.27073)
-
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())
-
inline void set_step(ode_step<ubvector, ubvector, ubvector, ode_funct> &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.
Results
-
int nlevels¶
The number of levels.
-
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.
The meson and photon fields and field equations (protected)
-
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)¶
Isubvector density RHS.
-
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 inxi2
and the second is inxi1
in the function meson_iter() . The part ofxi2
at \( r^{\prime}=0 \) is stored inxi20
.
-
int meson_solve()¶
Solve for the meson profiles.
Density information (protected)
-
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)
-
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)
-
bool init_called¶
True if init() has been called.
-
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.
Gauss-Legendre integration points and weights
-
double x12[6]¶
-
double w12[6]¶
-
double x100[50]¶
-
double w100[50]¶
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 run_nucleus() 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.
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 at r=0.
-
double omega0¶
Vector field at r=0.
-
double rho0¶
Isubvector field at r=0.
-
double A0¶
Coulomb field at r=0.
-
double fermi_radius¶
The radius for which the fields are half their central value.
-
double fermi_width¶
The “width” of the Fermi-Dirac function.
-
double sigma0¶
-
struct odparms¶
A convenient struct for the solution of the Dirac equations.
-
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¶
Desc.
-
double eigenc¶
Desc.
-
int nodes¶
Number of nodes in the wave function.
-
int twojp1¶