Class Documentation

sna_thermo

class crust::sna_thermo

Thermodynamics in the single-nucleus approximation.

Different parts of the free energy

double part1

double part2

double part3

double part4

Public Functions

sna_thermo(ldrop_crust &lc, o2scl::eos_had_temp_base &he)

Create with the specified mass model and EOS.

int free_energy_sna(matter &m, double T)

Compute the free energy density in fm^{-4} in the single nucleus approximation.

This is faster than free_energy_dist() for one species since no iteration is required.

int free_energy_sna_fix_nb_nn(double nb, matter &m, double T)

For fixed nb and nn, determine nnuc and compute the free energy.

Used in crust_driver::compute_sna().

int free_energy_sna_fix_nb_nnuc(double nb, matter &m, double T)

For fixed nb and nnuc, determine nn and compute the free energy.

Used by sna_thermo::free_press_sna_fixed_ye .

int check_free_energy_sna(dist_thermo &dt, o2scl::test_mgr &t)

Check free_energy_sna() by computing the free energy a different way.

Public Members

double mag_field

The magnetic field in Gauss.

o2scl::eos_had_temp_base *het

The hadronic eos.

This needs to be public so it can be changed by crust_driver::model() .

Protected Functions

double get_phi_sna(matter &m, double T)

Compute the total fractional volume occupied by the the neutrons in all nuclei.

This function computes

\[ \phi \equiv \sum_i \phi_i \equiv \sum_i \frac{4}{3} \pi R_{n,i}^3 n_i \]

Used in rxns and in crust(), but can probably be called directly from gibbs_fixp_neutron() and/or gibbs_fixp().

int baryon_density_sna(matter &m, double T)

Compute total baryon number density (with dripped neutrons)

Protected Attributes

bool inc_nuc_trans

o2scl::fermion_rel relf

Electron thermodynamics.

o2scl::classical cla

Classical thermodynamics.

o2scl::fermion_mag_zerot mfz

Electron in a magnetic field.

o2scl::fermion elec_B

Electron.

o2scl::convert_units &cng

Convert units.

ldrop_crust *lda

The mass model.

class free_energy_sna_fix_nb_nnhat

To test minimization with respect to \( {\hat n}_n \) instead of \( n_n \).

Used at the end of crust::crust_driver::compute_sna().

Inherits from crust::sna_thermo::snat_funct_base

Public Types

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

Public Functions

free_energy_sna_fix_nb_nnhat(sna_thermo &st, double nb, matter &m, double &T)

Specify the member function pointer.

virtual ~free_energy_sna_fix_nb_nnhat()

virtual double operator()(double nnhat)

Compute the function at point x, with result y.

Protected Attributes

o2scl::mroot_hybrids gmh

The solver.

ubvector x

Hold the solution for the neutron and nuclear densities.

solve_nn_ni_fixed_nb_nnhat mfna

The object to solve.

double nb_

Baryon density.

class free_energy_sna_neut

Class to minimize the single-nucleus free energy at fixed baryon density as a function of the neutron number density.

Used in crust_driver::compute_sna().

Inherits from crust::sna_thermo::snat_funct_base

Public Functions

free_energy_sna_neut(sna_thermo &st, matter &m, double T, double nb)

double operator()(double nn) const

Free energy as a function of neutron drip density.

Public Members

double nb_

Baryon number density.

class free_press_sna_fixed_ye

Class to give the properties of matter as a function of baryon density at fixed electron fraction in the single nucleus approximation.

Inherits from crust::sna_thermo::snat_funct_base

Public Functions

free_press_sna_fixed_ye(sna_thermo &st, matter &m, double T, dist_thermo &dt, bool return_fr)

double operator()(double nb) const

Free energy as a function of neutron drip density.

Public Members

dist_thermo *dtp

Dist thermo object for gibbs energy.

double Ye_

Electron fraction.

bool return_fr_

If true, return free energy. Otherwise, return Gibbs energy.

class snat_funct_base

A base class for function objects.

Subclassed by crust::sna_thermo::free_energy_sna_fix_nb_nnhat, crust::sna_thermo::free_energy_sna_neut, crust::sna_thermo::free_press_sna_fixed_ye, crust::sna_thermo::solve_nn_ni_fixed_nb_nnhat

Public Functions

snat_funct_base(sna_thermo &st, matter &m, double &T)

Create with specified.

Public Members

sna_thermo *stp

Object of type sna_thermo.

matter *mp

Matter object.

double T_

Temperature.

class solve_nn_ni_fixed_nb_nnhat

Return the free energy as a one-dimensional function of the neutron number density.

Used in crust::crust_driver::compute_sna2().

Inherits from crust::sna_thermo::snat_funct_base

The parameters

double nb_

double nnhat_

Public Types

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

Public Functions

solve_nn_ni_fixed_nb_nnhat(sna_thermo &st, matter &m, double &T)

void set(sna_thermo &st, double nb, double nnhat, matter &m, double &T)

Specify the member function pointer.

virtual int operator()(size_t nv, const ubvector &x, ubvector &y)

Compute the function at point x, with result y.

dist_thermo

class crust::dist_thermo

An object to compute the thermodynamic properties of a distribution.

The variable \( \chi \) is the fractional volume occupied by the neutrons relative to the rest of the Wigner-Seitz cell

