ClearAll["Global`*"];
Clear["Global`*"];
CleanSlate[];
ClearInOut[];

(* An effective method to accurately calculate the phase space factors for (SuperMinus[β])(SuperMinus[β])decay 
 Andrei Neacsu and Mihai Horoi - 2015 
 Preprint available at: http://arxiv.org/abs/1510.00882 *)

(* Input example for 76Ge *)
a=76; (* the mass number *)
(Subscript[z, final])=34; (* the tomic number of the final nucleus *)
sf=94.5; (* the value of the "screening factor"; see 1510.00882, Table IV *)
J=0; (* J for the final state: valid options are J=0 for the (SuperPlus[0]) states, and J=2 for the (SuperPlus[2]) states *)
q=2.039; (* the (Subscript[Q, ββ]) value for this final state *)
(* End input *)
(* ------------------------------------------------- *)

z=(Subscript[z, final])*sf/100; (* atomic number multiplied by the "sceening factor" *)
hbc=197.33; (* ℏc *)
α=1/137.036; (* fine structure constant *)
(Subscript[m, e])=0.511; (* electron mass in MeV *)
t=q/(Subscript[m, e]); (* the kinematic energy, T=((Subscript[Q, ββ])/(Subscript[m, e]))*)
(Subscript[R, a])=(Subscript[m, e])*1.2*a^(1/3)/hbc;(* the nuclear radius multiplied by the electron mass, as commonly used in the formulae *)
(* ----------------------------------------------------- *)
(* the variables are either common to both electrons, or defined for each one of the two *)
γ[k_]:=Sqrt[k*k-(α*z)^2];  (* Eq. (A4) *)
ϵ2[ϵ1_]:=t+2-ϵ1; (* the energy of the "second" electron as function of the energy of the "first" *)
p1[ϵ1_]:=Sqrt[ϵ1^2-1];  
p2[ϵ1_]:=Sqrt[ϵ2[ϵ1]^2-1];
y1[ϵ1_]:=α*z*ϵ1/p1[ϵ1];
y2[ϵ1_]:=α*z*ϵ2[ϵ1]/p2[ϵ1];
(* -------------------------------------------------------- *)
(* The 9 kinematic factors of Eqs. (A5) *)
(Subscript[X, ββ]):=3*α*z+(Subscript[R, a])*(t+2);
b[1,ϵ1_]:=1;
b[2,ϵ1_]:=((ϵ1-ϵ2[ϵ1])^2)*0.5*(ϵ1*ϵ2[ϵ1]-1)/(ϵ1*ϵ2[ϵ1]);
b[3,ϵ1_]:=((ϵ1-ϵ2[ϵ1])^2)/(ϵ1*ϵ2[ϵ1]);
b[4,ϵ1_]:=(2/9)*(ϵ1*ϵ2[ϵ1]-1)/(ϵ1*ϵ2[ϵ1]);
b[5,ϵ1_]:=(4/3)*((Subscript[X, ββ])*(t+2)/(2*(Subscript[R, a])*ϵ1*ϵ2[ϵ1])-(ϵ1*ϵ2[ϵ1]+1)/(ϵ1*ϵ2[ϵ1]));
b[6,ϵ1_]:=4*(t+2)/((Subscript[R, a])*ϵ1*ϵ2[ϵ1]);
b[7,ϵ1_]:=(16/(3*(Subscript[R, a])*ϵ1*ϵ2[ϵ1]))*((ϵ1*ϵ2[ϵ1]+1)/(2*(Subscript[R, a]))*(Subscript[X, ββ])-t-2);
b[8,ϵ1_]:=(2/(9*(Subscript[R, a])^2*ϵ1*ϵ2[ϵ1]))*((ϵ1*ϵ2[ϵ1]+1)*((Subscript[X, ββ])^2+4*(Subscript[R, a])^2)-4*(Subscript[R, a])*(Subscript[X, ββ])*(t+2));
b[9,ϵ1_]:=(8/(Subscript[R, a])^2)*(ϵ1*ϵ2[ϵ1]+1)/(ϵ1*ϵ2[ϵ1]);
(* ------------------------------------------------ *)
(* The functions *)
ω0[ϵ1_]:=p1[ϵ1]*p2[ϵ1]*ϵ1*ϵ2[ϵ1]; (* from Eq. (A1) *)
(* Fermi function of Ec. (A3) for p1, y1 *)
(Subscript[F, 1])[k_,ϵ1_]:=(Gamma[2*k+1]/(Gamma[k]*Gamma[2*γ[k]+1]))^2*(2*p1[ϵ1]*(Subscript[R, a]))^(2*γ[k]-2*k)*(Abs[Gamma[γ[k]+I y1[ϵ1]]])^2*Exp[π*y1[ϵ1]]; 
(* Fermi function of Ec. (A3) for p2, y2 *)
(Subscript[F, 2])[k_,ϵ1_]:=(Gamma[2*k+1]/(Gamma[k]*Gamma[2*γ[k]+1]))^2*(2*p2[ϵ1]*(Subscript[R, a]))^(2*γ[k]-2*k)*(Abs[Gamma[γ[k]+I y2[ϵ1]]])^2*Exp[π*y2[ϵ1]]; 
(* The expressions that get integrated in Ec. (A1) and Ec. (B1) *)
func[0,0,n_,ϵ1_]:=b[n,ϵ1]*ω0[ϵ1]*(Subscript[F, 1])[1,ϵ1]*(Subscript[F, 2])[1,ϵ1]; (* func[0,n,ϵ1] is for 0νββ->g.s. *)
(* Notice that for 2νββ it's (Subscript[F, 1]) used twice, and not (Subscript[F, 2]). 
This is because in Eq. (B2) we use (Subscript[ϵ, 2]) defined between 0 and T+2−(Subscript[ϵ, 1]) and it is no loger a function of (Subscript[ϵ, 1]). 
As there are two intgrals for these PSF, they are significantly slower to calculate than those of the 0νββ. *)
func[1,0,n_,ϵ1_]:=p1[ϵ1]*ϵ1*(Subscript[F, 1])[1,ϵ1]*NIntegrate[(Subscript[F, 1])[1,ϵ2]*p1[ϵ2]*ϵ2*(t+2-ϵ1-ϵ2)^5,{ϵ2,1.0001,t+1.9999-ϵ1}]; (* func[2,0,ϵ1] is for 2νββ->g.s. and excited (SuperPlus[0]) *)
func[1,2,n_,ϵ1_]:=p1[ϵ1]*ϵ1*(Subscript[F, 1])[1,ϵ1]*NIntegrate[(Subscript[F, 1])[1,ϵ2]*p1[ϵ2]*ϵ2*((ϵ1-ϵ2)^2)*(t+2-ϵ1-ϵ2)^(5+2),{ϵ2,1.0001,t+1.9999-ϵ1}];(* func[2,2,ϵ1] is for 2νββ->excited (SuperPlus[2]) *)
g0[0]=2.8*10^(-22); (* (Subscript[g, 0]) for 0νββ *)
g0[1]=3.78*10^(-24)/Sqrt[7^J]; (* (Subscript[g, 0]) for 2νββ *)
(* (Subscript[a, 0ν]) of Eq. (A2) *)
PSF[n_,type_]:=(g0[type]/(Subscript[R, a])^(2-2*type))*NIntegrate[func[type,J,n,ϵ1],{ϵ1,1.0001,t+0.9999}];
(* --------------------------------------------- *)
mint=If[J<2,0,1];
Do[Print["0νββ decay -> (Subscript[J, final])=",J,"\t ",Subscript[HoldForm@"(Subscript[(G^(0  ν)), 0])",n],"=",PSF[n,type],"(yr^(-1))"],{type,mint,0},{n,1,9}]
Print["-----------------------------------------"]
Quiet[Print["2νββ decay -> (Subscript[J, final])=",J,"\t (G^(2  ν))=",PSF[1,1],"(yr^(-1))"]]

(* Output provides (Subscript[G, 01]^(0ν))-(Subscript[G, 09]^(0ν)) and (G^(2ν)) for J=0, and only (G^(2ν)) for J=2: *)