Dla struktury gwiazd, i ogólnie w astrofizyce, bardzo ważne jest aby uwzględnić możliwość rozbicia cząstek złożonych (atomów np: H,He;  molekuł:  (Subscript[N, 2]), (Subscript[H, 2])O; jąder atomowych 4He) na składowe. Proces ten określamy mianem jonizacji lub dysocjacji. Może on zachodzić w odwrotnym kierunku, przy wyższych temperaturachprowadzac do kreacji cząstek. Najprostszy opis tego zjawiska opiera się o równowagowe obliczenia termodynamiczne. Zaczynamy od wypisania zachodzącej reakcji, na przykładzie jonizacji wodoru atomowego, który może być rozbity na składowe: elektron  i proton. Rozważamy najprostszy model tego procesu:

.

Powyższa reakcja jest równoważna równości odpowiednich potencjałów chemicznych μ

.

Będziemy potrzebować wzoru na potencjał chemiczny, na szczęście wystarczy na razie ten dla gazu doskonałego. Zastosujemy wzór na liczbę cząstek obsadzonych zgodnie z rozkładem Boltzmana

.

$Assumptions=kB>0&&m>0&&T>0

ne=Inactive[Integrate][g Exp[-(En-μ)/(kB T)] 4 Pi p^2/h^3,{p,0,Infinity}]

(*ne=Inactive[Integrate][g Exp[-En/(kB T)] Exp[μ/(kB T)]4 Pi p^2/h^3,{p,0,Infinity}]*)

Wystarczy użyć nierelatywistycznego wzoru na energię kinetyczną.

En=p^2/(2 m)

Mathematica oczywiście bez problemu policzy poniższą całkę Gaussowską, ale prawidłowa procedura polega na wyodrębnieniu istotnej części wymiarowej (stałe fizyczne) i drugorzędnej części matematycznej (wartość całki), na ogół ,,rzędu jedności’’. Dochodzimy do tego zamieniając zmienną całkowania (pęd).

calka=IntegrateChangeVariables[ne, ξ, p == Sqrt[2*m*kB*T]*ξ]

Mathematica niestety nie wyciąga przed symbol całki stałych, co dla studenta fizyki jest oczywiste. Częściowo można to zrobić jedną z poniższych, wygenerowanych przez AI reguł:

pullConst=Inactive[Integrate][c_.*f_,r__]/;FreeQ[c,First@{r}[[1]]]:>c Inactive[Integrate][f,r];(* ChatGPT o3 *)

pullConstantsRule=Inactive[Integrate][factor_*rest_,{var_Symbol,lim___}]:>factor*Inactive[Integrate][rest,{var,lim}]/;FreeQ[factor,var]; (* Gemini 2.5 Pro *)

calka//.pullConstantsRule

calka//.pullConst

Obliczamy całkę (można było to zrobić od razu).

Activate[calka]

Standardowy trik upraszczający wynikowy wzór to wprowadzenie termicznej długości fali de’Broglie:

ne=%/.h->λ Sqrt[2 Pi m kB T]

Otrzymujemy znany wzór na potencjał chemiczny gazu doskonałego. Liczba g stanów ,,wewnętrznych’’ cząstki gazu jednoatomowego to g=1; dla elektronów i protonów g=2 (dwa ustawienia spinu).

μRozw=Solve[n==ne,μ,Reals][[1]]

NA potrzeby wyprowadzenia równania Saha na stopień jonizacji, wygodnie użyć jest relatywistycznej konwencji na potencjał chemiczny μ, uwzględniającej energię spoczynkową cząstek.  O ile typowe energie jonizacji są dla wodoru o rzędy wielkości mniejsze niż energie spoczynkowe, ich różnice reprezentują energię wiązania, zgodnie ze słynnym wzorem Einsteina .  Można użyć jednej z  2 wersji równania równowagi chemicznej/kinetycznej:

. (Potencjały chemiczne zawierają )

. (Potencjały chemiczne nie zawierają , Q jest energią wiązania atomu wodoru/ energią jonizacji)

Refine[μ/.μRozw, n>0&&g>=1&&λ>0]

mu=m c^2 + %

Jesteśmy obecnie gotowi dop zapisania wyjściowego równania, które będziemy następnie chcieli rozwiązać.

Clear[ne]

μH=mH (c^2) +kB T Log[((nH ((λH)^3))/gH)];
μe=me (c^2) +kB T Log[((ne ((λe)^3))/ge)];
μp=mp (c^2) +kB T Log[((np ((λp)^3))/gp)];

Z poniższego równania chcemy wyliczyć liczbę neutralnych atomów wodoru (), protonów (jonów dodatnich ) oraz elektronów (jonów ujemnych ). Dwa brakujące równania dostarcza zasada zachowania ładunku elektrycznego (obojętność elektryczna plazmy - dogmat astrofizyki współczesnej) oraz zasada zachowania liczby barionowej.

μH==μe+μp

Zastosujemy kilka założeń upraszczających:

μH-(μe+μp) /. {λp->λH, gH->1, ge->2,gp->2, mH ->mp +me - Q/c^2}

Simplify[%]

Collect[%,{kB,T} ]

eq1=%==0

eq2=ne==np (* Neutralność ładunkowa *)

eq3=nH+np==nB (* Zachowanie liczby barionowej *)

sys={eq1,eq2,eq3}

FullSimplify[sys/.{np->x nB, nH->xH nB, ne->Ye nB},nB>0&&λH>0]