\[ \chi_i \equiv \left( \frac{R_{n,i}}{R_{\mathrm{WS},i}} \right)^3 \]

Different parts of the free energy

double part1

double part2

double part3

double part4

Testing

int check

const int check_none =0

const int check_mass_density =1

const int check_free_energy_dist =3

const int check_ldrop_derivs =5

const int check_free_energy_cell =8

Public Types

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

Public Functions

dist_thermo(o2scl::convert_units &conv)

Create with the specified unit conversion class.

int check_pressure(matter &m, double T)

Check the pressure computed in gibbs_energy_dist()

int free_energy_dist_sna(double nb, matter &m, double T)

Given nb and nn, determine nnuc and compute the free energy.

int free_energy_dist_alpha(matter &m, double &T, double alpha)

For fixed \( \alpha=n_n (1-\phi) \), solve for n_n and compute the free energy.

This function uses iteration.

int free_energy_dist(matter &m, double T, bool eval_chem_pots = false)

Compute the free energy density in \( \mathrm{fm}^{-4} \).

Also computes the baryon density and the rest mass energy density.

If eval_chem_pots is true, this function evaluates the chemical potentials and stores them in matter::mu and matter::mun.

If check is equal to check_free_energy_dist, then this function checks the free energy calculation by computing the free energy using a different method and comparing the results.

int gibbs_energy_dist(matter &m, double T)

Compute the Gibbs energy density from the free energy density.

This function calls free_energy_dist() and then computes the gibbs energy density afterwards from a sum of the product of densities and chemical potentials. It also computes the free energy, pressure, rest mass energy density, and baryon density as a by-product.

int mass_density(matter &m, double T)

Compute the rest mass density in \( \mathrm{g}/\mathrm{cm}^3 \) (including the nuclear binding energy) at temperature T.

Results are stored in matter::rho .

If check is equal to check_mass_density, then this function computes the rest mass density with an alternative approach as a check.

int free_energy_cell(matter &m, double T, int ix, double &fr, double &vol, double &Nprime)

Compute the total free energy of one cell.

Compute the free energy of cell of species ix in \( \mathrm{fm}^{-1} \), the volume in \( \mathrm{fm}^{3} \), and the value of \( N^{\prime} \equiv N - \phi n_{n,\mathrm{out}} / n_i \).

double get_phi(matter &m, double T)

Compute the total fractional volume occupied by the the neutrons in all nuclei.

This function computes

\[ \phi \equiv \sum_i \phi_i \equiv \sum_i \frac{4}{3} \pi R_{n,i}^3 n_i \]

It is used in crust::rxns, crust::pyc_rates and in crust::crust_driver::compute_sna().

int baryon_density(matter &m, double T)

Compute total baryon number density (with dripped neutrons)

int gibbs_fixp(double pr_target, matter &m_new, double T)

Create a new configuration m_new with a fixed pressure of pr_target at temperature T.

This function modifies the composition in m_new by multiplying all of the nuclear number densities and the external neutron density (in m.n.n) by the same factor until the pressure at temperature T is equal to that specified in pr_target.

Because this function operates on m.n.n and not on the quantity m.n.n*(1-phi), it is typically only used with m.n.n is zero.

Uses gibbs_pressure_func().

Currently, this function uses a local o2scl::root_cern object to solve for pressure equality, but there is some legacy code for using o2scl::root_brent_gsl to do the solution.

int gibbs_fixp_neutron(double pr_old, double nn_full_old, matter &m_new, double T, double nn_shift)

Create a new configuration m_new with pressure of pr_target at temperature T and a fixed value of \( n_n (1-\phi) \) .

Given the old neutron density and value of phi in nn_old and phi_old, this modifies the composition of m_new until its pressure at temperature T and at fixed \( \alpha = n_n (1-\phi) + n_{n,\mathrm{shift}} \) is the same as that given in pr_target.

This uses the function gibbs_press_neut_func which has two unknowns, the new external neutron number density, m_new.n.n and the correction factor corr to multiply the nuclear densities by. It solves for pressure equality, and for

\[ (1-\phi^{[2]})n^{[2]}_{n} = \mathrm{corr} (1-\phi^{[1]})n^{[1]}_{n} \]

which ensures that the reduced number of external neutrons is properly scaled by the same factor.

int gibbs_energy_per_baryon_cell(matter &m, double T, double P, int ix, double &gpb)

Compute the gibbs energy per baryon of cell with index ix at pressure P.

Uses free_energy_cell() to compute the free energy, etc.

int gibbs_energy_cell(matter &m, double T, double P, int ix, double &g_cell)

The total Gibbs energy of one cell.

Public Members

ldrop_crust *lda

The mass model.

o2scl::eos_had_temp_base *het

The hadronic eos (default points to sk)

double mag_field

The magnetic field in Gauss.

Protected Functions

double dist_sna_solve(double nb, matter &m, double T)

Function to solve for a fixed baryon density.

double gibbs_pressure_func(double pr0, double n_fact, matter &m_new, double T)

Function to solve to ensure the pressure is equal to pr0.

The variable n_fact is the factor that one should multiply all the number densities of the nuclei by in order to ensure that the new pressure is equal to pr0.

This function returns (pr-pr0)/pr0

