Neutron Star Outer Crust

The class o2scl::eos_crust computes the composition and EOS of the outer crust of a neutron star from $ 10^{6}~\mathrm{g}/\mathrm{cm}^3 $ up to the neutron drip density of $ 4 \times 10^{11}~\mathrm{g}/\mathrm{cm}^3 $ . It uses a simple expression for the lattice energy and works with any nuclear mass formula specified as an object of type o2scl::nucmass .

Outer crust example

This example computes the outer crust EOS from the o2scl::eos_crust class and then compares it with the numerical values from the original paper in Baym71tg as stored in the file examples/ex_bps_input.txt. Then it computes several crusts from Shen11b, Steiner12, and Newton13.

/* Example: ex_eos_crust.cpp
-------------------------------------------------------------------
Compute the Baym-Pethick-Sutherland equation of state
*/
#include <fstream>
#include <o2scl/test_mgr.h>
#include <o2scl/lib_settings.h>
#include <o2scl/table_units.h>
#include <o2scl/eos_crust.h>
#include <o2scl/eos_tov.h>
#include <o2scl/hdf_io.h>
using namespace std;
using namespace o2scl;
using namespace o2scl_const;
using namespace o2scl_hdf;
// A simple function to load the BPS data without
// HDF support, in case O2scl was compiled without it.
int bps_load(table_units<> &bps);
int main(void) {
cout.setf(ios::scientific);
// Energy, pressure, etc.
thermo th;
// Proton and mass numbers
int Z, A;
// The baryon number density
double barn;
// Energy density from BPS table in units of g/cm^3
double ed_bps;
// Pressure from BPS table in units of dyn/cm^2
double pr_bps;
// Read the original BPS results
bps_load(bps);
// Ensure that this works without GNU units
cu.use_gnu_units=false;
// Cubic spline interpolation doesn't do very well here, so we use
// linear interpolation to interpolate the BPS table
bps.set_interp_type(itp_linear);
// Tabulate the o2scl results using eos_crust::calc_pressure() and
// compare them to the original table. The energy density
// and pressure are output in MeV/fm^3, the baryon number density
// is given in fm^{-3}. The last two columns are the fractional
// deviations from the original BPS table.
cout << "Default mass formula, using calc_pressure()." << endl;
cout << "eden pres barn A Z ";
cout << "ed_err pr_err" << endl;
for(double pr=8.9e-4/hc_mev_fm;pr>=2.12e-10/hc_mev_fm;pr/=2.0) {
th.pr=pr;
be.calc_pressure(th,barn,Z,A);
// Output energy and pressure in MeV/fm^3
cout << th.ed*hc_mev_fm << " " << th.pr*hc_mev_fm
<< " " << barn << " ";
cout.width(3);
cout << A << " ";
cout.width(3);
cout << Z << " ";
// Convert the BPS results from g/cm^3 and dyn/cm^2 to MeV/fm^3
ed_bps=cu.convert("g/cm^3","MeV/fm^3",
bps.interp("nb",barn*1.0e39,"rho"));
pr_bps=cu.convert("dyne/cm^2","MeV/fm^3",
bps.interp("nb",barn*1.0e39,"P"));
cout.setf(ios::showpos);
cout << (ed_bps-th.ed*hc_mev_fm)/ed_bps << " "
<< (pr_bps-th.pr*hc_mev_fm)/pr_bps << endl;
cout.unsetf(ios::showpos);
t.test_rel((ed_bps-th.ed*hc_mev_fm)/ed_bps,0.0,2.0e-3,"Energy");
t.test_rel((pr_bps-th.pr*hc_mev_fm)/pr_bps,0.0,1.5e-1,"Pressure");
}
cout << endl;
// Tabulate the o2scl results using eos_crust::calc_density() and
// compare them to the original table
cout << "Default mass formula, using calc_density()." << endl;
cout << "eden pres barn A Z ";
cout << "ed_err pr_err" << endl;
for(barn=4.19e-4;barn>=6.39e-9;barn/=2.0) {
be.calc_density(barn,th,Z,A);
ed_bps=cu.convert("g/cm^3","MeV/fm^3",
bps.interp("nb",barn*1.0e39,"rho"));
pr_bps=cu.convert("dyne/cm^2","MeV/fm^3",
bps.interp("nb",barn*1.0e39,"P"));
cout << th.ed*hc_mev_fm << " " << th.pr*hc_mev_fm
<< " " << barn << " ";
cout.width(3);
cout << A << " ";
cout.width(3);
cout << Z << " ";
cout.setf(ios::showpos);
cout << (ed_bps-th.ed*hc_mev_fm)/ed_bps << " "
<< (pr_bps-th.pr*hc_mev_fm)/pr_bps << endl;
cout.unsetf(ios::showpos);
t.test_rel((ed_bps-th.ed*hc_mev_fm)/ed_bps,0.0,2.0e-3,"Energy");
t.test_rel((pr_bps-th.pr*hc_mev_fm)/pr_bps,0.0,1.5e-1,"Pressure");
}
cout << endl;
// ---------------------------------------------------------
// Create a table comparing the crust EOSs from eos_tov_interp
// Read the APR EOS table
string name;
hf.open("ex_eos_had_apr_nstar.o2");
hdf_input(hf,tab,name);
hf.close();
// Send the EOS to the eos_tov_interp object
te.read_table(tab,"ed","pr","nb");
// Define a grid of pressures over which to evaluate the EOS
size_t ngrid=200;
uniform_grid_log_end<double> pr_grid(2.0e-14,4.0e-3,ngrid-1);
// Create a table to store the data
table_units<> crust_comp;
crust_comp.line_of_names("pr ed_NVBPS ed_SHO ed_PNM_L40 ed_PNM_L100");
crust_comp.line_of_names("ed_J35_L40 ed_J35_L100 ed_SLy4 ed_APR ed_Rs");
crust_comp.line_of_names("nb_NVBPS nb_SHO nb_PNM_L40 nb_PNM_L100");
crust_comp.line_of_names("nb_J35_L40 nb_J35_L100 nb_SLy4 nb_APR nb_Rs");
crust_comp.set_unit("pr","1/fm^4");
crust_comp.set_unit("ed_NVBPS","1/fm^4");
crust_comp.set_unit("ed_SHO","1/fm^4");
crust_comp.set_unit("ed_PNM_L40","1/fm^4");
crust_comp.set_unit("ed_PNM_L100","1/fm^4");
crust_comp.set_unit("ed_J35_L40","1/fm^4");
crust_comp.set_unit("ed_J35_L100","1/fm^4");
crust_comp.set_unit("ed_APR","1/fm^4");
crust_comp.set_unit("ed_SLy4","1/fm^4");
crust_comp.set_unit("ed_Rs","1/fm^4");
crust_comp.set_unit("nb_NVBPS","1/fm^3");
crust_comp.set_unit("nb_SHO","1/fm^3");
crust_comp.set_unit("nb_PNM_L40","1/fm^3");
crust_comp.set_unit("nb_PNM_L100","1/fm^3");
crust_comp.set_unit("nb_J35_L40","1/fm^3");
crust_comp.set_unit("nb_J35_L100","1/fm^3");
crust_comp.set_unit("nb_APR","1/fm^3");
crust_comp.set_unit("nb_SLy4","1/fm^3");
crust_comp.set_unit("nb_Rs","1/fm^3");
// Energy density and baryon density
double ed, nb;
cout << "NVBPS crust." << endl;
crust_comp.set_nlines(200);
for(size_t i=0;i<200;i++) {
crust_comp.set("pr",i,pr_grid[i]);
te.ed_nb_from_pr(cu.convert("1/fm^4","Msun/km^3",pr_grid[i]),ed,nb);
ed=cu.convert("Msun/km^3","1/fm^4",ed);
crust_comp.set("ed_NVBPS",i,ed);
crust_comp.set("nb_NVBPS",i,nb);
if (i==199) {
cout << pr_grid[i] << " " << ed << " " << nb << endl;
}
}
cout << "SHO crust." << endl;
for(size_t i=0;i<200;i++) {
crust_comp.set("pr",i,pr_grid[i]);
te.ed_nb_from_pr(cu.convert("1/fm^4","Msun/km^3",pr_grid[i]),ed,nb);
ed=cu.convert("Msun/km^3","1/fm^4",ed);
crust_comp.set("ed_SHO",i,ed);
crust_comp.set("nb_SHO",i,nb);
if (i==199) {
cout << pr_grid[i] << " " << ed << " " << nb << endl;
}
}
cout << "NGL13, L=40 crust." << endl;
for(size_t i=0;i<200;i++) {
crust_comp.set("pr",i,pr_grid[i]);
te.ed_nb_from_pr(cu.convert("1/fm^4","Msun/km^3",pr_grid[i]),ed,nb);
ed=cu.convert("Msun/km^3","1/fm^4",ed);
crust_comp.set("ed_PNM_L40",i,ed);
crust_comp.set("nb_PNM_L40",i,nb);
}
cout << "NGL13, L=100 crust." << endl;
te.ngl13_low_dens_eos(100.0);
for(size_t i=0;i<200;i++) {
crust_comp.set("pr",i,pr_grid[i]);
te.ed_nb_from_pr(cu.convert("1/fm^4","Msun/km^3",pr_grid[i]),ed,nb);
ed=cu.convert("Msun/km^3","1/fm^4",ed);
crust_comp.set("ed_PNM_L100",i,ed);
crust_comp.set("nb_PNM_L100",i,nb);
}
cout << "NGL13, L=40 crust." << endl;
te.ngl13_low_dens_eos(40.0,"J35");
for(size_t i=0;i<200;i++) {
crust_comp.set("pr",i,pr_grid[i]);
te.ed_nb_from_pr(cu.convert("1/fm^4","Msun/km^3",pr_grid[i]),ed,nb);
ed=cu.convert("Msun/km^3","1/fm^4",ed);
crust_comp.set("ed_J35_L40",i,ed);
crust_comp.set("nb_J35_L40",i,nb);
}
cout << "NGL13, L=100 crust." << endl;
te.ngl13_low_dens_eos(100.0,"J35");
for(size_t i=0;i<200;i++) {
crust_comp.set("pr",i,pr_grid[i]);
te.ed_nb_from_pr(cu.convert("1/fm^4","Msun/km^3",pr_grid[i]),ed,nb);
ed=cu.convert("Msun/km^3","1/fm^4",ed);
crust_comp.set("ed_J35_L100",i,ed);
crust_comp.set("nb_J35_L100",i,nb);
}
cout << "S12, SLy4 crust." << endl;
te.s12_low_dens_eos("SLy4");
for(size_t i=0;i<200;i++) {
crust_comp.set("pr",i,pr_grid[i]);
te.ed_nb_from_pr(cu.convert("1/fm^4","Msun/km^3",pr_grid[i]),ed,nb);
ed=cu.convert("Msun/km^3","1/fm^4",ed);
crust_comp.set("ed_SLy4",i,ed);
crust_comp.set("nb_SLy4",i,nb);
}
cout << "S12, APR crust." << endl;
te.s12_low_dens_eos("APR");
for(size_t i=0;i<200;i++) {
crust_comp.set("pr",i,pr_grid[i]);
te.ed_nb_from_pr(cu.convert("1/fm^4","Msun/km^3",pr_grid[i]),ed,nb);
ed=cu.convert("Msun/km^3","1/fm^4",ed);
crust_comp.set("ed_APR",i,ed);
crust_comp.set("nb_APR",i,nb);
}
cout << "S12, Rs crust." << endl;
te.s12_low_dens_eos("Rs");
for(size_t i=0;i<200;i++) {
crust_comp.set("pr",i,pr_grid[i]);
te.ed_nb_from_pr(cu.convert("1/fm^4","Msun/km^3",pr_grid[i]),ed,nb);
ed=cu.convert("Msun/km^3","1/fm^4",ed);
crust_comp.set("ed_Rs",i,ed);
crust_comp.set("nb_Rs",i,nb);
}
cout << endl;
hf.open_or_create("ex_crust_comp.o2");
hdf_output(hf,crust_comp,"crust_comp");
hf.close();
t.report();
return 0;
}
// End of example
o2scl::uniform_grid_log_end
thermo_tl< double >::ed
double ed
o2scl_const
hdf_output
void hdf_output(hdf_file &hf, o2scl::uniform_grid< double > &h, std::string name)
o2scl::eos_crust
Baym-Pethick-Sutherland equation of state.
Definition: eos_crust.h:84
o2scl::eos_tov_interp::read_table
void read_table(table_units<> &eosat, std::string s_cole, std::string s_colp, std::string s_colnb="")
Specify the EOS through a table.
o2scl::table_units::set_unit
void set_unit(std::string scol, std::string unit)
o2scl::test_mgr::report
bool report() const
thermo_tl< double >
o2scl::lib_settings_class::get_convert_units
convert_units< double > & get_convert_units()
o2scl_hdf::hdf_file::close
void close()
hdf_input
void hdf_input(hdf_file &hf, o2scl::table< vec_t > &t, std::string name)
o2scl::table_units
o2scl_hdf::hdf_file::open
int open(std::string fname, bool write_access=false, bool err_on_fail=true)
o2scl_settings
lib_settings_class o2scl_settings
o2scl::test_mgr::test_rel
bool test_rel(data_t result, data_t expected, data_t rel_error, std::string description)
hc_mev_fm
const double hc_mev_fm
o2scl::convert_units::convert
virtual fp_t convert(std::string from, std::string to, fp_t val)
thermo_tl< double >::pr
double pr
o2scl::test_mgr
o2scl_hdf
Additional functions to read and write EOS data to HDF5 files.
o2scl::eos_crust::calc_density
virtual int calc_density(double barn, thermo &th, int &Z, int &A)
Calculate the equation of state as a function of the baryon number density barn.
o2scl::eos_crust::calc_pressure
virtual int calc_pressure(thermo &th, double &barn, int &Z, int &A)
Calculate the equation of state as a function of the pressure.
o2scl::test_mgr::set_output_level
void set_output_level(int l)
o2scl::eos_tov_interp::ed_nb_from_pr
virtual void ed_nb_from_pr(double pr, double &ed, double &nb)
Given the pressure, produce the energy and number densities.
o2scl_hdf::hdf_file
o2scl_hdf::hdf_file::open_or_create
void open_or_create(std::string fname)
o2scl::eos_tov_interp
An EOS for the TOV solver using simple linear interpolation and an optional crust EOS.
Definition: eos_tov.h:741
o2scl::convert_units< double >
o2scl::eos_tov_interp::sho11_low_dens_eos
void sho11_low_dens_eos()
Crust EOS from Shen11b.
o2scl::convert_units::use_gnu_units
bool use_gnu_units
barn
const double barn
o2scl::eos_tov_interp::s12_low_dens_eos
void s12_low_dens_eos(std::string model="SLy4", bool external=false)
Crust EOS from Steiner12.
o2scl::eos_tov_interp::ngl13_low_dens_eos
void ngl13_low_dens_eos(double L, std::string model="PNM", bool external=false)
Crust EOS from Newton13 given L in MeV.

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).