Na ćwiczeniach pojawił się problem z przejściem od równania Chandrasekhara struktury białych karłów do równania Lane-Emdena z n=3/2. Podręczniki podchodzą do tej sprawy mętnie. Dla przykładu, Kippenhahn pisze:
In fact (37.8) becomes the Emden equation for
indices n= 3 and n= 3/2 if we go to the limits z->∞, and z->0,  respectively.
O ile limit n=3 jest faktycznie oczywisty ((Subscript[z, c]) jawnie znika z równania (37.8) ) to ten drugi przypadek jest wątpliwy, co poniżej wykażę. Nie zgadzają się np. warunki początkowe i brzegowe, więc aby dostać (Subscript[w, 3/2]) z φ((Subscript[z, c])->1) trzeba co najmniej dokonać przeskalowania funkcji φ, a być może także jej argumentu. Cox&Guili w swoim podręczniku podchodzą do sprawy bardziej poważnie, ale ich rozumowania nie udało się nam odtworzyć. Ciekawe co na ten temat pisze sam Chandrasekhar.

Zaczynamy od wczytania napisanych przeze mnie numerycznych implementacji f . Lane - Emdena oraz f . Chandrasekhara .

Import["http://th.if.uj.edu.pl/~odrzywolek/homepage/students/aa/Lane-Emden.m"]

Import["https://th.if.uj.edu.pl/~odrzywolek/homepage/students/stellar_evolution/Chandra.m"]

Funkcje używany następująco. Aby przyspieszyć obliczenia kod zwraca od razu całą funkcję i ją zapamiętuje.

Plot[{w[3/2][x],φ[3/2][x]},{x,0,zero[3/2]}]

Zaczniemy od weryfikacji 2 wzorów podręcznikowych: z funkcją Lane-Emdena i z f. Chandrasekhara. Powinny dać prawie identyczną strukturę białego karła.

C1=Pi me^4 c^5/3/h^3 (* Kippenhahn, równanie (15.31)*)

C2=8Pi me^3 c^3/3/h^3 μe mu(* Kippenhahn, równanie (15.15)*)

Kipp=1/20(3/Pi)^(2/3)h^2/me/(μe mu)^(5/3)//PowerExpand (* Kippenhahn, równanie (19.24)*)

Teraz sprawdzimy czy powyższy współczynnik K  dostaniemy wprost z teorii Chandrasekhara

8Integrate[ξ^4/Sqrt[1+ξ^2],{ξ,0,x}]//Expand  (* Kippenhahn, całka (15.12) oraz wzór (15.14) tak jak ją oblicza Mathematica; wzory nieco inne ale równoważne. UWAGA na 1/8 !!!! *)

Asymptotic[%,x->0] (* Przybliżenie ,,nierelatywistyczne''. *)

Wyliczenie wsp. K poprzez eliminację parametru x z parametrycznej formy równania stanu (ang. EOS), Kippenahh (37.1).

Eliminate[{P==C1 8/5 x^5,ρ==C2 x^3},x]

Solve[%,P]

K=8/5 C1 C2^(-5/3)//PowerExpand

K/Kipp (* Oba wyniki na stałą politropową K są identyczne *)

Clear[K]

G=C1=C2=1 (*Na początek pomiajmy stałe fizyczne*)

K=8/5 C1 C2^(-5/3)//PowerExpand

ρC=1/11 (* Limit niezdegenerowany odpowiada małej gęstości w centrum (Subscript[ρ, c])->0 *)

4Pi G/K/(n+1)ρC^((n-1)/n)/.n->3/2//Sqrt (* Kippenhahn, wzory (19.9) dla n=3/2*)

A=Sqrt[8Pi G/5/K ρC^(1/3)]

%/%%

ρLE=ρC w[3/2][A*r]^(3/2) (* Kippenhahn, wzór (19.11)*)

w1=Plot[Legended[Evaluate@ρLE,"Lane-Emden"],{r,0, zero[3/2]/A},PlotStyle->Red,PlotRange->{0,ρC}]

