Class Documentation¶

seminf_nr¶

class seminf_nr¶

Semi-infinite nuclear matter for nonrelativistic models in the Thomas-Fermi approximation.

This discussion is a combination of Jim Lattimer’s notes and the material in Steiner05ia .

Given a energy density functional separated into bulk and gradient terms

\[ \varepsilon(n_n,n_p,\nabla n_n,\nabla n_p ) = \varepsilon_{\mathrm{bulk}}(n_n,n_p) + \frac{Q_{nn}}{2} \left( \nabla n_n \right)^2 + \frac{Q_{pp}}{2} \left( \nabla n_p \right)^2 + Q_{np} \nabla n_n \nabla n_p \]

one can minimize the thermodynamic potential

\[ \Omega = \int \left(-P\right) dV = \int~dV \left[ f_{\mathrm{bulk}} - \mu_n n_n - \mu_p n_p + + \frac{Q_{nn}}{2} \left( \nabla n_n \right)^2 + \frac{Q_{pp}}{2} \left( \nabla n_p \right)^2 + Q_{np} \nabla n_n \nabla n_p \right] \]

to obtain the neutron and proton density profiles. The Euler-Lagrange equations are

\[ Q_{nn} \nabla^2 n_n + Q_{np} \nabla^2 n_p - \frac{1}{2} \left[ \frac{\partial Q_{nn}}{\partial n_n} \left(\nabla n_n\right)^2 + \frac{\partial Q_{np}}{\partial n_n} \nabla n_n \nabla n_p + \frac{\partial Q_{pp}}{\partial n_n} \left(\nabla n_p\right)^2 + \right] = \frac{f_{\mathrm{bulk}}}{n_n} - \mu_{n,0} \]

and a second equation with \( n \leftrightarrow p \), where \( f_{\mathrm{bulk}} = \varepsilon_{\mathrm{bulk}} - T s_{\mathrm{bulk}} \) and \( \mu_{n,0} \) is a constant. Multiplying the first equation by \( \nabla n_n \) and the second by \( \nabla n_p \), adding, and integrating gives

\[ \frac{1}{2} \left[ Q_{nn} (\nabla n_n)^2 + Q_{pp} (\nabla n_p)^2 + 2 Q_{np} \nabla n_n \nabla n_p \right] = f_{\mathrm{bulk}}(n_B) - f_{\mathrm{bulk}}(n_{B,drip}) - \mu_{n,0} ( n_n - n_{n,\mathrm{drip}}) - \mu_{p,0} ( n_p - n_{p,\mathrm{drip}}) \]

Then the surface thermodynamic potential is

\[ \Omega + P_0 V = 2 \int \left[ f_{\mathrm{bulk}}(n_B) - f_{\mathrm{bulk}}(n_{B,\mathrm{drip}}) - \mu_{n,0} ( n_n - n_{n,\mathrm{drip}}) - \mu_{p,0} ( n_p - n_{p,\mathrm{drip}}) \right]~dV \]

In the semi-infinite matter approximation with a single coordinate, \( x \), the E-L equations become

\[\begin{split}\begin{eqnarray} Q_{nn} n_n^{\prime \prime} + Q_{np} n_p^{\prime \prime} = \mu_n - \mu_{n,0} \nonumber \\ Q_{np} n_n^{\prime \prime} + Q_{pp} n_p^{\prime \prime} = \mu_p - \mu_{p,0} \end{eqnarray}\end{split}\]

The surface tension is

\[\begin{split}\begin{eqnarray*} \omega &=& \int_{-\infty}^{\infty} \left( f - \mu_{n,0} n_n - \mu_{p,0} n_p \right) dx \nonumber \\ \omega &=& 2 \int_{-\infty}^{\infty} \left( f_{\mathrm{bulk}} - \mu_{n,0} n_n - \mu_{p,0} n_p \right) dx \end{eqnarray*}\end{split}\]

with the second equality because the bulk and gradient contributions to the surface tension are exactly the same.

In semi-infinite matter one can define a radius by

\[ n_{i,LHS}(R_i - L) = \int_{-L}^{\infty} n_i~dz \, , \]

and using this defintion, the neutron skin thickness in semi-infinite matter is

\[ t \equiv \int_{-\infty}^{\infty} \left( \frac{n_n(z)} {n_{n,\mathrm{LHS}}} - \frac{n_p(z)}{n_{p,\mathrm{LHS}}} \right)~dz \]

and converting to a spherical geometry gives a factor of \( \sqrt{3/5} \) so the neutron skin thickness is \( R_n - R_p = t \sqrt{3/5} \) .

Derive:

\[ \omega(\delta=0) = \sqrt{Q_{nn} + Q_{np}} \int_{-\infty}^{\infty} \left[ f_{\mathrm{bulk}}(n_B,\alpha=0)+n_B B \right]^{1/2} dn_B \]

In neutron-rich matter, the proton density vanishes while the neutron density extends out to \( x \rightarrow \infty \). In order to handle this, the algorithm separates the solution into two intervals: the first where the proton density is non-zero, and the second where the proton density is zero. These two solutions are then pasted together to create the final table.

12/15/03 - Implemented a new solution method for Qnn!=Qnp for when nn_drip becomes > 0. It seems that in this case, it is better not to automatically adjust the scale of the x-axis, but use a gradual broadening similar to what we have done using the relativistic code. This allows calculation down to really low proton fractions! We need to see if this is consistent with the analytic code for both Qnn<Qnp (hard) and Qnn>Qnp (probably works).

Note: The first proton fraction cannot be 0.5, as the algorithm cannot handle later values different from 0.5

Storage for the solution

si_vector_t xstor¶

si_matrix_t ystor¶

Store initial guesses for next proton fraction

si_vector_t xg1¶

