Nel precedente articolo sulle stelle a neutroni abbiamo studiato numericamente l’inviluppo della stella e abbiamo scoperto che per le stelle più massicce l’inviluppo è irrilevante sia in termini di massa che di raggio. Per migliorare il limite superiore alla massa di una stella a neutroni dobbiamo quindi concentrarci sul nucleo, cioè la zona centrale dove $\rho > 5 \cdot 10^{17} \mathrm{kg/m^3}$.
Abbiamo già ottenuto un limite alla massa del nucleo utilizzando metodi puramente analitici: $5M_\odot$. Questo limite è valido qualunque sia l’equazione di stato del nucleo. Tuttavia l’equazione di stato del nucleo non può essere un’equazione qualsiasi, deve invece soddisfare una serie di requisiti. In particolare, la velocità del suono non può superare la velocità della luce, ovvero:
$$\dv{p}{\rho}\leq 1$$
dove abbiamo posto come al solito $c=G=1$. Conoscendo questo requisito e l’equazione per la massa totale del nucleo
$$m_N = \int_0^{r_N} 4\pi r^2 \rho dr$$
possiamo impostare un problema variazionale: vogliamo cioè cercare l’equazione di stato $p(\rho)$ che massimizzi la massa del nucleo. L’analisi è abbastanza complicata, al punto che non esiste tuttora un’analisi completa del problema; potete trovarne l’argomentazione nell’articolo che la presenta originariamente, Rhoades & Ruffini (1974), Maximum Mass of a Neutron Star. Viene fuori che la massa massima si ottiene quando la velocità del suono è massima, ovvero:
$$\dv{p}{\rho} = 1$$
Risolvendo l’equazione otteniamo $p = \rho + a$, dove la costante $a$ è determinata imponendo la continuità di $\rho$ e $p$ al bordo tra il nucleo e l’inviluppo, dove avremo una certa altra equazione di stato (nel nostro caso quella dei neutroni liberi).
Sostituendo nelle equazioni Tolman – Oppenheimer – Volkoff otteniamo:
\begin{align*}
\frac{dm}{dr} &= 4\pi r^2 \rho\\
\frac{d\rho}{dr} &= -(2\rho+a) \frac{m+4\pi r^3 (\rho + a)}{r(r-2m)}
\end{align*}
Queste equazioni sono ancora troppo complicate per poterle risolvere analiticamente e vanno quindi risolte numericamente.
La procedura è simile alla precedente. Scegliamo un certo valore di $\rho_N$ (di solito $\rho_N = 5 \cdot 10^{17} \mathrm{kg/m^3}$) ma dobbiamo variarlo per ottenere tutte le possibilità, dato che non è un valore esatto. A questo punto data l’equazione di stato per $\rho < \rho_N$ che abbiamo usato in precedenza e applicando la continuità al bordo del nucleo otteniamo $a = p_N -\rho_N$.
Come condizioni al contorno per l’integrazione utilizziamo $m(0)=0$ e $\rho(0)=\rho_c$ dove $\rho_c$, la densità centrale, è un parametro ignoto. Dobbiamo quindi provare un certo numero di valori per $\rho_c$ e osservarne il comportamento. L’unica condizione che abbiamo su $\rho_c$ è che sia maggiore di $\rho_N$, poiché sappiamo che la densità della stella decresce verso l’esterno. Il comportamento tipico della massa del nucleo al variare della densità centrale è evidenziato nel grafico seguente:
La massa del nucleo cresce inizialmente con la densità centrale e poi decresce verso un valore asintotico. Pertanto la procedura che seguiremo sarà la seguente: variando $\rho_N$ come nell’articolo precedente, integriamo le equazioni TOV per molti valori della densità centrale, trovando il valore massimo della massa nucleare in funzione della densità al bordo. Otteniamo un grafico del genere:
In questo caso la massa aumenta notevolmente al diminuire della densità al bordo; già nella formula analitica per la massa massima avevamo infatti ottenuto $m \propto 1/\sqrt{\rho_N}$. Si tratta quindi di decidere fino a che punto ci fidiamo dell’equazione di stato per basse densità. Scegliendo $\rho_N \approx 1,5 \cdot 10^{17} \mathrm{kg/m^3}$, che è anche circa il valore di inizio del plateau nel caso dell’articolo precedente, otteniamo per la massa massima $\approx 4,70\,\mathrm{km}$ ovvero $\boxed{\approx 3,2 M_\odot}$ che è quindi la massa massima di una stella a neutroni.
Il risultato può essere migliorato conoscendo l’equazione di stato per densità maggiori di $\rho_0\approx 1,5 \cdot 10^{17} \mathrm{kg/m^3}$. Ad esempio conoscendo l’equazione di stato tra $\rho_0$ e $\rho_1$, procederemmo all’integrazione numerica con l’equazione di stato lineare per $\rho > \rho_1$, integreremmo con l’equazione di stato nota tra $\rho_1$ e $\rho_0$ e utilizzeremmo il gas libero di elettroni per $\rho < \rho_0$. Ciò è stato fatto ad esempio da Kalogera & Baym (1996), dove il limite superiore è abbassato a circa $2,9 M_\odot$.
Tuttavia si può probabilmente fare di meglio: questo perché gli astronomi tipicamente non osservano stelle a neutroni più grandi di circa $2,1 M_\odot$: perciò questo numero è probabilmente più vicino al limite effettivo.