Used in gibbs_fixp().

int gibbs_press_neut_func(double pr, double nnf, double alpha, double nn_new, matter &m_new, double T, double &y1, double &y2)

Function to solve to ensure the pressure is equal to pr and that the “full” neutron density is nnf.

Protected Attributes

o2scl::fermion_rel relf

For the electron thermodynamics.

o2scl::classical cla

For the thermodynamics of the nuclei.

o2scl::fermion_mag_zerot mfz

Thermodynamics in the magnetic field.

o2scl::fermion elec_B

Electron in a magnetic field.

o2scl::convert_units &cng

Convert units.

class dt_funct_base

A base class for function objects.

Subclassed by crust::dist_thermo::free_energy_dist_neut, crust::dist_thermo::funct_dist_sna, crust::dist_thermo::funct_solve_pressure, crust::dist_thermo::funct_solve_pressure2, crust::dist_thermo::funct_solve_pressure3

Public Functions

dt_funct_base(dist_thermo &dt, matter &m, double &T)

Create a function object.

Public Members

dist_thermo *dtp

The base pointer.

matter *mp

Matter data.

double Tv

Temperature.

class free_dist_deriv_alpha

For computing the derivative of the free energy wrt alpha in dist_thermo::check_pressure()

Public Functions

free_dist_deriv_alpha(dist_thermo &dt, matter &m, double &T, size_t pt = 0)

virtual ~free_dist_deriv_alpha()

virtual double operator()(double x)

Desc.

Public Members

dist_thermo *dtp

matter *m_

double T_

size_t pt_

class free_dist_deriv_nuc

For computing the derivative of the free energy wrt the number density of one of the nuclei in dist_thermo::check_pressure()

Public Functions

free_dist_deriv_nuc(dist_thermo &dt, matter &m, double &T, double alpha, size_t i, size_t pt = 0)

virtual ~free_dist_deriv_nuc()

virtual double operator()(double x)

Desc.

Public Members

size_t i_

dist_thermo *dtp

matter *m_

double T_

double alpha_

size_t pt_

class free_energy_dist_neut

Return the free energy as a one-dimensional function of the neutron number density.

Used in crust_driver::compute_sna_dist().

Inherits from crust::dist_thermo::dt_funct_base

Public Functions

free_energy_dist_neut(dist_thermo &dt, matter &m, double &T, double nb)

Specify the member function pointer.

virtual ~free_energy_dist_neut()

virtual double operator()(double nn) const

Compute the function at point x, with result y.

Protected Attributes

double nbx

The baryon density.

class funct_dist_sna

Return the free energy as a one-dimensional function of the nuclear number density.

Used in dist_thermo::free_energy_dist_sna() but currently commented out .

Inherits from crust::dist_thermo::dt_funct_base

Public Functions

funct_dist_sna(dist_thermo &dt, matter &m, double T, double nb)

Specify the member function pointer.

virtual ~funct_dist_sna()

virtual double operator()(double ni) const

Compute the function at point x, with result y.

Protected Attributes

double nb_

Baryon density.

class funct_solve_pressure

Solve for the pressure.

Used in dist_thermo::gibbs_fixp()

Inherits from crust::dist_thermo::dt_funct_base

Public Functions

funct_solve_pressure(dist_thermo &dt, matter &m, double &T, double pr0)

Specify the member function pointer.

virtual ~funct_solve_pressure()

virtual double operator()(double alf) const

Compute the function at point x, with result y.

Protected Attributes

double ppr0

The pressure to fix.

class funct_solve_pressure2

Solve for the pressure.

Used in dist_thermo::gibbs_fixp_neutron()

Inherits from o2scl::mm_funct< double > &) >, crust::dist_thermo::dt_funct_base

The parameters

double ppr0

double nnfull

Public Functions

funct_solve_pressure2(dist_thermo &dt, matter &m, double T, double pr0, double nn_full)

Specify the member function pointer.

virtual ~funct_solve_pressure2()

virtual int operator()(size_t nv, const ubvector &x, ubvector &y)

Compute the function at point x, with result y.

class funct_solve_pressure3

Solve for the pressure.

Used in dist_thermo::gibbs_fixp_neutron()

Inherits from o2scl::mm_funct< double > &) >, crust::dist_thermo::dt_funct_base

The parameters

double ppr0

double nnfull

Public Functions

funct_solve_pressure3(dist_thermo &dt, matter &m, double T, double pr0, double nn_full)

Specify the member function pointer.

virtual ~funct_solve_pressure3()

virtual int operator()(size_t nv, const ubvector &x, ubvector &y)

Compute the function at point x, with result y.

multi_zone

class crust::multi_zone

Multi-zone stability test.

Public Types

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

Public Functions

multi_zone()

void read_zones(size_t nz, o2scl::table_units<> &t, size_t fz)

Read table into a vector of matter objects.

int solve_fun(size_t nv, const ubvector &x, ubvector &y)

Desc.

double min_fun(size_t nv, const ubvector &x)

Desc.

void test_stability(double T, dist_thermo &dt, size_t nz, o2scl::table_units<> &t, size_t fz)

Public Members

double schwarz_km

Schwarzchild radius.

o2scl::convert_units cng

Object to convert units.

o2scl::nucmass_info nmi

