Neutron freeze out

(* fundamental constants *)
Clear[fund,c,ℏ,G,J,Q];
fund={
(* speed of light *)
c->299792458(*m/s*),
(* Planck's constant / 2π *)
ℏ->1.05457180 (10^(-34))(*J s*),
(* Newton's gravitational constant *)
G->6.67384 (10^(-11))(* (m^3)(kg^(-1))(s^(-2)) *),
(* Joule in MeV *)
J->6.241509126 (10^18)/(10^6)(*MeV*),
(* neutron-proton mass difference *)
Q->1.293332(*MeV*)
}

(* radiation-dominated Hubble parameter at T = Q, in (s^(-1)) *)
Clear[H];
H=(Sqrt[((gs 4(π^3)G(((Q/J))^4))/(45(c^5)(ℏ^3)))])/.gs->2+(7/8)10/.fund

(* check units of H *)
Simplify[(Sqrt[((((*G*)(m^3)(kg^(-1))(s^(-2)))((((*Q*)kg (m^2)(s^(-2))))^4))/(((((*c*)m (s^(-1))))^5)((((*ℏ*)kg (m^2)(s^(-1))))^3)))]),s>0]

Clear[amin,amax,τn,λnp];
(* minimum and maximum values of cosmic scale factor, normalized to a = 1 at T = Q *)
amin=(10^(-2));amax=(10^2);
(* neutron lifetime in (s^(-1)) *)
τn=877.7;
(* n -> p rate in (s^(-1)) *)
λnp[a_]=(255/(τn (a^5)))(12+6a+(a^2));

λnp[1]

(* solve differential equation for neutron fraction Xn == ((Subscript[n, n])/((Subscript[n, p])+(Subscript[n, n]))) *)
Clear[Xn,XnSoln];
XnSoln=NDSolve[{(∂_a)Xn[a]==((a λnp[a])/H)(-Xn[a]+(1-Xn[a])(E^(-a))),Xn[amin]==(1/2)},Xn[a],{a,amin,amax}][[1]]

(* plot solution for neutron fraction *)
LogLogPlot[Xn[a]/.XnSoln,{a,amin,amax},PlotRange->{.1,.5}]

(* final value of neutron fraction *)
Xn[a]/.XnSoln/.a->amax