Przechodzimy do wzorów wynikających z teorii Chandrasekhara. Wstawiając do wzoru (37.10) φ=1 wynikające z warunków początkowych (37.9) dostajemy relację pomiędzy (Subscript[z, c]) oraz (Subscript[ρ, c]):

zC=Sqrt[1+(ρC/C2)^(2/3)]

zC//N

α=Sqrt[2 C1/Pi/G]/C2/zC (* Kippenhahn, wzór (37.6) *)

ρCh=C2 zC^3(φ[zC][r/α]^2-1/zC^2)^(3/2) (* Kippenhahn, wzór (37.10) *)

zero[3/2]/A (* To jest powierzchnia białego karła *)

zeroCh[zC]*α

w2=Plot[Legended[Evaluate@ρCh,"Chandrasekhar"],{r,0,zero[3/2]/A},PlotRange->{0,ρC}]

Wykres poniżej pokazuje, że struktura białego karła obliczona z r. Lane-emdena i r. Chandrasekhara (tak jak powinno być!) jest podobna dla (Subscript[ρ, c])->0. Wynikiem referencyjnym (bardziej bliskim rzeczywistości) jest oczywiście wynik Chandrasekhara.

Show[w1,w2]

(Astro)fizyka działa, POZOSTAJE matematyczne pytanie: czy funkcję φ Chandrasekhara można wyrazić przez funkcję Lane-Emdena (Subscript[w, 3/2])dla (Subscript[z, c])->(SuperPlus[1])? Innymi słowy, czy potrafimy sprowadzić równanie (37.8) do równania (19.10) z n=3/2? Zaczniemy od odtworzenia Table 37.1 z podręcznika Kippenhahn, czyli Table 25.1 str. 886 z Cox&Guili. Zwracam uwagę na idiotycznie zapisaną w obu podręcznikach ostatnią kolumnę pokazującą promień białego karł. Wynik jest oczywiście w TYSIĄCACH kilometrów. Kippenhahh całkowicie błędnie używa kropki dziesiętnej, a z kolei Cox niby traktuje , jako separator, ale w przypisie podaje dla porównania promień Ziemi BEZ separatora, a powinno być w tej samej formie czyli 6,731 km.

Clear[C1,C2]

me=QuantityMagnitude[ParticleData["Electron","Mass"],"kg"]

c=QuantityMagnitude["SpeedOfLight","m/s"]

h=QuantityMagnitude["PlanckConstant","J*s"]

mu=QuantityMagnitude["AtomicMassUnit","kg"]

G=QuantityMagnitude["GravitationalConstant","m^3/s^2/kg"]

C1=Pi me^4 c^5/3/h^3 (* Kippenhahn, równanie (15.31)*)

C2=8Pi me^3 c^3/3/h^3 μe mu(* Kippenhahn, równanie (15.15)*)

invZcTBL={0,0.01,0.02,0.05,0.1,0.2,0.3,0.5,0.8,1.0}; (* Pierwsza kolumna tabeli *)

zcTBL=1/Sqrt[invZcTBL]

xcTBL=Sqrt[zcTBL^2-1]

Tabela w pełnej okazałości.

header={"1/((Subscript[z, c])^2)\n","(Subscript[x, c])","Przecięcie z 1/(Subscript[z, c])","Nachylenie","Gęstość w centrum/(Subscript[μ, e])","Masa * ((Subscript[μ, e])^2)", "Promień białego karła R*(Subscript[μ, e]) w metrach!"};