To interpret element labels in accreted crust table.

ubvector nN

Number of nuclei.

ubvector nn

Number of neutrons.

int sdebug

Desc.

ubvector x_solve

Desc.

size_t n_min

Desc.

ubvector min_scale

Desc.

ubvector solve_scale

Desc.

double nn_base

Desc.

double conv

Conversion between.

ubvector nN_base

Target number of nuclei.

ubvector gm

Gravitational mass grid.

double gmtot

Gravitational mass at outer boundary.

double nntot

Total number of neutrons over all zones.

ubvector P

Pressure grid (in \( M_{\odot}/km^3 \))

ubvector r

Radial grid (in km)

ubvector eps

Energy density grid (in \( M_{\odot}/km^3 \))

size_t n_solve

Number of variables to solve for.

size_t n_zones

Number of zones.

size_t first_zone

The index of the first zone in the table.

bool first_time

If true, set the target number of nuclei in solve_fun()

std::vector<matter *> zones

Storage for the zones.

dist_thermo *dtp

To compute thermodynamics of each zone.

double Tptr

Temperature.

double Pbottom

Desc.

double Ptop

Desc.

o2scl::mroot_hybrids gmr

Solver.

o2scl::mmin_simp2 gms

Minimizer.

pyc_rates

class crust::pyc_rates

Compute pycnonuclear reaction rates.

The S-factor here is a function of \( Z_1, Z_2, A_1, A_2 \) and the energy in the center of mass, \( E \).

Ensure Z’s and A’s are even by increasing Z and decreasing A. Rates tend to decrease with Z and increase with A, so this choice minimizes fusion.

Beard10 gives computes the S-factors for several reactions and fits each to the same functional form

\[ S(E) = \exp \left\{ B_1 + B_2 E + B_3 E^2 + \frac{C_1 + C_2 E + C_3 E^2 + C_4 E^3} {1 + \exp\left[\left(E_C - E\right)/D\right]}\right\} \]

these rates are tabulated for reactions among even-even isotopes of C, O, Ne, and Mg up to a maximum value of A ( \( A_{\mathrm{max}} \)) which depends on which nuclei are being fused. S-factors for \( A > A_{\mathrm{max}} \) are assumed to be equal to those for \( A > A_{\mathrm{max}} \) (which increases fusion slightly for heavier isotopes).

Yakovlev10 gives a fit in Eqs. 5 and 6 in terms of four quantities, \( S_0, E_C, \delta \) and \( \xi \):

\[ S(E) = S_0 \exp \left[ 2 \pi \eta + \Phi(E) \right] \]

for \( E \leq E_C \) and

\[ S(E) = S_0 \exp (2 \pi \eta) \sqrt{E/E_C} \left[ 1+\xi (E-E_C)/E\right] \]

for \( E > E_C \) where \( \eta \equiv \sqrt(E_R/E) \), \( E_R \equiv \alpha^2 \mu /( 2 \hbar^2) \), \( \mu \) is the reduced mass, \( \alpha \equiv Z_1 Z_2 e^2 \), and \( \Phi(E) \equiv \) …

The fit and Eqs. 19 and 20 given in Table II of gives quantities

\[\begin{split}\begin{eqnarray*} E_C &=& \alpha/R_C^{(0)} \nonumber \\ R_C^{(0)} &=& R + \Delta R_1 |A_1 - A_{10}| + \Delta R_2 | A_2 - A_{20}| \nonumber \\ \xi &=& \xi_0 + \xi_1 (A_1 + A_2) \nonumber \end{eqnarray*}\end{split}\]

where the table gives \( R, \Delta R_{1 a}, \Delta R_{2 a}, \Delta R_{1 b}, \Delta R_{2 b}, \delta, S_0, \xi_0 \), and \( \xi_1 \) where \( A_{10} \equiv 2 Z_1 \) and \( A_{20} \equiv 2 Z_2 \) and the quantities \( \Delta R_1 \) and \( \Delta R_2 \) are given by

\[\begin{split}\begin{eqnarray*} \Delta R_1 &=& \theta(A_1-A_{10}) \Delta R_{1a} + \theta(A_{10}-A_1) \Delta R_{1b} \nonumber \\ \Delta R_2 &=& \theta(A_2-A_{20}) \Delta R_{2a} + \theta(A_{20}-A_2) \Delta R_{2b} \nonumber \end{eqnarray*}\end{split}\]

The default values of Cpyc, Cexp, and Cpl are from the “optimal” model in Table II in Yakovlev06 .

Public Functions

pyc_rates()

Initialize with “Optimal” rate.

int load_data(std::string dir)

Load the S-factor data.

This is called by the constructor to set b10.

int set_fit()

Set fit coefficients from Table in Yakovlev et al. 2010.

This is called by the constructor to set fit.

int test_Sfactor()

Test the S-factor data.

int fusion_times(dist_thermo &dt, matter &m, double T, double &t_acc, double &t_c_fusion, double &t_o_fusion, double &t_ne_fusion, double &t_mg_fusion)

Compute the fusion times for C, O, Ne, and Mg assuming one-component composition.

double Sfactor(size_t Z1, size_t Z2, size_t A1, size_t A2, double E)

Compute the S-factor in units of \( \mathrm{MeV} \cdot \mathrm{barn} \).

