The case of an interacting relativistic gas is too difficult to discuss here. An ideal gas with a relativistic temperature is considered in the same way as in the non-relativistic case.
The distribution function is given by the standard expression (isotropy is assumed),
f(p) dp = C p^2 exp(-E/T) dp.
Here, E is the energy, C is constant, and dp is a scalar. The difference with the non-relativistic case is the relation between the energy E and the momentum p:
E = √ (p^2 c^2 + m^2 c^4).
One can rewrite f in terms of the Lorentz gamma-factor, γ = 1/√(1-v^2), p = mc v γ (velocity is in units of c):
f(γ) d γ = C γ^2 v exp(-γ/θ) d γ,
where
θ=T/mc^2.
The distribution function is normalized by the condition ∫ f d γ=1 (to avoid nonessential constants, I take a unit particle density). To calculate the integral, make the substitution
γ = cosh(x), v= √(γ^2-1)/γ = tanh(x).
Then
γ^2 v d γ = cosh(x) sinh^2(x) =
cosh^3(x) - cosh(x) = (1/4) (cosh(3x)-cosh(x)).
Use the integral formula for the Bessel function K_n(z):
K_n(z) = exp( -z cosh(x) ) cosh(n x) dx,
the integration limits are from 0 to infinity, and apply the recursion relation for the Bessel function
(2n/z) K_n (z)= K_{n+1} (z) - K_{n-1}
with n=2. This gives
C = 1/( θ K_2(1/θ) ).
The kinetic energy is equal to W = mc^2 (γ-1),
and the average kinetic energy is
= mc^2 C ∫ (γ-1) γ^2 v exp(-γ/θ) d γ
The first integral contains the factor
γ^3 v d γ = cosh^2(x) sinh^2(x) =
(1/4) sinh^2(2x) = (1/8) (cosh (4x) - 1).
Using above formulas, after some algebra we get
= mc^2 ( (3/4) K_3 - K_2 + (1/4) K_1 ) / K_2.
Arguments of Bessel functions K_n are (1/θ).
Check the low temperature limit, θ << 1. Asymptotics of Bessel functions at z>>1 are
K_n (z) -> √(π/2z) exp(-z) [ 1 + (4n^2-1)/(8z) ].
In this limit, -> mc^2 (θ/8) [ (3/4)*35 - 15 +(1/4)*3 ] =
(T/8) [ 27 - 15] = 12 T/8 = 3 T/2 ,
which coincides with the non-relativistic case.