Join[{header},Table[{1/zc^2,Sqrt[zc^2-1],zeroCh[zc],(-zeroCh[zc]^2 (φ[zc]')[zeroCh[zc]]),C2*(zc^2-1)^(3/2)/μe,μe^2 4Pi/C2^2(2 C1/Pi/G)^(3/2)(-zeroCh[zc]^2 (φ[zc]')[zeroCh[zc]]), μe zeroCh[zc]*Sqrt[2 C1/Pi/G]/C2/zc//EngineeringForm},{zc,zcTBL}]]//TableForm

Normalnie, średnia waga molekularna wynosi (Subscript[μ, e])=2 dla materii symetrycznej, t.j. składającej się z jąder atomowych posiadających tyle samo protonów i neutronów. Ważne astrofizyczne przykłady to 4He oraz 12C i 16O.

Msun=QuantityMagnitude["SolarMass","kg"]

REarth=QuantityMagnitude["EarthRadius","m"]

header2={"1/((Subscript[z, c])^2)\n","(Subscript[x, c])","Przecięcie z 1/(Subscript[z, c])","Nachylenie","Gęstość w centrum [kg/(m^3)]","Masa [(Subscript[M, ⊙])]", "Promień białego karła [(Subscript[R, ⊕])]"};

Block[{μe=2},Join[{header2},Table[{1/zc^2,Sqrt[zc^2-1],zeroCh[zc],(-zeroCh[zc]^2 (φ[zc]')[zeroCh[zc]]),C2*(zc^2-1)^(3/2), 4Pi/C2^2(2 C1/Pi/G)^(3/2)(-zeroCh[zc]^2 (φ[zc]')[zeroCh[zc]])/Msun,  zeroCh[zc]*Sqrt[2 C1/Pi/G]/C2/zc/REarth},{zc,zcTBL}]]
]//TableForm

Poniższy manipulator ilustruje fakt, że dla (Subscript[z, c])->∞ f. Ch. dąży do (Subscript[w, 3]), natomiast dla (Subscript[z, c])->(SuperPlus[1]) NIE DĄŻY do (Subscript[w, 3/2])!!!

Manipulate[Plot[{Legended[φ[zc][x],"F. Chandrasekhara"],Legended[1/zc,"1/(Subscript[z, c])"],Legended[w[3][x],"(Subscript[w, 3])"],Legended[w[3/2][x],"(Subscript[w, 3/2])"]},{x,0,1.15*zeroCh[zc]},PlotRange->{1-2(1-1/zc),1},GridLines->{{zeroCh[zc]},None},PlotStyle->{Red,{Black,DotDashed},{Blue,Dashed},{Green,Dotted},{Black,Dashed}}],{{zc,2},1.01,11,0.2}]

Wiemy, że odpowiednie przeskalowanie f. Lane-Emdena (Subscript[w, 3/2]) (n=3/2), daje f. Chandrasekchara φ[(Subscript[z, c])]. Aby zilustrować nietrywialność tego skalowana porównajmy zero miejsca zerowe. Zachowanie jest niemonotoniczne, i nie dąży do żadnej stałej dla (Subscript[z, c])->(SuperPlus[1])!

Plot[{zeroCh[ZC],zero[3/2]},{ZC,1,4}]

z[ZC_?NumericQ]:=zeroCh[ZC]-zero[3/2]

FindRoot[z,{2.2}]

FindRoot[z,{1.2}]

Również skalowanie wsp . A  i α jest nietrywialne, i jakoś kasuje powyższe  Zmienna ζ=ρ/(Subscript[C, 2]).

LogLogPlot[{Legended[2/(ζ^(-1/3)+ζ^(1/3)),"(A α)^2"],2 ζ^(1/3),2 ζ^(-1/3)},{ζ,1/100,100},PlotRange->All]

Mając na uwadze powyższe trudności, wszelkie przekształcenia analityczne będę weryfikował numerycznie. Na początek sprawdzimy, czy przekształcenie
(φ^2)-1/((Subscript[z, c])^2)=(φ-1/(Subscript[z, c]))(φ+1/(Subscript[z, c])), biorąc pod uwagę, że zarówno φ->1 jak i (Subscript[z, c])->1, czyli po prostu drugi składnik to φ+1/(Subscript[z, c])=1+1=2, mamy
(φ^2)-1/((Subscript[z, c])^2)≃2(φ-1/(Subscript[z, c])). Równanie Chandrasekhara przyjmuje przybliżoną postać:

Δφ + ((((φ^2)-1/((Subscript[z, c])^2)))^(3/2)) ≃ Δφ + ( 2(φ-1/(Subscript[z, c])) )^(3/2)=0

Bardzo ważne jest, że uproszczone równanie nadal pozostaje liniowe za sprawą nienaruszonego operatora Laplace’a Δ.

Sprawdzamy numerycznie, czy przybliżenie nic nie zepsuło:

F=f/.Block[{zc=1+1/111,ϵ=2^-26},NDSolve[{f''[x]+2/x f'[x]+Sqrt[2]^3(f[x]-1/zc)^(3/2)==0,f[ϵ]==1,f'[ϵ]==0},f,{x,ϵ,8}]][[1]]

Block[{zc=1+1/111},Plot[{φ[zc][x],F[x]},{x,0,7},PlotStyle->{Blue,Red,{Black,Dashed},Green,Purple},GridLines->{None,{1/zc}},PlotRange->{1-2(1-1/zc),1}]
]

Na razie jest idealnie, przybliżone równanie ma takie samo rozwiązanie  jak oryginał. Do f. Lane-Emdena jednak nadal daleko. Przesuwany więc funkcję o 1/(Subscript[z, c]): 
φ->φ-1/(Subscript[z, c]). Niestety, zmieniają się też warunki początkowe. Sprawdzmy efekt:

G=f/.Block[{zc=1+1/111,ϵ=2^-26},NDSolve[{f''[x]+2/x f'[x]+Sqrt[2]^3(f[x])^(3/2)==0,f[ϵ]==1-1/zc,f'[ϵ]==0},f,{x,ϵ,8}]][[1]]

Block[{zc=1+1/111},Plot[{φ[zc][x],F[x],G[x]+1/zc},{x,0,7},PlotStyle->{Blue,Red,{Black,Dashed},Green,Purple},GridLines->{None,{1/zc}},PlotRange->{1-2(1-1/zc),1}]
]

Nadal OK! Postać równania psuje czynnik ((Sqrt[2])^3) oraz złe warunki początkowe. Na początek dokonamy przeskalowania o czynnik λ, zamieniając zmienną x=λ y
Δφ/(λ^2)+((((Sqrt[2])))^3)φ=0  
czyli λ=(2^(3/4)).

H=f/.Block[{zc=1+1/111,ϵ=2^-26},NDSolve[{f''[x]+2/x f'[x]+f[x]^(3/2)==0,f[ϵ]==1-1/zc,f'[ϵ]==0},f,{x,ϵ,12}]][[1]]

Block[{zc=1+1/111},Plot[{φ[zc][x],F[x],G[x]+1/zc,H[x*2^(3/4)]+1/zc},{x,0,7},PlotStyle->{Blue,Red,{Black,Dashed},Green,Purple},GridLines->{None,{1/zc}},PlotRange->{1-2(1-1/zc),1}]
]

Ponieważ nadal nic nie udało się zepsuć, pozostaje wykonać kolejne przeskalowanie , tym razem zarówno funkcji jak i zmiennej 
φ->(Sqrt[1-(1/(Subscript[z, c]))])φ, y->κ t, co daje κ=( 1-(1/(Subscript[z, c])) )^(1/4). Ostateczny wynik to

φ(x) = (1-1/(Subscript[z, c]))(Subscript[w, 3/2])((2^(3/4))(((1-(1/(Subscript[z, c]))))^(1/4))x)+1/(Subscript[z, c])

Block[{zc=1+1/111},
Plot[
{φ[zc][x],
(1-1/zc)w[3/2][(1-1/zc)^(1/4)x*2^(3/4)]+1/zc
},{x,0,7},PlotStyle->{Red,{Black,Dashed},Green,Purple},GridLines->{None,{1/zc}},PlotRange->{1-2(1-1/zc),1}]
]