double Tp(size_t Z1, size_t Z2, double n12, double mu12)

Compute the plasma temperature \( T_{12}^{(p)} \) in 1/fm.

Density n12 should be in \( \mathrm{fm}^{-3} \) and reduced mass mu12 should be in 1/fm.

From Eq. 10 in Yakovlev06 this gives the plasma temperature as

\[ T_p = \sqrt{ \frac{4 \pi Z_1 Z_2 e^2 n_{12}}{2 \mu_{12}} } \]

where \( e^2 \) is replaced with gsl_num::fine_structure.

double rB(size_t Z1, size_t Z2, double m1, double m2, double mu12)

Compute \( r_{B12} \) in fm.

Reduced mass mu12 should be in 1/fm.

From Eq. 11 in Yakovlev06

\[ r_{B12} = \frac{1}{2 \mu_{12} Z_1 Z_2 e^2} \]

where \( e^2 \) is replaced with gsl_num::fine_structure and m1 and m2 are ignored.

int lambda_omega(size_t Z1, size_t Z2, double m1, double m2, double n1, double n2, double &lambda, double &omega)

Compute \( \lambda_{12} \) and \( \omega_{12} \).

Densities should be in \( \mathrm{fm}^{-3} \), masses should be in 1/fm, lambda is unitless, and omega is returned in 1/fm.

The quantity omega is computed from Tp() lambda is computed from Eq. 12 in Yakovlev06

\[ \lambda_{12} = r_{B12} \left( \frac{n12}{2} \right)^{1/3} \]

where \( n_12 \equiv \frac{3}{4 \pi a_{12}^3} \), \( a_{12} \equiv (a_1 + a_2)/2 \) , \( a_i \equiv [3/(4 \pi n_i)]^{1/3} \) and \( r_{B12} \) is given from rB().

int lambda_omega2(size_t Z1, size_t Z2, size_t A1, size_t A2, double m1, double m2, double n1, double n2, double rho, double XN, double aveZ, double aveA, double &lambda, double &omega)

Compute \( \lambda_{12} \) and \( \omega_{12} \).

Densities should be in \( \mathrm{fm}^{-3} \), masses should be in 1/fm, lambda is unitless, and omega is returned in 1/fm.

The quantity omega is computed from Tp() lambda is computed from Eq. 12 in Yakovlev06

\[ \lambda_{12} = \frac{A_1 + A_2}{A_1 A_2 Z_1 Z_2 \left(Z_1^{1/3} + Z_2^{1/3}\right)} \left( \frac{ \rho X_N \left<Z\right> } {\left<A\right> \rho_0 } \right)^{1/3} \]

where \( \rho_0 \equiv 1.3574 \times 10^{11} \mathrm{g/cm}^{3} \) and \( X_N \) is the total mass fraction contained in nuclei.

This function doesn’t work currently.

int test_rate()

Test the reaction rate.

int test_rate2()

Test the reaction rate.

double rate(size_t Z1, size_t Z2, size_t A1, size_t A2, double m1, double m2, double n1, double n2, double ntot, double avgZ, double avgA, double rho, double XN)

Return the rate in fm^{-3}*s^{-1}.

Densities n1, n2, and ntot should be in \( \mathrm{fm}^{-3} \), masses should be in 1/fm, rho should be in g/cm3,

As in Eq. 33 in Yakolev06 , we define

\[ R^{\mathrm{pyc}}_{ij} \equiv 8 \times 10^{7} C_{\mathrm{pyc}} \frac{\rho X_N x_i x_j A_i A_j \left<A\right> Z_i^2 Z_j^2}{(1+\delta_{ij}) (A_1+A_2)^2} S(E) {\tilde{\lambda}}^{3-C_{\mathrm{pl}}} \exp \left( -\frac{C_{\mathrm{exp}}} {\sqrt{\tilde{\lambda}}}\right)~\mathrm{fm}^{-3}~\mathrm{s} \]

int test_mixture()

Test a mixture to see that the reaction rates are reasonable.

bool is_allowed(size_t Z1, size_t Z2, size_t A1, size_t A2, double m1, double m2, double n1, double n2, double ntot, double avgZ, double avgA, double rho, double XN, double Mdot, double P)

Return true if the rate is allowed and false otherwise.

Densities n1, n2, and ntot should be in \( \mathrm{fm}^{-3} \), rho should be in g/cm3, Mdot should be in Msun/yr, mue should be in 1/fm

We always replace with a smaller rate, that is we increase Z by one if Z is odd and decrease A by one if A is odd.

If either of Z1 or Z2 are less than or equal to 4, this function always returns true.

Using the rate given in rate(), the associated fusion timescale is computed as

\[ t_{\mathrm{fusion}} = \frac{n_{\mathrm{tot}}}{3 R} \]

where \( n_{\mathrm{tot}} \) is the total number density of nuclei.

The accretion timescale is

\[ t_{\mathrm{acc}} = \frac{4 \pi R^2 \dot{M}}{\Sigma} \]

where \( \Sigma = P/g\) is the local column depth, \( P \) is the pressure and \( g \) is the surface gravity.