Eliminate[Eliminate[%,Ye],xH]

eqSaha=FullSimplify[%]

Poniżej standardowa forma lewej strony równania Saha. Równanie jest kwadratowe na x, stąd 2 identyczne wyrażenia.

x^2/(1-x)/.Solve[eqSaha,x]//Simplify

RHS=%[[1]]

Poniżej przystąpimy do wstawienia konkretnych wartości do równania Saha.

ElementData["Hydrogen","IonizationEnergies"] (* Jakiś błąd w Wolfram α ... )

ElementData["Helium",(EntityProperty["ionization energies", EntityProperty["Element","IonizationEnergies"], "EntityProperty[\"Element\", \"IonizationEnergies\"]"])]

Ponieważ Wolfram α ma jakiś problem, obliczymy energię jonizacji wodory ze stałęj Rydberga.

Clear[E0,me,mu,h,ρ,kB,T]

h=Quantity[1,"PlanckConstant"];

c=Quantity[1,"SpeedOfLight"];

R∞=Quantity[1,"RydbergConstant"];

E0=UnitConvert[R∞*h c,"eV"] (* Energia wiązania atomu wodoru *)

me=Quantity[1,"ElectronMass"]

mu=Quantity[1,"AtomicMassUnit"] (* Różnica pomiędzy mu, masą protonu czy atomu wodoru jest nieistontna *)

kB=Quantity[1,"BoltzmannConstant"]

λe=h/Sqrt[2 Pi me kB T]

A=RHS/.{Q->E0,nB->ρ/mu}

4 mu/ρ/λ^3 Exp[-E0/kB/T] /. {ρ->Quantity[150, "g/cc"],T->Quantity[10^7,"Kelvins"]}

A/.{ρ->Quantity[150, "g/cc"],T->Quantity[10^7,"Kelvins"]}

{xp^2/xH==A,xp+xH==1}

Clear[T,ρ]

rozwXanal=xH/.Solve[{xp^2/xH==AA,xp+xH==1},{xH,xp}][[2]]//Expand

rozwX=rozwXanal/.AA->A/.{ρ->Quantity[ρ,"kg/m^3"],T->Quantity[T,"Kelvins"]}

rozwX/.{T->10^2,ρ->10^1}

w=LogLogPlot[{rozwX/.ρ->10^5,rozwX/.ρ->1000,rozwX/.ρ->1},{T,1000,10^8},PlotRange->{10^-6,1},PlotStyle->{Red,Green, Blue},PlotLegends->{"ρ=(10^5) kg/(m^3)","ρ=1000 kg/(m^3)","ρ=1 kg/(m^3)"},AxesLabel->{"T [K]","(Subscript[X, H])"},ImageSize->300]

(*Export["JonizacjaH.png",w,ImageResolution->1200,ImageSize->1600]*)

Block[{T=10^7},LogLogPlot[rozwX,{ρ,100,10^9},PlotRange->{10^-3,1},PlotLabel->"T=(10^7)K"]]

Wpływ jonizacji na równanie stanu gazu doskonałego ilustruje wykres odwrotności średniej wagi molekularnej  w zależności od stopnia jonizacji. Wodór atomowy liczy się do ciśnienia ze współczynnikiem 1. Po jonizacji, zamienia się na 2 cząstki: proton i elektron, które dawałyby 2x większe ciśnienie przy identycznej temperaturze i gęstości gazu. 

 = 1+x

ContourPlot[(2-rozwX)/.{ρ->10^lgρ,T->10^lgT},{lgρ,-2,2},{lgT,3,6},Contours->{1.001,1.01,1.1,1.2,1.3,1.4,1.5,1.6,1.8,1.9,1.99,1.999}]

Poniżej niedokończony nieco bardziej ambitny rachunek uwzględniający występowanie w niskich temperaturach wodoru molekularnego w postaci dwuatomowej cząsteczki .

μH2=mH2 (c^2) +kB T Log[((nH2 ((λH2)^3))/gH2)];
μH=mH (c^2) +kB T Log[((nH ((λH)^3))/gH)];
μe=me (c^2) +kB T Log[((ne ((λe)^3))/ge)];
μp=mp (c^2) +kB T Log[((np ((λp)^3))/gp)];

Z poniższego równania chcemy wyliczyć liczbę neutralnych atomów wodoru (), protonów (jonów dodatnich ) oraz elektronów (jonów ujemnych ). Dwa brakujące równania dostarcza zasada zachowania ładunku elektrycznego (obojętność elektryczna plazmy - dogmat astrofizyki współczesnej) oraz zasada zachowania liczby barionowej.

reactions={μH==μe+μp,μH2==2 μH}

Zastosujemy kilka założeń upraszczających:

reactions /. {λH2->λH/Sqrt[2],λp->λH, gH->1, ge->2,gp->2,mH2->2 mH-Q2/c^2,mH ->mp +me - Q/c^2}

Simplify[%]

Collect[%,{kB,T} ]

eq1=%

eq2=ne==np (* Neutralność ładunkowa *)

eq3=2nH2+nH+np==nB (* Zachowanie liczby barionowej *)

sys=Join[eq1,{eq2,eq3}]

FullSimplify[sys/.{np->x nB, nH->xH nB, ne->Ye nB,nH2->X/2 nB},nB>0&&λH>0]

Eliminate[Eliminate[%,Ye],xH]

eqSaha2=FullSimplify[%]