Particle density

n[m_,T_,a_,g_] := g T m^2/(2 Pi^2) Sum[a^(k+1) / k  BesselK[2,k m /T],{k,1,20}]   (* a=1 for bosons and -1 for fermions *)

Ratios of particles

mHe3= 2808.4;
mHe4= 3727.4;
md=1875.6 ; (* Masses and temperature are in MeV *)
T=156.5 ;

Deuteron to He3 ratio

dHe3 = 3n[md,T,1,1]/(2n[mHe3,T,-1,1])

He3 to He4 ratio

He3He4 = 2 n[mHe3,T,-1,1]/n[mHe4,T,1,1]

Particle density vs. mass

(* We compute the particle density per degree of freedom and we use three "methods" *)

Quantum statistics  (bosons and fermions)

Plot[{n[m,T,-1,1],n[m,T,+1,1]},{m,0,1000},AxesLabel-> {"m (MeV)","n/g"},PlotLegends->{Fermion, Boson}]

Boltzmann approximation

nboltz[m_,T_,g_] := g /(2 Pi^2) m^2 T BesselK[2,m/T];

Plot[nboltz[m,T,1],{m,0,1000},AxesLabel-> {"m (MeV)","n/g"},PlotLegends->{"Boltzmann app."}]

Boltzmann approximation with large-argument approximation at leading order

nlargem[m_,T_,g_]:=  g /(2 Pi^2) m^2 T Sqrt[T Pi/(2m)] Exp[-m/T] ;

Plot[nlargem[m,T,1],{m,0,1000},AxesLabel-> {"m (MeV)","n/g"},PlotLegends->{"Large mass app."}]

All together now!

Plot[{n[m,T,-1,1],n[m,T,+1,1],nboltz[m,T,1],nlargem[m,T,1]},{m,0,1000},AxesLabel-> {"m (MeV)","n/g"},PlotLegends->{Fermion, Boson,"Boltzmann app.", "Large mass app."}]