void get_times(size_t Z1, size_t Z2, size_t A1, size_t A2, double m1, double m2, double n1, double n2, double ntot, double avgZ, double avgA, double rho, double Xn, double Mdot, double P, double &t_fusion, double &t_acc)

Get fusion and accretion times.

Densities n1, n2, and ntot should be in \( \mathrm{fm}^{-3} \), rho should be in g/cm3, Mdot should be in Msun/yr, mue should be in 1/fm

We always replace with a smaller rate, that is we increase Z by one if Z is odd and decrease A by one if A is odd.

If either of Z1 or Z2 are less than or equal to 4, this function always returns true.

Using the rate given in rate(), the associated fusion timescale is computed as

\[ t_{\mathrm{fusion}} = \frac{n_{\mathrm{tot}}}{3 R} \]

where \( n_{\mathrm{tot}} \) is the total number density of nuclei.

The accretion timescale is

\[ t_{\mathrm{acc}} = \frac{4 \pi R^2 \dot{M}}{\Sigma} \]

where \( \Sigma = P/g\) is the local column depth, \( P \) is the pressure and \( g \) is the surface gravity.

Public Members

bool loaded

double Cpyc

Overall coefficient (default 3.9)

double Cexp

Coefficient in exponent (default 2.638)

double Cpl

Coefficient controlling power of \( \lambda \) (default 1.25)

o2scl::table b10

Table of S-factor coefficients from Beard10.

o2scl::convert_units &cng

Convert units.

bool debug

If true, use debug mode (default false)

bool allow_highZ

Allow fusions involving high Z nuclei(default true)

bool add14

If true, add fusions from Z=14 nuclei (default false)

bool use_fit

If true, use fit from Yakovlev10 (default false)

double fit[10][13]

Fit coefficients from Table II in Yakovlev10.

crust_driver

class crust::crust_driver

Main crust driver class.

Initial values for equilibrium crust

int feq_Z

Proton number.

int feq_N

Neutron number.

double feq_nn

Neutron drip density.

Baryon density grid for equilibrium crust

double feq_nb

Initial.

double feq_end

Final.

double feq_dnb

Stepsize.

The neutron and proton objects used by the ldrop_mass object

o2scl::fermion n_lda

o2scl::fermion p_lda

Tabulated EOS

table3d_eos t3d

Check parameter (default \ref check_none)

int check

const int check_none =0

const int check_mass_density =1

const int check_free_energy_sna =2

const int check_free_energy_dist =3

const int check_nm_from_dist =4

const int check_ldrop_derivs =5

const int check_mixture =6

const int check_sfactor =7

const int check_free_energy_cell =8

const int check_pressure =9

const int check_rate2 =10

const int check_mass_fit =11

const int check_feq_numbers =12

const int check_acc_numbers =13

const int check_hz_single =14

Variables needed by the reactions object \ref rn

double ec_heating

Electron capture heating fraction (default 0.25, typically 1/6 to 1/4)

bool simple_pyc

If true, use simplified pycnonuclear rates (default true)

pyc_rates pyc

Compute pycnonuclear rates.

Public Types

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

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

Public Functions

int prune_distribution(std::vector<o2scl::nucleus> &dist)

Remove nuclei in the distribution which have vanishing contribution.

Used by dist_switch_gb() and in crust::rxns for summary functions.

int check_fun(std::vector<std::string> &sv, bool itive_com)

Perform a check.

int delta_ZN(int &Z, int &N, int &deltaZ_lo, int &deltaZ_hi, int &deltaN_lo, int &deltaN_hi)

Compute loop sizes for compute_sna() or crust::rxns::gen_reaction()

bool pycno_allowed(double Z1, double Z2, double rho)

The limit for pycnonuclear fusion (old)

This function returns the upper limit in proton number for fusion given a fixed mass density

int update_flow(int Z1, int N1, int Z2, int N2, std::string type, double dheat, double rho)

Update reaction flow with a reaction between \( (Z_1,N_1) \rightarrow (Z_2,N_2) \).

This function is called by the crust::rxns object .

int check_free_energy_cell_fun()

Check for consistency of the functions free_energy_cell() and free_energy_dist()

bool dist_switch_gb(matter &m, matter &m_new, int &cnt, bool &restart, double &heat, double T, std::string type, int debug = 0)

Switch distributions and compute the heating, assuming the pressure was fixed.

crust_driver()

Constructor.

~crust_driver()

int run(int argc, char *argv[])

Basic wrapper handling command-line (called from main())

Public Members

o2scl::convert_units &cng

Convert units.

tov_shear tsh

Compute shear modulus in a neutron star crust.

dist_thermo dt

Compute the thermodyanmics of a matter distribution.

nm_thermo nmt

Properties of nuclear matter.

sna_thermo snat

Properties of crust matter with a single nucleus.

ldrop_crust lda

The mass model.

Protected Functions

int init_dist(std::string mode, matter &m)

Initialize the distribution (called by acc() )

int compute_sna(double nb, double T, matter &m, bool debug = false)

Compute the single nucleus approximation at the baryon density nb.

Called by full_eq(). Uses the given values of Z and N as initial guesses, updates these values and also returns rho and fr. Uses sna_thermo::free_energy_sna_neut and dist_thermo::mass_density().

Sets temperature to zero automatically.

int compute_sna_dist(double nb, matter &m, double T, bool debug = false)

