/* GROUND STATE COMPOSITION AND STRUCTURE OF THE OUTER CRUST OF MAGNETARS */
/* THE CODE COMPUTES THE FOLLOWING QUANTITIES: */
/* - EQUILIBRIUM NUCLEI CONSTITUTING EACH CRUSTAL LAYER, */
/* - TRANSITION PRESSURES, */
/* - DENSITIES OF ADJACENT LAYERS, */
/* - RELATIVE NUCLEAR ABUNDANCES, */
/* - DEPTHS RELATIVE TO NEUTRON-DRIP TRANSITION. */
/* REMARKS: */
/* -ONLY LANDAU-RABI QUANTIZATION OF ELECTRON MOTION INCLUDED, */
/* -ELECTRONS ARE ASSUMED TO BE ULTRARELATIVISTIC, */
/* -ELECTRONS ARE SUPPOSED TO BE IN THE LOWEST LANDAU-RABI LEVEL. */
/* execution: smagcrust <atomic_mass_data_file> */
/* options: */
/* -w Wigner-Seitz approximation */
/* The format of <atomic_mass_data_file> should be the following: */
/* atomic number / mass number / mass excess in MeV */
/* help: smagcrust -h */
/* code written by Nicolas Chamel */
/* version 22.05.2020 */
/* C headers */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
/* --- fundamental constants from NIST CODATA 2018 --- */
/* https://physics.nist.gov/cuu/Constants/ */
/* atomic mass unit in MeV/c^2 */
const double ATOMIC_MASS=931.49410242 ;
/* electron mass in MeV/c^2 */
const double ELECTRON_MASS=0.51099895 ;
/* neutron mass in MeV/c^2 */
const double NEUTRON_MASS=939.56542052 ;
/* neutron magnetic moment to nuclear magneton ratio */
const double MUN=-1.91304273 ;
/* dielectric permittivity of vacuum in units of 10^12 F/m */
const double EPS0=8.8541878128 ;
/* elementary electric charge in units of 10^-19 C */
const double E1=1.602176634 ;
/* hbar c in units MeV.fm */
const double HBARC=197.3269804 ;
/* subroutine to empty the buffer */
void empty_buffer(void)
{
int c;
while((c=getchar()) != EOF && c != '\n');
}
int main(int argc, char* argv[]) {
/* useful combinations of constants */
const double LAMBDAE=HBARC/ELECTRON_MASS ; /* electron Compton wave length */
const double E2=E1*100.0/(4.0*M_PI*EPS0) ;
const double ALPHA=E2/HBARC ; /* fine structure constant */
const double BREL=sqrt(ELECTRON_MASS/ALPHA/(LAMBDAE*LAMBDAE*LAMBDAE)*E1*1e33) ; /* characteristic magnetic field in G */
double PRESS_CONST ;
double DENS_CONST ;
/* body-centered cubic lattice constant (1) or Wigner-Seitz approximation (0) */
char BCC=1 ;
/* variables */
short NDATA ;
short *Z, *A ;
double *EXCESS ;
double MASS1,MASS2,NUCLEON_MASS0 ;
double xe,sqtxe,mue12,gam12,F12,F12bar ;
double P120,n1max0,n2min0,P122,n1max2,n2min2,P12,n1max,n2min ;
double u, gameq, gam0, gam2, gammae ;
short Z1, A1, Z2, A2 ;
double Pmin, g12, Peq, gs, gdrip ;
short i, imin ;
double y1,y2 ;
FILE *data_file ;
int read ;
char nuc=1 ;
char sol=0 ;
double xi=0.0 ;
double sum=0.0 ;
double BSTAR ;
char reliability=0 ;
int input ;
double Cion ;
double depth ;
/* check the number of arguments */
if((argc<2)||(argc>3))
{
fprintf(stderr,"The number of arguments is invalid. Use option -h for help.\a\n" );
return EXIT_FAILURE ;
}
/* print help */
if(strcmp(argv[1],"-h") == 0 || strcmp(argv[1],"-help") == 0)
{
printf("Usage:\n");
printf("%s <filename> options\n",argv[0]);
printf("where <filename> is the name of a file containing the atomic mass table.\n");
printf("Options:\n");
printf("-w Wigner-Seitz approximation\n") ;
printf("The format of the mass table should be the following:\n") ;
printf("atomic number / mass number / mass excess in MeV.\n") ;
return EXIT_SUCCESS;
}
/* read arguments */
if (argc==3)
{
if(strcmp(argv[2],"-w") == 0)
BCC=0 ;
}
/* open file */
data_file=fopen(argv[1], "r") ;
if(data_file==NULL)
{
fprintf(stderr,"The file does not exist!\a\n") ;
return EXIT_FAILURE ;
}
/* determine number of nuclei and check for errors */
i=0 ;
read=fscanf(data_file, "%hd %hd %lf",&Z1,&A1,&MASS1) ;
while(read!=EOF)
{
if((read!=3)||(Z1<1)||(A1<1))
{
fprintf(stderr,"Error when reading file at line %d!\a\n",i) ;
exit(EXIT_FAILURE) ;
}
read=fscanf(data_file, "%hd %hd %lf",&Z1,&A1,&MASS1) ;
i++ ;
}
NDATA=i ;
rewind(data_file) ;
/* allocate memory for the arrays */
Z=(short *)malloc(NDATA*sizeof(short)) ;
A=(short *)malloc(NDATA*sizeof(short)) ;
EXCESS=(double *)malloc(NDATA*sizeof(double)) ;
/* read mass excesses */
for (i=0; i<NDATA; i++)
{
if (fscanf(data_file, "%hd %hd %lf",&Z[i],&A[i],&EXCESS[i])!=3)
{
fprintf(stderr,"Error when reading mass excess at line %d!\a\n",i) ;
exit(EXIT_FAILURE) ;
}
}
fclose(data_file) ;
printf("Magnetic field strength in units of critical field %4.3le G?\n",BREL) ;
input=scanf("%lf",&BSTAR) ;
/* check the parameter */
while((input!=1)||(BSTAR<0))
{
fprintf(stderr,"Invalid argument!\a\n") ;
empty_buffer() ;
input=scanf("%lf",&BSTAR) ;
}
/* useful constants */
DENS_CONST=BSTAR/(2.0*M_PI*M_PI*LAMBDAE*LAMBDAE*LAMBDAE) ;
PRESS_CONST=0.5*DENS_CONST*ELECTRON_MASS ;
if (BCC==0)
/* Wigner-Seitz approximation */
Cion=-0.9*cbrt(4.0*M_PI/3.0) ;
else
/* body-centered cubic lattice constant */
/* from Baiko et al., Phys. Rev. E64,057402 (2001) */
Cion=-0.895929255682*cbrt(4.0*M_PI/3.0) ;
printf("\f") ;
printf("----------------------------------------------------------------------------------\n") ;
printf("GROUND-STATE COMPOSITION OF THE OUTER CRUST OF A COLD NONACCRETING MAGNETAR \n") ;
printf("WITH A MAGNETIC FIELD B=%4.3le G\n", BSTAR*BREL) ;
printf("STRONGLY QUANTIZING REGIME FOR ELECTRONS IS ASSUMED.\n") ;
printf("----------------------------------------------------------------------------------\n") ;
printf("Atomic mass table: %s\n",argv[1]) ;
printf("%hd nuclei included.\n",NDATA) ;
if (BCC==0)
printf("Wigner-Seitz approximation adopted.\n") ;
else
printf("Body-centered cubic lattice assumed.\n") ;
printf("\n") ;
/* composition of the surface layers */
Z1=26 ;
A1=56 ;
y1=(double)Z1/(double)A1 ;
MASS1=-60.607+(double)A1*ATOMIC_MASS+0.0000144381*pow(Z1,2.39)+1.55468e-12*pow(Z1,5.35) ; /* subtract electron binding energy */
NUCLEON_MASS0=MASS1/A1 ; /* mean nucleon mass in iron-56 */
Pmin=1e-2 ;
xe=pow(Z1*Z1*pow(-Cion*ALPHA,3)*BSTAR/(2.0*M_PI*M_PI),1.0/5.0) ;
gameq=sqrt(1.0+xe*xe) ;
gs=MASS1/(double)A1+y1*((gameq-1.0)*ELECTRON_MASS+4.0/3.0*Cion*E2*cbrt(Z1*Z1*DENS_CONST*xe)) ;
g12=gs ;
Peq=0.0 ;
gdrip=NEUTRON_MASS ;
/* composition of the layers below neutron drip */
printf("Z1 A1 Z2 A2 \t xe \t n1max[fm-3] \t n2min[fm-3] \t P12[MeV fm-3] \t mue12[MeV] \t g12[MeV]\t abundance \t depth\n") ;
while((g12<gdrip)&&(gameq<sqrt(1.0+2.0*BSTAR)))
{
for (i=0; i<NDATA; i++)
{
sol=0 ;
Z2=Z[i] ;
A2=A[i] ;
y2=(double)Z2/(double)A2 ;
MASS2=EXCESS[i]+(double)A2*ATOMIC_MASS+0.0000144381*pow(Z2,2.39)+1.55468e-12*pow(Z2,5.35) ; /* subtract electron binding energy */
/* rule out transitions such that Z1/A1=Z2/A2 */
/* see Appendix A of Phys.Rev.C94,065802(2016) */
if (y1!=y2)
{
/* analytical solution for the transition between (A1,Z1) and (A2,Z2) */
mue12=(MASS2/(double)A2-MASS1/(double)A1)/(y1-y2)+ELECTRON_MASS ;
gam12=mue12/ELECTRON_MASS ;
F12=((4.0*y1/3.0-y2/3.0)*cbrt(Z1*Z1)-y2*cbrt(Z2*Z2))/(y1-y2) ;
F12bar=F12*Cion/3.0*ALPHA*cbrt(BSTAR/(2.0*M_PI*M_PI)) ;
if ((F12bar>0.0)&&(gam12>0.0))
{
sol=1 ;
u=0.5*gam12/pow(sqrt(F12bar),3.0) ;
gammae=8.0*pow(sqrt(F12bar),3.0)*pow(sinh(asinh(u)/3.0),3.0) ;
}
if (F12bar<0.0)
{
u=0.5*gam12/pow(sqrt(fabs(F12bar)),3.0) ;
if (u>=1.0)
{
gammae=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cosh(acosh(u)/3.0),3.0) ;
sol=1 ;
}
if ((u<1.0)&&(u>=0.0))
{
gammae=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cos(acos(u)/3.0),3.0) ;
sol=1 ;
}
if ((u>=-1.0)&&(u<0.0))
{
gam0=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cos(acos(u)/3.0),3.0) ;
gam2=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cos(acos(u)/3.0+4.0*M_PI/3.0),3.0) ;
sol=2 ;
}
}
/* handle special cases for which two real positive mathematical solutions exist */
if (sol==2)
{
sol=0 ;
xe=sqrt(gam0*gam0-1.0) ;
n1max0=xe/y1 ;
n2min0=xe/y2*(1.0+Cion*ALPHA/3.0*(cbrt(Z1*Z1)-cbrt(Z2*Z2))*cbrt(BSTAR/(2.0*M_PI*M_PI))*gam0/pow(xe,5.0/3.0)) ;
xe=sqrt(gam2*gam2-1.0) ;
n1max2=xe/y1 ;
n2min2=xe/y2*(1.0+Cion*ALPHA/3.0*(cbrt(Z1*Z1)-cbrt(Z2*Z2))*cbrt(BSTAR/(2.0*M_PI*M_PI))*gam2/pow(xe,5.0/3.0)) ;
/* check mechanical stability of each solution separately */
if ((n2min0>n1max0)&&(n1max0>0.0)&&(gam0>=1.0))
{
sol=1 ;
/*printf("k=0 solution\n") ;*/
gammae=gam0 ;
}
if ((n2min2>n1max2)&&(n1max2>0.0)&&(gam2>=1.0))
{
/*printf("k=2 solution\n") ;*/
gammae=gam2 ;
sol++ ;
}
/* if the two solutions are still admissible, select the one with the lowest transition pressure */
if (sol==2)
{
/*printf("Two real solutions found!\n") ;*/
xe=sqrt(gam0*gam0-1.0) ;
P120=PRESS_CONST*(xe*gam0-log(xe+gam0))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
xe=sqrt(gam2*gam2-1.0) ;
P122=PRESS_CONST*(xe*gam2-log(xe+gam2))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
if (P120<P122)
gammae=gam0 ;
else
gammae=gam2 ;
sol=1 ;
}
}
/* if a solution to the transition between (A1,Z1) and (A2,Z2) exists */
if (sol==1)
{
xe=sqrt(gammae*gammae-1.0) ;
n1max=DENS_CONST*xe/y1 ;
n2min=DENS_CONST*xe/y2*(1.0+Cion*ALPHA/3.0*(cbrt(Z1*Z1)-cbrt(Z2*Z2))*cbrt(BSTAR/(2.0*M_PI*M_PI))*gammae/pow(xe,5.0/3.0)) ;
/* check mechanical stability */
if ((n2min>n1max)&&(n1max>0.0)&&(gammae>=1.0))
{
P12=PRESS_CONST*(xe*gammae-log(xe+gammae))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
/* find lowest transition pressure */
/* ensure the pressure increases with density */
if ((P12<Pmin)&&(P12>=Peq))
{
Pmin=P12 ;
gameq=gammae ;
imin=i ;
}
}
}
}
}
Z2=Z[imin] ;
A2=A[imin] ;
y2=(double)Z2/(double)A2 ;
MASS2=EXCESS[imin]+(double)A2*ATOMIC_MASS+0.0000144381*pow(Z2,2.39)+1.55468e-12*pow(Z2,5.35) ;
mue12=(MASS2/(double)A2-MASS1/(double)A1)/(y1-y2)+ELECTRON_MASS ;
gam12=mue12/ELECTRON_MASS ;
xe=sqrt(gameq*gameq-1.0) ;
n1max=DENS_CONST*xe/y1 ;
n2min=DENS_CONST*xe/y2*(1.0+Cion*ALPHA/3.0*(cbrt(Z1*Z1)-cbrt(Z2*Z2))*cbrt(BSTAR/(2.0*M_PI*M_PI))*gameq/pow(xe,5.0/3.0)) ;
P12=PRESS_CONST*(xe*gameq-log(xe+gameq))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
g12=MASS1/(double)A1+y1*((gameq-1.0)*ELECTRON_MASS+4.0/3.0*Cion*E2*cbrt(Z1*Z1*DENS_CONST*xe)) ;
if ((g12<gdrip)&&(gameq<sqrt(1.0+2.0*BSTAR)))
{
nuc++ ;
sum+=xi ;
xi=P12-sum ;
Peq=P12 ;
depth=(pow(g12/gs,2)-1.0)/(pow(NEUTRON_MASS/gs,2)-1.0) ;
/* check reliability of the solution */
if (gam12<=0.0)
{
reliability=1 ;
printf("%hd %hd %hd %hd \t %4.2f \t %4.3le \t %4.3le \t %4.3le \t %4.2f* \t %6.3f \t %3.2le \t %5.3f\n",Z1,A1,Z2,A2,xe,n1max,n2min,P12,mue12,g12,xi,depth) ;
}
/* normal case */
if (gam12>0.0)
printf("%hd %hd %hd %hd \t %4.2f \t %4.3le \t %4.3le \t %4.3le \t %4.2f \t\t %6.3f \t %3.2le \t %5.3f\n",Z1,A1,Z2,A2,xe,n1max,n2min,P12,mue12,g12,xi,depth) ;
Z1=Z2 ;
A1=A2 ;
y1=y2 ;
MASS1=MASS2 ;
Pmin=1e-2 ;
}
}
/* check strongly quantizing regime */
if (gameq>=sqrt(1.0+2.0*BSTAR))
{
printf("Calculation stopped: lowest Landau-Rabi level completely filled.\a\n") ;
printf("Abundances cannot be normalized!\n") ;
}
if (gameq<sqrt(1.0+2.0*BSTAR))
{
/* properties of the neutron-drip transition */
mue12=(-MASS1+(double)A1*gdrip)/(double)Z1+ELECTRON_MASS ; /* always positive ?? */
gam12=mue12/ELECTRON_MASS ;
F12bar=4.0/3.0*Cion/3.0*ALPHA*cbrt(Z1*Z1*BSTAR/(2.0*M_PI*M_PI)) ; /* always negative */
u=0.5*gam12/pow(sqrt(fabs(F12bar)),3.0) ;
if (u>=1.0)
gameq=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cosh(acosh(u)/3.0),3.0) ;
if ((u<1.0)&&(u>=0.0))
gameq=8.0*pow(sqrt(fabs(F12bar)),3.0)*pow(cos(acos(u)/3.0),3.0) ;
xe=sqrt(gameq*gameq-1.0) ;
n1max=DENS_CONST*xe/y1 ;
P12=PRESS_CONST*(xe*gameq-log(xe+gameq))+Cion*E2/3.0*cbrt(pow(DENS_CONST*xe,4)*Z1*Z1) ;
sum+=xi ;
xi=P12-sum ;
depth=1.0 ;
printf("%hd %hd - - \t %4.2f \t %4.3le \t - \t\t %4.3le \t %4.2f \t\t %6.3f \t %3.2le \t %5.3f\n",Z1,A1,xe,n1max,P12,mue12,gdrip,xi,depth) ;
printf("The outer crust is stratified into %d layers.\n",nuc) ;
printf("Warning: nuclear abundances are not normalized. Their sum is not 1 but %4.3le.\n", sum+xi) ;
printf("Depths are relative to the neutron-drip point.\n") ;
}
if(reliability==1)
printf("* Warning: this solution is not expected to be accurate.\a\n") ;
free(Z) ;
Z=NULL ;
free(A) ;
A=NULL ;
free(EXCESS) ;
EXCESS=NULL ;
return EXIT_SUCCESS ;
}