(* Express Equation 14 in the natural basis *)

EqA1[e0_,e1_,ψ_]:=EqA2[e0,λ0]EqA2[e1,λ1]/(2Pi)
EqA2[e_,λ_]:=2e Exp[-e^2/λ]/λ

(* trivially integrate away the phase *)

f[e0_,e1_]:=Integrate[EqA1[e0,e1,ψ],{ψ,0,2Pi}]

(* Verify that total probability is unity *)

$Assumptions=e0∈Reals&&e1∈Reals&&e0 > 0&&e1>0&&ψ∈Reals&&λ0∈Reals&&λ1∈Reals&&λ0>0&&λ1>0
Integrate[f[e0,e1],{e0,0,Infinity},{e1,0,Infinity}]

(* Follow Fercher & Steeger (1981) and compute the joint probability of e0e1 where:
Λ = 2 e0 e1
d = e0
Jacobian = del(e0,e1)/del(a,d) = 1/(2d) *)

$Assumptions=Λ∈Reals&&d∈Reals&& Λ > 0&& d>0 &&e0∈Reals&&e1∈Reals&&ψ∈Reals&&λ0∈Reals&&λ1∈Reals&&λ0>0&&λ1>0 && e0 > 0&& e1>0

EqA3[Λ_]=Integrate[f[d,Λ/(2d)]/(2d),{d,0,Infinity}]

(* Now project onto one of the axes using:
s2 = Λ cosψ
b = cosψ 
Jacobian = del(Λ,cosψ)/del(s2,b) = 1/b *) 

$Assumptions=s2∈Reals&&s2>0&&b∈Reals&&b>0&&λ0∈Reals&&λ1∈Reals&&λ0>0&&λ1>0 
fb[b_]:=1/(Pi Sqrt[1-b^2])
EqA4a[s2_]=Integrate[EqA3[s2/b]fb[b]/b,{b,0,1}]

(* Note that S = 2 (Sqrt[λ0 λ1])
By symmetry, a similar result is found for s2 < 0 *)

(* Verify that total probability is unity, mean is zero, and variance is S^2/2 *)

EqA4[s2_]:=Exp[-2Abs[s2]/S]/S
$Assumptions=s∈Reals&&S∈Reals&&S>0

Integrate[EqA4[s],{s,-Infinity,Infinity}]

mean_s2 = Integrate[EqA4[s]s,{s,-Infinity,Infinity}]

var_s2 = Integrate[EqA4[s]s^2,{s,-Infinity,Infinity}]