Compute the single nucleus approximation at the baryon density nb.

Called by full_eq2(). Uses the given values of Z and N as initial guesses, updates these values and also returns rho and fr. Uses dist_thermo::free_energy_dist_sna() and dist_thermo::mass_density().

int model(std::vector<std::string> &sv, bool itive_com)

Select a model for the EOS.

int make_table(std::vector<std::string> &sv, bool itive_com)

Make an EOS table for interpolation.

int full_eq(std::vector<std::string> &sv, bool itive_com)

Compute crust in full_equilibrium.

int full_eq2(std::vector<std::string> &sv, bool itive_com)

Compute crust in full_equilibrium (V2)

This version uses free_energy_dist() instead of the “sna” functions

int full_eq_dist(std::vector<std::string> &sv, bool itive_com)

Compute crust in full_equilibrium with a distribution.

int acc(std::vector<std::string> &sv, bool itive_com)

Compute accreted crust.

int test_ndrip(std::vector<std::string> &sv, bool itive_com)

Test neutron drip.

int set_tptr(std::vector<std::string> &sv, bool itive_com)

Set temperature (takes argument in Kelvin and sets the value of Tptr in \( \mathrm{fm}^{-1} \).

Protected Attributes

bool model_set

True if the skyrme model has been set.

multi_zone mz

The object to perform multizone calculations.

double heat_full

Full heating without neutrino losses.

double rho_summary

Density at which to compute the summary objects (default \( 10^{12} \))

std::string flow_fn

Output filename for flow.

bool more_reactions

Ensure all reactions happen e.g. capturing after fusing (default true)

crust_fit cf

Object which fits the associated mass formula to experiment.

bool use_pasta

If true, include pasta deformation (default false - does not yet work)

rxns rn

Test for reactions.

bool inc_nuc_trans

If true, include the translational free energy for nuclei (default false)

std::string dist_type

Initial distribution type (default “iron”)

Either “iron”, “heavy”, “schatz” or “ashes” .

Option “ashes” gives the X-ray burst ashes from Horowitz07

bool calc_heating

If true, compute the heating rate of matter (default true)

o2scl::eos_had_temp_base *het

The hadronic eos (default points to sk)

o2scl::eos_had_skyrme sk

The homogeneous EOS (the default is SLy4)

o2scl::eos_had_apr apr

APR EOS.

o2scl::eos_had_rmf rmf

RMF EOS.

double Tptr

Temperature (in \( \mathrm{fm}^{-1} \)) (the default is \( 10^{8} \) K.

o2scl::nucmass_info nmi

To convert elements to strings.

o2scl::min_cern cm

Minimizer for full equilibrium crust (used in compute_sna() and compute_sna_dist()).

int verbose

Verbosity (default 1)

bool allow_palpha

Allow emission of protons and alpha particles (default false)

bool feq_fix_mode

(default false)

int fix_Z

Proton number for feq_fix_mode.

int fix_N

Neutron number for feq_fix_mode.

double acc_inc_factor

Factor to increase pressure or density by for accreted crust (default 1.01)

rxns

class crust::rxns

Apply nuclear reactions to a matter distribution.

Public Types

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

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

Public Functions

rxns()

int elec_capture(crust_driver *a, matter &m, matter &m_new, double T, int &cnt, double &heat)

Electron capture.

Tries to electron capture for every nucleus in the distribution.

int ec_summary(crust_driver *a, matter &m, matter &m_new, double T, ubmatrix &gpb_store, ubmatrix &m_store)

Summarize the electron capture thresholds.

int add_missing(std::vector<o2scl::nucleus> &trial, int newZ, int newN, size_t &ix)

Add a nucleus to the distribution if it’s missing, and return its index in ix.

int beta_decay(crust_driver *a, matter &m, matter &m_new, double T, int &cnt, double &heat)

Beta decay.

int emit_neutron(crust_driver *a, matter &m, matter &m_new, double T, int &cnt, double &heat)

Neutron emission.

int en_summary(crust_driver *a, matter &m, matter &m_new, double T, ubmatrix &gpb_store, ubmatrix &m_store)

Summarize the neutron emission thresholds.

int neut_capture(crust_driver *a, matter &m, matter &m_new, double T, int &cnt, double &heat)

Neutron capture.

int gen_reaction(crust_driver *a, matter &m, matter &m_new, double T, int &cnt)

Generalized reaction.

int pyc_fusion(crust_driver *a, matter &m, matter &m_new, double T, int &cnt, double &heat)

Pycnonuclear fusion.

int emit_fragment(crust_driver *a, matter &m, matter &m_new, double T, int &cnt, double &heat)

Emit proton or alpha particle.

Public Members

bool ec_one_cell

Compute electron captures by considering only one cell (default false)

bool fused_6

True if carbon has been fused.

bool fused_8

True if oxygen has been fused.

bool fused_10

True if neon has been fused.

bool fused_12

True if magnesium has been fused.

double delta_n

The fractional shift in the total number of nuclei (default 0.01)

double mdot

Mass accretion rate in solar masses per year (default \( 10^{-10} \))

shear_eigen

class crust::shear_eigen

Compute eigenfrequencies from a crust table.

Steiner & Watts

\[ \xi^{\prime \prime} + A_{\mathrm{sw}} \xi^{\prime} + B_{\mathrm{sw}} = 0 \]

\[ A_{\mathrm{sw}} \equiv \frac{\mu^{\prime}}{\rho} \frac{1}{v_s^2+v_A^2} \]

\[ B_{1,\mathrm{sw}} \equiv \omega^2 \left( 1+v_A^2\right) \]

\[ B_{2,\mathrm{sw}} \equiv \frac{-(\ell-1)(\ell+2) v_s^2}{R^2} \]

\[ B_{\mathrm{sw}} \equiv \frac{B_{1,\mathrm{sw}}+B_{2,\mathrm{sw}}} {v_s^2+v_A^2} \]

Term from Samuelsson and Andersson

\[ A_{\mathrm{sa}} \equiv \frac{d}{dr} \left\{ \log \left[ r^4 e^{\nu-\lambda} (\varepsilon + p_t) v_s^2 \right] \right\} \]

\[ B_{1,\mathrm{sa}} \equiv e^{2 \lambda-2\nu} \omega^2 \]

\[ B_{2,\mathrm{sa}} \equiv -e^{2 \lambda} v_s^2 \frac{\left(\ell-1\right)\left(\ell+2\right)}{r^2} \]

\[ B_{\mathrm{sa}} \equiv \frac{B_{1,\mathrm{sa}}+B_{2,\mathrm{sa}}}{v_s^2} \]

New expressions for Diebel, Steiner, and Brown

\[ \xi^{\prime \prime} + A_{\mathrm{dsb}} \xi^{\prime} + B_{\mathrm{dsb}} = 0 \]

\[ A_{\mathrm{dsb}} \equiv \frac{v_s^2}{v_s^2+v_A^2} \left \{ \log \left[ r^4 e^{\nu-\lambda} \left(\varepsilon+p_t\right) v_s^2 \right] \right\}^{\prime} \]

\[ B_{1,\mathrm{dsb}} \equiv \omega^2 e^{2 \lambda-2 \nu} (1+v_A^2) \]

\[ B_{2,\mathrm{dsb}} \equiv \frac{-(\ell-1)(\ell+2) v_s^2 e^{2 \lambda}} {r^2} \]

\[ B_{\mathrm{dsb}} \equiv \frac{B_{1,\mathrm{dsb}}+B_{2,\mathrm{dsb}}} {v_s^2+v_A^2} \]

Public Types

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

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

Public Functions

shear_eigen()

void store_table(std::string fname)

Store output table.

int derivs(double r, size_t nv, const ubvector &y, ubvector &dydx)

Shear frequency differential equation.

int eigen_solve(size_t nv, const ubvector &x, ubvector &y)

Solve for the eigenvalue problem.

int compute_freq(o2scl::table_units<> *shear, double Rc, double R, double ell, double &freq_lo, double freq_hi, int verbose = 0)

Compute frequency.

Public Members

bool freq_gr

If true, include GR effects in frequency calculation (default false)

o2scl::table_units *shtab

Pointer to shear table.

o2scl::ode_iv_solve ois

ODE solver.

ubvector r_grid

Grid of radii.

ubmatrix xi

Dispacement.

ubvector xi_start

Dispacement at r=R_c.

ubmatrix dxi

Derivatives.

ubmatrix xi_err

Displacement uncertainty.

double Rad

Radius.

double Radc

Radius of bottom edge of crust.

size_t ntab

Table size.

int nodes

The number of nodes.

o2scl::interp<ubvector> si

Interpolation.

tov_shear

class crust::tov_shear

Compute neutron star and shear modulus.

Public Functions

tov_shear()

int tov(std::vector<std::string> &sv, bool itive_com)

Compute crust in full_equilibrium with a distribution.

int shear(std::vector<std::string> &sv, bool itive_com)

Compute shear properties and magnetar oscillation frequencies.

Public Members

int verbose

Verbosity parameter.

double mag_field

The magnetic field in Gauss.

std::string in_dir

Desc.

o2scl::convert_units cng

Convert units.

bool freq_gr

If true, include GR effects in frequency calculation (default false)

double freq_mass

Mass for frequency calculation (default 1.4)

int hd_flag

High-density EOS flag (default 0)

Used in tov().

0 is most probable, 1 is plus 1 sigma, 2 is plus 2 sigma, 3 is minus 1 sigma, and 4 in minus 2 sigma

table3d_eos

class crust::table3d_eos

EOS from a table3d for the NS crust.

The table, tab, is created in crust_driver::make_table().

Inherits from o2scl::eos_had_temp_base

Public Functions

table3d_eos()

virtual ~table3d_eos()

virtual int calc_e(o2scl::fermion &ne, o2scl::fermion &pr, o2scl::thermo &th)

Equation of state as a function of density.

virtual int calc_temp_e(o2scl::fermion &ne, o2scl::fermion &pr, double T, o2scl::thermo &th)

A wrapper for calc_e() since this is only zero temperature.

virtual int calc_p(o2scl::fermion &ne, o2scl::fermion &pr, o2scl::thermo &th)

virtual int calc_temp_p(o2scl::fermion &ne, o2scl::fermion &pr, double T, o2scl::thermo &th)

Public Members

o2scl::table3d tab

The table which stores the EOS.