L’hamiltoniana per una particella libera relativistica è data da:
$$H = \sqrt{p_x^2 + p_y^2 + p_z^2 +m^2}$$
Per ricavare l’energia di un gas di particelle relativistico dovremmo calcolare la funzione di partizione
$$Z = \int d\Omega e^{-\beta H}$$
ma la forma dell’hamiltoniana rende estremamente difficile il calcolo di questo integrale. Invece utilizziamo il teorema di equipartizione, per cui data una variabile canonica $\xi$ abbiamo:
$$\langle \xi \pdv{H}{\xi}\rangle=k_B T$$
Il teorema può essere applicato ai tre momenti $p_x$, $p_y$ e $p_z$ separatamente, ottenendo:
$$\langle \frac{p_x^2}{H}\rangle=k_B T$$
e altre due analoghe espressioni per $p_y$ e $p_z$. Sommando le tre otteniamo:
$$\langle \frac{p_x^2+p_y^2+p_z^2}{H}\rangle=3k_B T \,\,\,\,\,\,\,\,\, \implies \,\,\,\,\,\,\,\,\, \langle H\rangle -m^2\langle \frac{1}{H}\rangle=3k_B T$$
Ora se la massa è trascurabile, ovvero nel caso cosiddetto ultrarelativistico, possiamo porre $m=0$ e quindi $\langle H\rangle=3k_B T$. Altrimenti dobbiamo ottenere un’espressione per il termine orrendo $\langle \frac{1}{H}\rangle$. Possiamo provarci calcolandone la derivata:
\begin{align*}
\pdv{}{\beta}\langle \frac{1}{H}\rangle &= \pdv{}{\beta}\bqty{\frac{1}{Z}\int d\Omega \frac{1}{H} e^{-\beta H}}=\\
&= -\frac{1}{Z^2}\pdv{Z}{\beta}\int d\Omega \frac{1}{H} e^{-\beta H}-\frac{1}{Z}\int d\Omega e^{-\beta H}=\\
&= -\pdv{\log{Z}}{\beta} \frac{1}{Z}\int d\Omega \frac{1}{H} e^{-\beta H}-1=\\
&=\langle H \rangle \langle \frac{1}{H}\rangle-1=\\
&=m^2\langle \frac{1}{H}\rangle^2+3k_B T \langle \frac{1}{H}\rangle -1
\end{align*}
Ovvero, ponendo $x = \beta$ e $y=\langle \frac{1}{H}\rangle$ abbiamo l’equazione differenziale
$$y’ = m^2 y^2 +\frac{3}{x}y-1$$
Questa è un’equazione di Riccati che può essere risolta con metodi standard. Semplicemente la mettiamo dentro Wolfram Alpha e otteniamo la soluzione:
$$y = -\frac{1}{m^2} \frac{c\bqty{\frac{1}{2}imx^2\pqty{J_1-J_3}+2x J_2}-\frac{1}{2}imx^2\pqty{Y_1 -Y_3}+2x Y_2}{x^2\pqty{cJ_2 + Y_2}}$$
dove $J_n = J_n(imx)$, $Y_n = Y_n(-imx)$ e le $J_n$ e $Y_n$ sono le funzioni di Bessel, mentre $c$ è una costante di integrazione. La formula può essere manipolata nella forma seguente:
$$m^2 y = -\frac{2}{x}-\frac{1}{2}im\frac{\pqty{cJ_1-Y_1}-\pqty{cJ_3 -Y_3}}{cJ_2 + Y_2}$$
Rimane da eliminare la costante. $y$ deve essere reale e positiva. Per ottenere una condizione necessaria calcoliamone il limite all’infinito. Utilizzando le proprietà asintotiche delle funzioni di Bessel (o anche solo Wolfram Alpha, dividendo sopra e sotto per $Y_2$) abbiamo
$$m^2 y \to -\frac{1}{2}im\frac{\pqty{c-i}-\pqty{-c +i}}{ci+1}=-im\frac{c-i}{ci+1}=-m\frac{c-i}{c-i}$$
Se $c \neq i$ allora possiamo semplificare e otteniamo che il limite è negativo; ma sappiamo che $1/H$ è positivo, per cui ciò non può succedere. Pertanto $c=i$ e quindi otteniamo invece un’espressione valida:
$$m^2 y = -\frac{2}{x}+\frac{1}{2}m\frac{J_1-J_3+iY_1-iY_3}{iJ_2+Y_2}$$
Utilizzando le proprietà delle funzioni di Bessel (date da Wolfram Alpha), questa espressione può essere ulteriormente semplificata:
$$m^2 y =m\frac{K_1(mx)}{K_2(mx)}$$
dove le $K_n$ sono le funzioni di Bessel modificate. In particolare nel nostro caso poiché $mx$ è sempre reale e positivo, entrambe $K_1$ e $K_2$ sono anch’esse reali e positive. Abbiamo inoltre le seguenti espressioni integrali:
$$K_n (mx) = \int_0^\infty e^{-mx\cosh{t}} \cosh{n t}\,dt$$
Per cui otteniamo in totale:
$$\boxed{\langle H\rangle=3k_B T+m\frac{K_1(m /k_B T)}{K_2(m/k_B T)}}$$
Utilizzando le formule standard per le equazioni di Bessel, otteniamo che per $T\to 0$ $\langle H\rangle \to m$, cioè l’unica forma di energia a temperatura nulla è la massa a riposo. Invece per $T\to \infty$ abbiamo $\langle H\rangle = 3k_B T + \mathcal{O}\pqty{\frac{1}{T}}$, cioè ad alte temperature il gas si avvicina sempre di più al gas ultrarelativistico.