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 Prakash94. 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 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
- Todo:
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
- Idea for Future:
Improve the numerical methods
- Idea for Future:
Make the neutron and proton orbitals more configurable
- Idea for Future:
Generalize to \( m_n \neq m_p \) .
- Idea for Future:
Allow more freedom in the integrations
- Idea for Future:
Consider converting everything to inverse fermis.
- Idea for Future:
Convert to zero-indexed arrays (mostly done)
- Idea for Future:
Warn when the level ordering is wrong, and unoccupied levels are lower energy than occupied levels
- Idea for Future:
Connect with o2scl::nucmass ?
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().
-
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.
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()
-
std::shared_ptr<table_units<>>
get_profiles
()¶ Get the radial profiles.
The profiles are calculated each iteration by iterate().
-
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.
-
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.
Numeric configuration
-
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)
-
double
step_size
¶ The grid step size (default 0.04)
-
double
mnuc
¶ The nucleon mass (automatically set in init_fun())
-
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.
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]¶
Protected Static Attributes
-
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
-
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