si_vector_t xg2¶

si_matrix_t yg1¶

si_matrix_t yg2¶

Nucleon-nucleon interactions

eos_had_apr eosa¶

eos_had_potential eosg¶

eos_had_skyrme eoss¶

eos_had_schematic eosp¶

Public Functions

seminf_nr()¶

int calc(std::vector<std::string> &sv, bool itive_com)¶: Calculate nucleon profiles for a list of proton fractions.

int check(std::vector<std::string> &sv, bool itive_com)¶: Compute standard results and compare with stored output.

Public Members

double rhs_adjust¶: Adjustment for RHS boundary.

string model¶: Type of EOS model in use.

string out_file¶: Name of output file.

Protected Functions

int derivs(double sx, const si_vector_t &sy, si_vector_t &dydx)¶: The differential equations to be solved.

double solve_qnn_neq_qnp(double lmonfact)¶: Solve the ODEs when \( Q_{nn} = Q_{np} \).

int qnnqnpfun(size_t sn, const si_vector_t &sx, si_vector_t &sy)¶: Desc.

int ndripfun(size_t sn, const si_vector_t &sx, si_vector_t &sy)¶: Compute the boundary densities in the case of a neutron drip.

int pdripfun(size_t sn, const si_vector_t &sx, si_vector_t &sy)¶: Compute the boundary densities in the case of a proton drip.

double difeq(size_t ieq, double x, si_matrix_row_t &y)¶: Reformat differential equations for o2scl::ode_it_solve.

double left(size_t ieq, double x, si_matrix_row_t &y)¶: LHS boundary conditions for o2scl::ode_it_solve.

double right(size_t ieq, double x, si_matrix_row_t &y)¶: RHS boundary conditions for o2scl::ode_it_solve.

Protected Attributes

double mun¶: Neutron chemical potential.

double mup¶: Proton chemical potential.

double nn_left¶: Neutron number density on the LHS.

double np_left¶: Proton number density on the LHS.

double qnn¶: Isoscalar gradient term for neutrons.

double qpp¶: Isoscalar gradient term for protons.

double qnp¶: Isovector gradient term.

double n0half¶: Saturation density of isospin-symmetric matter.

bool show_relax¶: Desc.

double barn¶: Desc.

double dmundn¶: Desc.

double dmundp¶: Desc.

double nsat¶: Saturation density of matter with current proton fraction.

double protfrac¶: Current proton fraction (in high-density region)

int pf_index¶: Proton fraction index.

eos_had_base *eos¶: EOS pointer.

fermion neutron¶: Neutron.

fermion proton¶: Proton.

thermo hb¶: Thermodynamic functions.

bool relax_file¶: If true, output relaxation iterations (default false)

bool qnn_equals_qnp¶: If true, \( Q_{nn}=Q_{np} \).

bool low_dens_part¶: If true, match to exponential decay on the RHS.

double min_err¶: Minimum err.

double monfact¶: Desc.

double expo¶: Desc.

double nn_drip¶: Neutron drip density.

double np_drip¶: Proton drip density.

bool flatden¶: Desc (default false)

mroot_hybrids<mm_funct, si_vector_t, si_matrix_t, jac_funct> nd¶: Nonlinear equation solver.

double nnrhs¶: Neutron drip density.

double nprhs¶: Proton drip density.

double rhs_length¶: Size of the low-density surface region.

bool goodbound¶: If true, we have good boundary conditions.

double first_deriv¶: Initial derivative for densities on LHS (typically negative)

double big¶: Desc (default 3.0)

int ngrid¶: The grid size (default 100)

double relax_converge¶: Desc.

double rhs_min¶: Desc.

double initialstep¶: The step size factor for constructing the initial guess (default 7.0)

inte_qagiu_gsl gl2¶: Integrator to compute integrals of final solution.

table_units tab_prof¶: Store the profiles.

table_units tab_summ¶: Store results across several profiles.

deriv_gsl df¶: To compute derivatives (currently unused)

seminf_rel¶

class seminf_rel¶

Semi-infinite nuclear matter for relativistic mean-field models in the Thomas-Fermi approximation.

RMF models in the semi-infinite matter approximation were first computed in Boguta77 and Stocker91 . This code was developed for Steiner05ia.

In the semi-infinite nuclear matter approximation, the meson field equations are

\[\begin{split}\begin{eqnarray*} \frac{d^2 \sigma}{d z^2} &=& m_{\sigma}^2 \sigma - g_{\sigma} \left( n_{s n} + n_{s p} \right) + b M g_{\sigma}^3 \sigma^2 + c g_{\sigma}^4 \sigma^3 - g_{\rho}^2 \rho^2 \frac{\partial f}{\partial \sigma} \nonumber \\ \frac{d^2 \omega}{d z^2} &=& m_{\omega}^2 \omega - g_{\omega} \left(n_n+n_p\right) + \frac{\zeta}{6} g_{\omega}^4 \omega^3 + g_{\rho}^2 \rho^2 \frac{\partial f}{\partial \omega} \nonumber \\ \frac{d^2 \rho}{d z^2} &=& m_{\rho}^2 \rho + \frac{1}{2} g_{\rho} \left(n_n-n_p\right) + 2 g_{\rho}^2 \rho f + \frac{\xi}{6} g_{\rho}^4 \rho^3 \end{eqnarray*}\end{split}\]

in the same notation as o2scl::eos_had_rmf .

The algorithm works by solving first over a small interval in coordinate space and slowly spreading the solution out until the derivatives of the meson fields at the left boundary drop below a specified tolerance.

Note: wd2int() is broken since fesym() tends to fail at low densities. For this reason, it wasn’t used in Steiner05ia .

The meson fields at the LHS

double sigma_left¶

double omega_left¶

double rho_left¶