(* Express Equation 18 in the natural basis using cylindrical coordinates *)

f[Λ_,θ_,s1_]:=Λ Exp[-2(S0 s0 - S1 s1)/Ssq]/(Pi Ssq s0)
Ssq:=S0^2-S1^2
s0:=Sqrt[Λ^2+s1^2]
$Assumptions=Λ∈Reals&&θ∈Reals&&s1∈Reals&&Λ>0&&S0∈Reals&&S1∈Reals&&0<S1<S0

(* Verify that total probability is unity *)

Integrate[f[Λ,θ,s1],{s1,-Infinity,Infinity},{Λ,0,Infinity},{θ,-Pi,Pi}]

(* Equation 20: the marginal distribution of the natural polarization s1 *)

Eq20[s1_]=Integrate[f[Λ,θ,s1],{Λ,0,Infinity},{θ,-Pi,Pi}]

(* Verify that total probability is unity, mean is S1, and variance is one half of Euclidean norm *)

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

mean_s1 = Integrate[Eq20[s]s,{s,-Infinity,Infinity}]

var_s1 = Simplify[Integrate[Eq20[s]s^2,{s,-Infinity,Infinity}]-S1^2]