Recap on I. to IV.¶
V. Neutron stars with protons and electrons, fermi gas equation of state¶
Why must here be $p$ and $e^-$ in a neutron star?
- free $n$ undergo weak decay $n \rightarrow p + e^- + \overline{\nu}_e$ with a lifetime of ~15 min
- presence of $p$ and $e^-$ needed to prevent decay of "neutron star"
- Protons, how?
- Decay products have energy ~0.78 MeV, mostly carried by $e^-$ and $\overline{\nu}_e$.
- Once low-energy proton levels are filled, Pauli principle inhibits the decay.
- Electrons, how?
- To cancel out the positive charge of protons, hence $$ k_{F, p} = k_{F, e} \tag{50} $$
- Weak interaction equilibrium: $p + e^- \rightarrow n + \nu_e$
In other words, we expect the following equilibrium: $$ \mu_n = \mu_p + \mu_e \tag{49} $$
((const.m_n - const.m_p - const.m_e) * const.c**2).to('MeV')
We may rewrite the chemical potentials in terms of Fermi momentums and masses, and work out ($c \equiv 1$) $$\begin{aligned} k_{F, p}(k_{F, n}) &= \frac{1}{2}\sqrt{\frac{(k_{F, n}^2 + m_n^2 - m_e^2)^2 - 2m_p^2(k_{F, n}^2 + m_n^2 + m_e^2) + m_p^4}{k_{F, n}^2 + m_n^2}} \\ &\approx \frac{k_{F, n}^2 + m_n^2 - m_p^2}{2\sqrt{k_{F, n}^2 + m_n^2}} \hspace{20pt}\mathrm{for}~\frac{m_e}{k_{F, n}}\rightarrow 0 \ . \end{aligned} \tag{55}$$
In terms of $x \equiv k_F/mc$, it reads $$ x_p(x_n) = x_e(x_n) \approx \frac{x_n^2 + 1 - (m_p/m_n)^2}{2\sqrt{x_n^2 + 1}} \ . $$
def fermi_momentum_of_proton_or_electron(x):
return (x**2 + 1 - (const.m_p / const.m_n)**2) / (2 * np.sqrt(x**2 + 1))
Energy density and pressure are given by $$ \varepsilon_i(x_i) = \varepsilon_0\cdot3\int_0^{x_i} \left[u^2+(m_i/m_n)^2\right]^{1/2} u^2 \ du \equiv \varepsilon_0\overline{\varepsilon} \tag{57} $$ and $$ p_i(x_i) = \varepsilon_0\cdot\int_0^{x_i} \left[u^2+(m_i/m_n)^2\right]^{-1/2} u^4 \ du \equiv \varepsilon_0\overline{p} \tag{62} $$ where $i = n, p, e$ labels the particles, and $$\varepsilon_0 \equiv \frac{m_n^4c^5}{3\pi^2\hbar^3} \tag{58} $$
epsilon_0 = (const.m_n**4 * const.c**5) / (3 * np.pi**2 * const.hbar**3)
epsilon_0.to('erg/cm^3')
def energy_density(x, particle, include_constant=True):
particle_mass = const.__dict__[f'm_{particle}']
mass_ratio = (particle_mass / const.m_n).value
integrand = lambda u: (u**2 + mass_ratio**2)**0.5 * u**2
if particle != 'n':
x = fermi_momentum_of_proton_or_electron(x)
integral = 3*quad(integrand, 0, x)[0]
return epsilon_0 * integral if include_constant else integral
def pressure(x, particle, include_constant=True):
particle_mass = const.__dict__[f'm_{particle}']
mass_ratio = (particle_mass / const.m_n).value
integrand = lambda u: (u**2 + mass_ratio**2)**(-0.5) * u**4
if particle != 'n':
x = fermi_momentum_of_proton_or_electron(x)
integral = quad(integrand, 0, x)[0]
return epsilon_0 * integral if include_constant else integral
The total energy density and total pressure are then given by summing of the three species of particles: $$ \varepsilon_\mathrm{tot}(x_n) = \sum_{i=n,p,e} \varepsilon_i(x_i) \tag{56} $$ $$ p_\mathrm{tot}(x_n) = \sum_{i=n,p,e} p_i(x_i) \tag{61} $$ Notice the dependency difference here. The L.H.S. depends on $x_n$, so when calculating the sum at the R.H.S., we need to express $x_p$ and $x_e$ in terms of $x_n$ according to the relation we worked out earlier: $$ x_p(x_n) = x_e(x_n) \approx \frac{x_n^2 + 1 - (m_p/m_n)^2}{2\sqrt{x_n^2 + 1}} \ . $$
def total_energy_density(x_n, include_constant=True):
x_p = fermi_momentum_of_proton_or_electron(x_n)
x_e = x_p
energy_densities = []
for x, particle in zip([x_n, x_p, x_e], ['n', 'p', 'e']):
energy_densities.append(energy_density(x, particle, include_constant))
return sum(energy_densities)
def total_pressure(x_n, include_constant=True):
x_p = fermi_momentum_of_proton_or_electron(x_n)
x_e = x_p
pressures = []
for x, particle in zip([x_n, x_p, x_e], ['n', 'p', 'e']):
pressures.append(pressure(x, particle, include_constant))
return sum(pressures)
Recall the "power law" $$ \overline{\varepsilon}(\overline{p}) = A_\mathrm{nonrel}\overline{p}^{3/5} + A_\mathrm{rel}\overline{p} \tag{44} $$
def power_law(p, A_nonrel, A_rel):
return A_nonrel * p**(3/5) + A_rel * p
# fit on dimensionless part only to avoid numerical issues of big exponents (>1e30)
x_plot = 10**np.linspace(-2, 2, 50)
p = np.array([total_pressure(x, include_constant=False) for x in x_plot])
eps = np.array([total_energy_density(x, include_constant=False) for x in x_plot])
param, _ = curve_fit(lambda p, A_nr, A_r: np.log10(power_law(p, A_nr, A_r)), p, np.log10(eps))
fig, ax = plt.subplots(dpi=160, figsize=(3.5, 3.0))
eps0 = epsilon_0
ax.errorbar(eps0*p, eps0*eps, fmt='o', color='black', ms=2, label='Numerical')
ax.errorbar(total_pressure(1.0).to('erg/cm^3').value,
total_energy_density(1.0).to('erg/cm^3').value,
fmt='x', ms=5, color='green', label=r'$x_n = 1$')
p_plot = np.linspace(1e-10, 1e7, 100)
ax.plot(eps0*p_plot, eps0*power_law(p_plot, A_nonrel=2.572, A_rel=2.891), color='blue', lw=1.0, label='Eq. (64)')
ax.plot(eps0*p_plot, eps0*power_law(p_plot, *param), color='red', ls='dashed', label='My fit')
ax.set_title(r'$x_n = 0.01$ to $x_n = 100$')
ax.set_xlabel(r'$p_\mathrm{tot}/\varepsilon_0$ (erg/cm$^3$)')
ax.set_ylabel(r'$\varepsilon_\mathrm{tot}/\varepsilon_0$ (erg/cm$^3$)')
ax.set_xscale('log')
ax.set_yscale('log')
ax.legend()
plt.show()
print(f'> My fit: {param}')
Result: $$ A_\mathrm{nonrel} = 2.572 \hspace{30pt} A_\mathrm{rel} = 2.891 \tag{64} $$ Compare with the pure neutron assumption in previous section (last week): $$ A_\mathrm{nonrel} = 2.4216 \hspace{30pt} A_\mathrm{rel} = 2.8663 \tag{46} $$
Some crust-related discussion authors made in the remarks:
VI. Introducing nuclear interactions¶
Two quantities we want to respect:
- Saturation density $n_0 = 0.16~/\mathrm{fm}^3$
- Binding energy per nucleon $\dfrac{\mathrm{BE}}{A} \approx -16~\mathrm{MeV}$
And two quantities we only know their rough constraints:
- Nuclear compressibility $K_0$ around $200$ ot $400~\mathrm{MeV}$
- Symmetry energy term ($Z = 0$) around $30~\mathrm{MeV}$ at $n_0$
Recall for fermions of spin $1/2$, $$ n = \frac{2}{(2\pi\hbar)^3}\int_0^{k_F} d^3k = \frac{k_F^3}{3\pi^2\hbar^3} \tag{7} \ , $$ which implies $$ k_F = \hbar\left(3\pi^2 n\right)^{1/3} $$
Hence the Fermi momentum of protons and neutrons (assuming symmetric nuclear matter, $N=Z$) is $$ k_F = \hbar\left(3\pi^2\cdot\frac{n_0}{2}\right)^{1/3} \approx 263~\mathrm{MeV}/c $$
n_0 = 0.16 / u.fm**3
k_F = const.hbar * (3 * np.pi**2 * 0.5 * n_0)**(1/3)
(k_F * const.c).to('MeV')
A. Symmetric nuclear matter¶
Goal: Relate $n_0$, $\mathrm{BE}$ and compressibility $K_0$ to energy density $\varepsilon(n)$. Symmetry energy term is reserved for next subsection.
Then we model the energy density per nucleon by $$ \frac{\varepsilon(n)}{n} = m_N + \frac{3}{5}\frac{k_F^2}{2m_N} + \frac{A}{2}u + \frac{B}{\sigma+1}u^\sigma \ . \tag{69} $$
- 1st term: rest mass energy
- 2nd term: kinetic energy; at $n_0$, it evaluates to $\langle E_F^0\rangle = 22.1~\mathrm{MeV}$, could be written as $\langle E_F^0\rangle u^{2/3}$
- Dimensionless $u\equiv n/n_0$
EF0 = ((3/5) * (k_F**2/(2 * 0.5*(const.m_p + const.m_n)))).to('MeV')
EF0
Result: $$ \sigma = \frac{K_0 + 2\langle E_F^0\rangle}{3\langle E_F^0\rangle - 9\mathrm{BE}} \tag{73} $$ $$ B = \frac{\sigma+1}{\sigma-1}\left[\frac{1}{3}\langle E_F^0\rangle - \mathrm{BE}\right] \tag{74} $$ $$ A = \mathrm{B.E.} - \frac{5}{3}\langle E_F^0\rangle - B \tag{75} $$
class PrakashEos:
def __init__(self, BE=-16.0, EF0=EF0.value, K0=400.0):
kwargs = locals(); kwargs.pop('self')
self.set_param(**kwargs)
def energy_density(self, u, include_rest_mass=True):
res = self.EF0 * u**(2/3)
res += 0.5 * self.A * u
res += (self.B * u**self.sigma) / (self.sigma + 1)
if include_rest_mass:
m_N = (0.5 * const.c**2 * (const.m_n + const.m_p)).to('MeV').value
res += m_N
return res
def pressure(self, u):
res = (2/3) * self.EF0 * u**(5/3)
res += 0.5 * self.A * u**2
res += self.B * self.sigma * u**(self.sigma+1) / (self.sigma + 1)
return 0.16 * res
def _sigma(self):
res = self.K0 + 2 * self.EF0
self.sigma = res / (3 * self.EF0 - 9 * self.BE)
def _B(self):
res = (self.sigma + 1) / (self.sigma - 1)
self.B = res * ((1/3) * self.EF0 - self.BE)
def _A(self):
self.A = self.BE - (5/3) * self.EF0 - self.B
def set_param(self, **kwargs):
self.__dict__.update(kwargs)
self._sigma()
self._B()
self._A()
eos = PrakashEos()
eos.A, eos.B, eos.sigma
At $K_0 = 400~\mathrm{MeV}$, the authors report a result of $$ A = -122.2~\mathrm{MeV} \hspace{30pt} B = 65.39~\mathrm{MeV} \hspace{30pt} \sigma = 2.112 $$
fig, ax = plt.subplots(ncols=2, dpi=160, figsize=(8, 3), constrained_layout=True)
u_plot = np.linspace(0, 2, 1000)
for K0, color in zip([200, 300, 400], ['red', 'green', 'blue']):
eos.set_param(K0=K0)
ax[0].set_title('Reproducing Fig. 7')
ax[0].set_ylabel(r'$E/A - m_N$ ($\mathrm{MeV}$)')
ax[0].plot(u_plot, eos.energy_density(u_plot, include_rest_mass=False),
color=color, label=f'{K0} MeV')
ax[1].set_title('Reproducing Fig. 8')
ax[1].set_ylabel(r'$p$ ($\mathrm{MeV}/\mathrm{fm}^3$)')
ax[1].plot(u_plot, eos.pressure(u_plot), color=color, label=f'{K0} MeV')
for axe in ax:
axe.axhline(0.0, color='black', lw=1.0)
axe.set_xlabel(r'$u \equiv n/n_0$')
axe.set_xlim(0, )
axe.legend(title=r'$K_0$')
plt.show()
K0_plot = np.linspace(200, 400, 100)
data = []
for K0 in K0_plot:
eos.set_param(K0=K0)
data.append([K0, eos.A, eos.B, eos.sigma])
df = pd.DataFrame(data, columns=['K0', 'A', 'B', 'sigma'])
fig, ax = plt.subplots(ncols=2, dpi=160, figsize=(8, 3), constrained_layout=True)
c = 0
ax[c].plot(df['K0'], df['A'], color='red', label='A')
ax[c].plot(df['K0'], df['B'], color='green', label='B')
c = 1
ax[c].plot(df['K0'], df['sigma'], color='blue', label=r'$\sigma$')
for axe in ax:
axe.set_xlabel(r'$K_0$ ($\mathrm{MeV}$)')
axe.legend()
plt.show()
B. Nonsymmetric nuclear matter¶
Introduce asymmetry $\alpha$ such that $$ n_n = \frac{1+\alpha}{2}n \hspace{30pt} n_p = \frac{1-\alpha}{2}n \tag{78} $$ Indeed, $$ \alpha = \frac{n_n - n_p}{n} \tag{79} $$
Contributions to the kinetic part of $\varepsilon$: $$\begin{aligned} \varepsilon_\mathrm{KE}(n, \alpha) &= \frac{3}{5}\frac{k_{F, n}^2}{2m_N}n_n + \frac{3}{5}\frac{k_{F, p}^2}{2m_N}n_p \\ &= n\langle E_F\rangle\cdot\frac{1}{2}[(1+\alpha)^{5/3} + (1-\alpha)^{5/3}] \end{aligned} \tag{81} $$ where $$ \langle E_F\rangle = \frac{3}{5}\frac{\hbar^2}{2m_N}\left(\frac{3\pi^2 n}{2}\right)^{2/3} $$ is the mean kinetic energy of symmetric matter at density $n$.
So the "excess kinetic energy" due to asymmetry of nuclear matter is $$\begin{aligned} \Delta\varepsilon_\mathrm{KE}(n, \alpha) &= \varepsilon_\mathrm{KE}(n, \alpha) - \varepsilon_\mathrm{KE}(n, 0) \\ &= n\langle E_F\rangle\cdot \left\{\frac{1}{2}[(1+\alpha)^{5/3} + (1-\alpha)^{5/3}] - 1\right\} \end{aligned} \tag{83}$$ For pure neutron matter, $\alpha = 1$, and $$ \Delta\varepsilon_\mathrm{KE}(n, \alpha) = n\langle E_F\rangle(2^{2/3} - 1) \tag{84} $$
Potential energy: $$ E(n, \alpha) = E(n, 0) + \alpha^2S(n) \tag{86} $$ where the authors assume $$ S(u) = (2^{2/3} - 1) \langle E_f^0\rangle (u^{2/3} - F(u)) + S_0F(u) \ . \tag{87}$$ Here, $F(u)$ is some functional on $u$ that we have to guess, and $S_0 = 30~\mathrm{MeV}$ is just the symmetry energy at saturation density.
def symmetry_energy(u, F=lambda u: u, S0=30.0):
return (2**(2/3)-1) * EF0.to('MeV').value * (u**(2/3) - F(u)) + S0*F(u)
u_plot = np.linspace(0, 2, 100)
eos.set_param(K0=400.0)
sym_ene = eos.energy_density(u_plot, include_rest_mass=False)
alpha = 1.0
asym_ene = alpha**2 * symmetry_energy(u_plot)
fig, ax = plt.subplots(dpi=160, figsize=(3.5, 3))
ax.plot(u_plot, sym_ene+asym_ene, color='blue', label=r'Fig. 9: $\alpha = 1$, $F(u)=u$')
ax.plot(u_plot, sym_ene, color='red', label=r'Fig. 7: $\alpha = 0$')
ax.set_title('Reproducing Fig. 9')
ax.set_xlabel(r'$u = n/n_0$')
ax.set_ylabel(r'$E(n, \alpha=1)/A-m_N$ ($\mathrm{MeV}$)')
ax.set_xlim(0, )
ax.axhline(0.0, color='black', lw=1.0)
ax.legend()
plt.show()
Similarly, the pressure is given by $$ p(u, \alpha) = p(u, 0) + n_0\alpha^2\left[\frac{1}{3}(2^{2/3}-1) \langle E_F^0\rangle (2u^{5/3}-3u^2) + S_0u^2\right] \tag{89}$$
def symmetry_pressure(u, alpha=1.0, S0=30.0):
res = (1/3) * (2**(2/3) - 1) * EF0.to('MeV').value * (2*u**(5/3) - 3*u**2) + S0 * u**2
return 0.16 * alpha**2 * res
u_plot = np.linspace(0, 8, 100)
alpha = 1.0
eos.set_param(K0=400.0)
sym_ene = eos.energy_density(u_plot, include_rest_mass=False)
asym_ene = alpha**2 * symmetry_energy(u_plot)
sym_pres = eos.pressure(u_plot)
asym_pres = symmetry_pressure(u_plot)
fig, ax = plt.subplots(dpi=160, figsize=(3.5, 3))
ax.plot(sym_ene+asym_ene, sym_pres+asym_pres, color='blue', label=r'Numerical')
ax.set_title('Reproducing Fig. 11')
ax.set_xlabel(r'$\varepsilon$ ($\mathrm{MeV}/\mathrm{fm}^3$)')
ax.set_ylabel(r'$p$ ($\mathrm{MeV}/\mathrm{fm}^3$)')
ax.set_xlim(0, )
ax.set_ylim(0, )
ax.legend()
plt.show()
Polytropic fit $$ p(\varepsilon) = \kappa_0 \varepsilon^\gamma$$ Found to be $$ \kappa = 3.548\times10^{-4} \hspace{30pt} \gamma = 2.1 $$
C. Does the speed of sound exceed that of light?¶
To recover causality, $$ V_\mathrm{nucl}(u, 0) = \frac{A}{2}u + \frac{B}{\sigma+1}\cdot\frac{u^\sigma}{1 + Cu^{\sigma-1}} \tag{94} $$
D. Pure neutron star with nuclear interaction¶
E. What about a cosmological constant?¶
