Application: Degenerate Electrons#
In white dwarfs, electron degeneracy pressure provides the dominant support against gravity. Electrons are fermions, and when the density of the star is high, they are packed closely together. Quantum mechanical effects come into play, and the electrons no longer behave as an ideal gas, but instead exhibit large pressure even at very low (or zero) temperature.
The origin of this is Fermi-Dirac statistics—no two electrons can occupy the same state, so as they get confined to a smaller and smaller volume, they need to take on larger momenta.
The distribution function describes the properties of the electrons, and has the following meaning: \(n(p) d^3 x d^3 p\) is the number of particles with momentum \(p\) in a volume \(d^3x\). Our distribution function is:
where \(\mathcal{E}\) is the kinetic energy of the electron, \(T\) is the temperature, and \(\Psi\) is the degeneracy parameter which is related to the chemical potential, \(\mu\), and the Fermi energy, \(\mathcal{E}_F\), as:
(Note: some sources use \(\eta\) instead of \(\Psi\) for the degeneracy parameter)
The electrons are degenerate when \(\Psi \gg 1\), which means that the Fermi energy is much greater than the thermal energy. The electrons behave as an ideal gas when \(-\Psi \gg 1\) (in this case, the \(+1\) in the denominator of \(n(p)\) is insignificant and we get the Maxwell-Boltzmann distribution).
To get the number density, we integrate over all momenta:
The pressure and specific internal energy, \(e\), are then found as:
Let’s look at the behavior of the distribution function. Consider:
let’s define
then we have
and for each choice of \(\Psi\) we have a unique distribution.
We note that \(\xi/\Psi = \mathcal{E}/\mathcal{E}_F\).
import numpy as np
import matplotlib.pyplot as plt
def F(xi, psi=100):
return 1.0 / (np.exp(xi - psi) + 1)
Let’s look at the case of \(\Psi > 1\)
fig, ax = plt.subplots()
# we'll work in terms of r = xi / Psi
r = np.linspace(0, 4, 200)
for psi in [1, 2, 10, 20, 100]:
ax.plot(r, F(psi * r, psi=psi), label=f"$\Psi = {psi}$")
ax.legend()
ax.set_ylabel("$F$")
ax.set_xlabel(r"$\mathcal{E}/\mathcal{E}_F$")
Text(0.5, 0, '$\\mathcal{E}/\\mathcal{E}_F$')
As we see, as \(\Psi \rightarrow \infty\), the distribution, \(F(\mathcal{E})\) becomes a step function. This is the limit of complete degeneracy, and we have:
We can alternately use the Fermi momentum, \(p_F\) to describe the location of the step.
With this form, the integral for number density is trivial:
where we switched to spherical coordinates in momentum space: \(d^3p \rightarrow 4\pi p^2 dp\).
For other values of \(\Psi\), we have a problem:
We don’t know what the value of \(\Psi\) is ahead of time. Instead, we typically know what the density of the star is, and we can find the number density of electrons as:
where \(Y_e\) is the electron fraction (typically around \(1/2\) for compositions heavier than hydrogen). This means that given \(\rho, T\), we can solve for \(\Psi\).
Once we have \(\Psi\), we can then find the pressure and energy of the gas.
Finding the trend of \(\Psi\) with \(\rho, T\)#
Let’s now implement an algorithm to find the \(\Psi\) that corresponds to an input \(\rho\) (or \(n_e\)) and \(T\).
We first need a function to compute \(n_e\) given a \(T\) and a guess for \(\Psi\). We’ll use our composite Simpson’s integration for this.
We’ll assume that we are non-relativistic, so
Our integral is then:
We can make a change of variables:
and then we have:
Note: that for \(\rho > 10^6~\mathrm{g~cm^{-3}}\) relativistic effects are important, and we would need to use a more general form for \(\mathcal{E}(p)\).
Let’s write a function that returns the integrand
def integrand(x, psi):
return x**2 / (np.exp(x**2 - psi) + 1.0)
and now lets plot this for various \(\Psi\)
fig, ax = plt.subplots()
x = np.linspace(0, 100, 10000)
for psi in [-100, -10, -2, -1, 0, 1, 2, 10, 100]:
ax.plot(x, integrand(x, psi), label=f"$\Psi = {psi}$")
ax.legend()
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_ylim(1.e-2, 1.e7)
/tmp/ipykernel_2899/3209928939.py:2: RuntimeWarning: overflow encountered in exp
return x**2 / (np.exp(x**2 - psi) + 1.0)
(0.01, 10000000.0)
Now let’s try scaling the integrand into \([0, 1]\). We’ll use the same technique we did for integrating the Planck function.
SMALL = 1.e-30
def zv(x, c):
""" transform the variable x -> z """
return x/(c + x)
def xv(z, c):
""" transform back from z -> x """
return c*z/(1.0 - z + SMALL)
fig, ax = plt.subplots()
z = np.linspace(0, 1, 100)
c = 2.0
x = np.linspace(0, 100, 10000)
for psi in [-100, -10, -2, -1, 0, 1, 2, 10, 100]:
ax.plot(z, integrand(xv(z, c), psi), label=f"$\Psi = {psi}$")
ax.legend()
ax.set_yscale("log")
ax.set_ylim(1.e-10, 100)
/tmp/ipykernel_2899/3209928939.py:2: RuntimeWarning: overflow encountered in exp
return x**2 / (np.exp(x**2 - psi) + 1.0)
(1e-10, 100)
This looks reasonable. We might need to worry about resolving the sharp drop for the very large \(\Psi\) (when we are completely degenerate). But we’ll produce an estimate of the error in the integral as we compute it.
First the fundamental constants we need
# CGS constants
h_planck = 6.63e-27
k_B = 1.38e-16
m_e = 9.11e-28
m_u = 1.67e-24
Now a Simpson’s rule integrator for computing our integral, using the techniques we developed earlier for integrating to infinity.
def fd_integral(N, psi):
"""compute the integral over the Fermi-Dirac distribution
using Simpsons rule with N intervals"""
assert N % 2 == 0
# we will transform from integrating over x to integrating over
# z = x / (c + x) with z = [0, 1]
c = 2.0
z = np.linspace(0.0, 1.0, N+1, endpoint=True)
I = 0.0
for n in range(0, N, 2):
fl = integrand(xv(z[n], c), psi) / (1.0 - z[n] + SMALL)**2
fc = integrand(xv(z[n+1], c), psi) / (1.0 - z[n+1] + SMALL)**2
fr = integrand(xv(z[n+2], c), psi) / (1.0 - z[n+2] + SMALL)**2
I += (1.0/3.0) * (z[n+1] - z[n]) * (fl + 4*fc + fr)
I *= c
return I
Now a function that returns \(n_e\). This does the integral, ensuring we meet some accuracy, and adds in the physical constants to give us \(n_e\).
def compute_ne(T, psi):
"""given a temperature and degeneracy parameter, psi,
compute the number density of electrons by integrating
the Fermi-Dirac distribution over all momenta"""
# we'll do the dimensionless integral of
# x**2 / (exp(x**2 - psi) + 1) using Simpson's rule
# we will pick a value of N and do the integration
# and then change N until the error is small
# we'll start with N = 50 and then keep doubling
# until the difference in I is small
tol = 1.e-8
N = 50
err = 1.e30
I_old = 1.e30
while err > tol:
I = fd_integral(N, psi)
err = np.abs(I - I_old) / np.abs(I_old)
I_old = I
N *= 2
ne = 8.0 * np.pi / h_planck**3 * (2.0 * m_e * k_B * T)**1.5 * I
return ne
Now we want to do a root find on this. We will start with bisection, although we note that there are better methods we could do.
def find_psi(rho, T, Ye=0.5, tol=1.e-6):
"""given rho, T, we want to find the degeneracy parameter, psi"""
ne_input = rho * Ye / m_u
psi_low = -100
psi_high = 1000
# we want to zero ne(T, psi) - ne_input)
ne_low = compute_ne(T, psi_low) - ne_input
ne_high = compute_ne(T, psi_high) - ne_input
if ne_low * ne_high > 0:
print("no root in this interval")
return None
err = 1.e10
psi_mid = 0.5 * (psi_low * psi_high)
while err > tol:
ne_mid = compute_ne(T, psi_mid) - ne_input
if ne_low * ne_mid > 0:
# the root is in the right half of the interval
psi_low = psi_mid
ne_low = ne_mid
else:
psi_high = psi_mid
ne_high = ne_mid
psi_mid = 0.5 * (psi_low + psi_high)
err = np.abs(psi_low - psi_high) / np.abs(psi_mid)
return psi_mid
This is an implementation of the secant method which is essentially Newton’s method with the derivative computed via a finite difference.
def find_psi2(rho, T, Ye=0.5, tol=1.e-6):
"""given rho, T, we want to find the degeneracy parameter, psi"""
ne_input = rho * Ye / m_u
psi_m1 = 10
psi_0 = 20
# we want to zero ne(T, psi) - ne_input)
ne_m1 = compute_ne(T, psi_m1) - ne_input
err = 1.e10
while err > tol:
ne_0 = compute_ne(T, psi_0) - ne_input
dne_dpsi = (ne_0 - ne_m1) / (psi_0 - psi_m1)
psi_m1 = psi_0
ne_m1 = ne_0
psi_0 -= ne_0 / dne_dpsi
err = np.abs(ne_0 / dne_dpsi) / np.abs(psi_0)
return psi_0
T = 1.e7
rhos = np.logspace(-4, 6, 50)
psis = []
for rho in rhos:
psi = find_psi2(rho, T)
psis.append(psi)
/tmp/ipykernel_2899/3209928939.py:2: RuntimeWarning: overflow encountered in exp
return x**2 / (np.exp(x**2 - psi) + 1.0)
fig, ax = plt.subplots()
ax.plot(rhos, psis)
ax.set_xlabel(r"$\rho$")
ax.set_ylabel(r"$\Psi$")
ax.set_xscale("log")
Improvements#
Simpson’s rule is not the best method for integrating this. In particular, we see that for \(\Psi \gg 1\), we need a lot of points because of the sharp step in the distribution function. A better method would be to not use uniform intervals when doing this integral.
For example, scipy.integrate.quad uses the QUADPACK library to do the integration, and uses a Gaussian quadrature for constructing the integral.
The secant method begins with an initial guess, and we could try to seed the initial guess for \(\Psi\) with the guess from the previous \(\rho\) solution. This would accelerate convergence.
Physically, we should extend \(\mathcal{E}(p)\) to be relativistic, which would allow us to probe higher densities.
SciPy implementation#
Here’s a version using the integration and root finding functions available in SciPy
from scipy import integrate, optimize
def ne_scipy(T, psi):
I, err = integrate.quad(integrand, 0, np.inf, args=(psi))
ne = 8.0 * np.pi / h_planck**3 * (2.0 * m_e * k_B * T)**1.5 * I
return ne
def find_psi_scipy(rho, T, Ye=0.5):
ne_input = rho * Ye / m_u
psi_low = -100
psi_high = 1000
psi = optimize.brentq(lambda psi: ne_input - ne_scipy(T, psi), psi_low, psi_high)
return psi
T = 1.e7
rhos = np.logspace(-4, 6, 50)
psis = []
for rho in rhos:
psi = find_psi_scipy(rho, T)
psis.append(psi)
fig, ax = plt.subplots()
ax.plot(rhos, psis)
ax.set_xlabel(r"$\rho$")
ax.set_ylabel(r"$\Psi$")
ax.set_xscale("log")
/tmp/ipykernel_2899/3209928939.py:2: RuntimeWarning: overflow encountered in exp
return x**2 / (np.exp(x**2 - psi) + 1.0)