Structure of a Zero-Temperature WD

Prelude

Mathematica Magic

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Structure

Physical Constants

A = 6.002*^22; (* dynes/(cm^2) *)
B = 9.736*^5; (* g/(cm^3) *)
G = 6.672*^-8; (* (cm^3)/(s^2)/g *)
Msun = 1.989*^33; (* g *)
Rsun = 6.955*^10; (* cm *)

Equations

ϵ = (10^(-8));

CGeqns[zc_]:={D[(ζ^2)Φ'[ζ],ζ] == - (ζ^2)(((((Φ[ζ])^2)-(1/(zc^2))))^(3/2)), Φ[ϵ] == 1-((ϵ^2)/6)(((1-(1/(zc^2))))^(3/2)), Φ'[ϵ] ==-(ϵ/3) (((1-(1/(zc^2))))^(3/2))}

Solver

ZTWD[vzc_]:= Module[{s, vζ1, vΦp1},
s = NDSolve[CGeqns[vzc], Φ, {ζ, ϵ, ∞}, Method->{"EventLocator", "Event"->Re[Φ[ζ]-(1/vzc)]}];
vζ1 =InterpolatingFunctionDomain[First[Φ/.s]][[1,-1]];
vΦp1 = Re[First[Φ'[vζ1]/.s]];
{zc-> vzc, ζ1-> vζ1, Φp1 -> vΦp1, Φ -> First[Φ/.s]}]

Sanity Check

This matches Table 25.1

zcvals = (1/(({0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.5, 0.8})^(1/2)));

Table[wd = ZTWD[zcvals[[i]]];{1/(zc^2),ζ1, -((ζ1)^2)Φp1}/.wd,{i,1,8}] //TableForm

Return to Physical Quantities

M= -((4 π )/Msun)(((((2 A)/(π G))))^(3/2))((Ye^2)/(B^2))* ((ζ1)^2)*Φp1;
R=(1/Rsun) (((((2 A)/(π G))))^(1/2))(Ye/(B zc)) * ζ1;
ρc =((B )/Ye)* ((((zc^2)-1))^(3/2));

Caclculate data table

Δ = (10^(-3));
myYe = 0.5;
data = Table[vzc=(1/(x^2));wd = ZTWD[vzc]; {vzc, ζ1,Φp1, ρc, M, R} /.wd/. Ye-> myYe, {x, Δ, 1-Δ, Δ}];

ρczc[vzc_]:= ρc /. {Ye -> myYe, zc -> vzc}
Mzc[vzc_]:= Interpolation[Transpose[{Transpose[data][[1]], Transpose[data][[5]]}]][vzc]
Rzc[vzc_]:= Interpolation[Transpose[{Transpose[data][[1]], Transpose[data][[6]]}]][vzc]
dlnMdlnρ[vzc_] := ((Mzc'[vzc] * ρczc[vzc])/(Mzc[vzc] * ρczc'[vzc]))
dlnRdlnρ[vzc_] := ((Rzc'[vzc] * ρczc[vzc])/(Rzc[vzc] * ρczc'[vzc]))
dlnMdlnR[vzc_] := ((Mzc'[vzc] * Rzc[vzc])/(Mzc[vzc] * Rzc'[vzc]))
derivs =  Table[vzc=(1/(x^2));{dlnMdlnρ[vzc], dlnRdlnρ[vzc], dlnMdlnR[vzc]}, {x, Δ, 1-Δ, Δ}];

dataset = Transpose[Join[Transpose[data], Transpose[derivs]]];
Export["/tmp/ZeroTemperatureWhiteDwarf.h5", dataset]