Table of Contents

• 1. Background

  • 1.1 Units

  • 1.2 Thermodynamics

• 2. Homogenous Nuclear Matter

  • 2.1 Electric Neutrality

  • 2.2 Chemical (Beta) Equilibrium

    • 2.2.1 Energy "Gap"

    • 2.2.2 Equilibrium with Nuclear Matter

    • 2.2.3 Energy-per-Particle and Chemical Potentials

    • 2.2.4 Pressure

• 3. Neutron Star Equation of State

• 4. Saturation

  • 4.1 Thermodynamic Equilibrium

  • 4.2 Coulomb Interaction

• 5. Surface Energy

• 6. Liquid Drop Model (LDM)

• 7. Problems

  • 7.1 Collapse of the Earth

1. Background

1.1 Units

In the following, we shall generally adopt "natural" units in which we choose our physical scales of length, time, temperature, mass, and charge to match the system under consideration making some of the dimensionful physical constants assume unit value. In particular:

• $k_B = 1$: The Boltzmann constant $k_B$ simply serves to convert energies to temperatures $E=k_BT$. By setting this to unity, we simply choose to measure temperatures in energy units (i.e. eV, MeV, or GeV in this context). Thus, if you see a "temperature" of 1MeV, simply divide by $k_B$ to obtain $10^{10}$K.

• $c = 1$: Since we will be dealing with relativistic energy scales, we can express all masses in terms of energy using Einstein's famous relationship $E=mc^2$. For nuclear physics this is extremely convenient. For example, an electron has a mass of 511keV. This means that if an electron annihilates with a positron producing two photons, each photon will also have an energy of 511keV conserving the total mass-energy of the system. Nucleons (protons and neutrons) are about 1GeV ($m_p \approx 0.938$GeV and $m_n \approx 0.940$GeV. Note that the neutron is heavier allowing it to decay into a proton, electron and anti-neutrino. Were the proton more massive, then they could decay into neutrons and all the earth would collapse... see below.)

• $\hbar = 1$: This is another common convention that allows one to relate energy and frequency, or momentum and distance.

• $\epsilon_0=1$: This is often used for convenience when discussing electromagnetism.

In addition, one sometimes specifies the mass scale:

• $m=1$: In some cases - such as discussing a single species - it is useful to use a natural mass unit which is the mass of the principle species. Here, however, we shall not adopt this convention since we discuss particles of several different masses.

Note: Since there are only 5 fundamental dimensions, one may only specify 5 dimensionful constants such as discussed above. These are not the only choices - when discussing the TOV equations for the structure of neutron stars for example, one might like to set $G=1$ (replacing $\hbar=1$) and express masses in terms of solar masses $M_\circ = 1$. In this case, one must release some of the previous constraints.

On the other hand, dimensionless constants such as $2\pi$ or $\alpha$ can never be adjusted unless they can be absorbed into a unit redefinition. For our purposes, the relevant dimensionless quantity is the fine-structure constant:

1.2 Thermodynamics

Here we present a brief review of the required thermodynamics for this problem.

Our starting point is the Helmholtz free energy $F = E - TS$. For a given quantum system with fixed boundary conditions (i.e. volume) and number of particles, this can be calculated from the spectrum (energy eigenvalues) of the system:

where $\beta = k_B/T$ and $k_B$ is the Boltzmann constant and $T$ is the temperature. In the last step we have expressed this directly in terms of the Hamiltonian ParseError: KaTeX parse error: Undefined control sequence: \op at position 1: \̲o̲p̲{H} of the system. For an explanation of this result, see Schrödinger's book Statistical Thermodynamics. Briefly, it comes from considering an arbitrary distribution (mixed state) over a complete set of states $\{\ket{n}\}$ which can be expressed in terms of a density matrix ParseError: KaTeX parse error: Undefined control sequence: \mat at position 1: \̲m̲a̲t̲{R}:

Formally, the equilibrium ensemble can be found by choosing ParseError: KaTeX parse error: Undefined control sequence: \mat at position 1: \̲m̲a̲t̲{R} to maximum the number of configurations with a given energy ParseError: KaTeX parse error: Undefined control sequence: \op at position 13: E = \braket{\̲o̲p̲{H}}_{\op{R}} =… -- this corresponds to the most probable configuration as discussed in Schrödinger's book. The entropy provides a measure of the (log of the) number of configurations:

To implement this constrained extremization, we introduce the energy constraint with a Lagrange multiplier $\beta$ which will later be interpreted as the inverse temperature $\beta = 1/k_B T$.

Performing the variation with respect to ParseError: KaTeX parse error: Undefined control sequence: \op at position 1: \̲o̲p̲{R}^T we obtain the following solution, where we have eliminated $\lambda$ to satisfy the normalization constraint, but left the Lagrange multiplier $\beta$ in place of the energy:

The quantity ParseError: KaTeX parse error: Undefined control sequence: \op at position 1: \̲o̲p̲{Z} is referred to as the "partition function" and the Helmholtz free energy is expressed in terms of this:

Here we will work at $T=0$, in which case $F=E$ and we will follow the convention of working with the energy $E$ rather than the free energy $F$.

This should be expanded with some intuitive discussions added, perhaps some plots of the Legendre transform etc.

For extensive systems, it is convenient to work with two intrinsic quantities: the energy density $\mathcal{E}$ and energy per particle $\varepsilon$. From these we may express the chemical potential $\mu$ and pressure $P$:

2. Homogenous Nuclear Matter

2.1 Electric Neutrality

We start by discussing homogeneous nuclear matter comprising neutrons and protons, along with leptons – electrons and possibly muons – required to maintain electrical neutrality. Bulk matter must overall be electrically neutral otherwise the energy of the system will not be extensive and will be huge.

We start by modeling the leptons as a free Fermi gas filling all states up the the Fermi surface which has Fermi momentum $p_F = \hbar k_F$. Here we consider a single component (for example spin-up electrons). The densities and chemical potentials involved will be sufficiently large that we must consider relativistic effects. The energy of a (plane wave) state with momentum $p=\hbar k$ is:

The density of the Fermi gas is thus obtained by counting all states in the Fermi sphere and the energy density is obtained by integrating $E(k)$ over the same region:

In neutron stars, the lepton chemical potential $\mu_e \sim \hbar c k_F \approx 10$ – $100$MeV $\gg m_ec^2 = 511$keV. Thus, as a first approximation, we can take the ultra-relativistic limit and assume that $0 \approx mc^2 \ll \mu_l$. This simplifies the integrals and we obtain:

For comparison, in the non-relativistic limit (such as in solid-state physics), one has $mc^2 \ll \mu_e$ so that ParseError: KaTeX parse error: Undefined control sequence: \order at position 64: …hbar^2k^2/2m + \̲o̲r̲d̲e̲r̲(\mu/mc^2)^2 so that

In natural units, $m=c=\hbar=1$, the full integral is:

Exercise: Use dimensional analysis to argue that the following result is correct including the omitted factors.

where in the last equation, we have assumed that $mc^2 \ll \hbar^2k_F^2$ so that the mass term can be neglected. In neutron stars, $E(k_F) = \mu_e \sim 10-100$MeV which is large compared with the mass of the electron $m_ec^2 \approx 511$keV$\approx 0.5$MeV, but not compared with the mass of the lightest muon $m_\mu c^2 \approx 100$MeV.

Here we check these results numerically:

Electrons and muons come with two spins, so to each of these equations we must add a degeneracy factor $g=2$ so that the final result is:

Solving these equations gives the energy-density of the (non-interacting) electron gas $\mathcal{E}_e(n_p)$ and muon gas $\mathcal{E}_\mu(n_\mu)$. Explicitly, for the electrons, when the approximations can be made, we have:

Note that the thermodynamic relationships are satisfied (by construction).

2.2 Chemical (Beta) Equilibrium

Physically, charge neutrality is maintained by $\beta$ decay, inverse $\beta$ decau, and muon decay processes:

Note that each reaction conserves electric charge, baryon charge, and lepton number.

Charge neutrality demands that the lepton number density equal the proton density:

Now consider the total lepton density $n_l$ as being fixed by nuclear physics (we will relax this in a moment). What electron fraction $x = n_e/n_l$ will minimize the energy density of the lepton gas:

Minimizing with respect to the electron fraction gives the equilibrium condition:

I.e. chemical equilibrium is achieved (via the muon-decay process) when the chemical potentials of the two species are equal. For a fixed lepton chemical potential $\mu_l = \mu_e = \mu_\mu$, we thus have the following total lepton density:

where the first form contains only electrons.

2.2.1 Energy "Gap"

What exactly does the previous expression mean? If $\mu_e < m_\mu c^2$, then there are no muons in the system, but the muonic equation of state does not appear to have a solution. What happens here is that the muonic spectrum has a "gap" (or more appropriately a "mass gap" in this case). For $-m_\mu c^2 < \mu < m_\mu c^2$ the muonic system stays in the same state – namely the vaccuum state $n_\mu = 0$. This is denoted by kink in the equation of state:

The chemical potential $\mu = \mathcal{E}'(n)$ is the slope here. The kink at $n=0$ implies that there is a range of chemical potentials over which the state does not change. This is the manifestation of a gap in the thermodynamic functions. To actually see the "gap" look at the single-particle energy spectrum $\mu = E(k)$:

The origin of this kink and gap is the zero-temperature limit $T=0$. For any finite temperature, one must be more careful and populate the various states with the appropriate Fermi distribution function that includes both particles and anti-particles:

These results can be expressed in terms of the Fermi-Dirac integrals and indeed smooth out the equation of state replacing the kink and gap with smooth functions over the energy scale $1/\beta = k_B T$:

As the temperature gets smaller and smaller, the distribution becomes sharper and sharper until at $T=0$, one has a step function at $k=k_F$. If ParseError: KaTeX parse error: Undefined control sequence: \abs at position 1: \̲a̲b̲s̲{\mu} \leq mc^2, then the distribution vanishes, but for finite $T$, there is a finite occupation for any positive $\mu$.

2.2.2 Equilibrium with Nuclear Matter

The total energy density of a gas of protons, neutrons, electrons, and muons is:

where $\mathcal{E}_n(n_n, n_p)$ is the nuclear equation of state, and we have used the condition of chemical equilibrium to determine the optimal ratio of electrons to muons as discussed in the previous section.

In general, one should consider as fixed quantities which are conserved. In this case, total baryon number (density) $n_B$ is a conserved quantity (at least in terms of the $\beta$-decay reactions), so it makes sense to describe the system in terms of $n_B$ and the proton fraction $x_p = n_p/n_B$ so that $n_p = x_pn_B$, and $n_n = (1-x_p)n_B$. Electric neutrality still requires $n_p = n_l$, so minimizing the energy density with respect to the proton fraction gives the $\beta$-equilibrium condition:

$$ \mu_p - \mu_n = \mu_l = \mu_e = \mu_\mu, \qquad \mu_{n,p} = \pdiff{\mathcal{E}_n(n_n, n_p)}{n_{n,p}}. \tag{$\betaParseError: KaTeX parse error: Expected 'EOF', got '}' at position 13: -equilibrium}̲ $

Again – the appropriate chemical equilibrium conditions apply.

Let's now take stock of what we have. First we need to fix the overall baryon density $n_B$: this is fixed by the external potential – gravity in the case of a neutron star – which sets the corresponding chemical potential $\mu_B = \partial\mathcal{E}/\partial n_B$. For this fixed density, one minimizes the ratio of the various other components satisfying $n_B = n_p + n_n$, and the neutrality condition $n_p = n_l = n_e + n_\mu$. These minimizations give the chemical equilibrium conditions that $\mu_p - \mu_n = \mu_l = \mu_e = \mu_\mu$. Now, the chemical potential difference $\mu_p - \mu_n = \partial \mathcal{E}_n(n_B, x_p)/\partial x_p$ is a function of $n_B$ and the proton fraction. Similarly, the lepton chemical potential $\mu_l = \mathcal{E}_{l}'(n_l)$ is a function of $n_l = n_p = x_pn_B$. Hence, the $\beta$-equilibrium condition provides an equation that determines the proton fraction at a fixed baryon densith:

Since we know the form of $\mu_e = \mu_l(n_p)$, we can express this as

The actual form of the last equation depends on the form of the nuclear equation of state, which we will discuss below.

2.2.3 Energy-per-Particle and Chemical Potentials

It is common to express the nuclear equation of state in terms of the energy-per-particle or energy-per-baryon $\varepsilon_n(n_B, x_p) = \mathcal{E}(n_n, n_p)/n_B$ since this is a minimum at saturation as we discuss below. Here are the various thermodynamic relationships for nuclear matter:

In equilibrium, the total energy density of homogeneous matter can be expressed as $\mathcal{E} = \mathcal{E}_n(n_n, n_p) + \mathcal{E}_l(n_p)$. In $\beta$ equilibrium, this will be minimized with respect to the proton fraction $\partial\mathcal{E}/\partial{x_p} = 0$ which allows us to define the baryon chemical potential $\mu_B = \partial\mathcal{E}/\partial{n_B}$:

The last relationship shows that, in $\beta$-equilibrium, $\mu_B = \mu_n$. The proton chemical potential differes $\mu_p = \mu_B - \mu_l \neq \mu_B$ due to the Coulomb interaction.

2.2.4 Pressure

We now address one final point: the pressure of this homogeneous state $P$.

Note that although the full formula includes $\mu_e n_e + \mu_\mu n_mu - \mathcal{E}_e(n_e) - \mathcal{E}_\mu(n_\mu)$ but chemical equilibrium requires $\mu_e = \mu_\mu = \mu_l$ so this can be expressed in terms of the single combination $\mu_l n_l - \mathcal{E}_{l}(n_l)$.

The form of this can be deduced from the energy per baryon as follows:

The last relationship holds because $\mu_n = \mu_B$ is the baryon chemical potential.

3. Neutron Star Equation of State

This notebook describes a simplified method of calculating the equation of state in the crust of a neutron star. It uses the so-called compressible liquid drop model (CLDM) to decribe nuclei in the crust and makes a few approximations that should be expanded in a realistic calculation.

1. We assume that the nuclei have equal numbers of protons and neutrons. As such, our nuclear equation of state depends only on the total "baryon" density $n_B = n_p + n_n$. In reality, there are more neutrons than protons and this is described by an additional dependence on the proton fraction $x=n_p/n_B$ which we fix here to $x=1/2$. In what follows, we thus simply refer to a single density $n = n_B = 2n_n = 2n_p$.

4. Saturation

The starting point for a thermodynamic exploration is to use the Hemlholtz free-energy $F = E - TS$. We will be considering low temperatures (on the nuclear scale) and so consider the limit $T=0$ in which case we can work with the energy. Thus, we shall talk about minimizing the energy when one should really generalize this to the free energy $F$ at finite temperature.

We start with the notion of "saturation". The nuclear equation of state has a minimum energy at a finite density. This means that, thermodynamically, the density inside a nucleus tends to a constant finite value called the "saturation density" $n_0 \approx 0.16$fm$^{-3}$. Thus, if the volume of the nucleus is $V$, then the atomic number $A = nV$. In the (non-compressible) liquid drop model (LDM) $n=n_0$ is constant and fixed and the energy density of nuclear matter at this density is expressed through the constant $\epsilon_0 \approx -16$MeV which is the binding energy per nucleon $\epsilon_0 = E(A)/A = \mathcal{E}_n(n_0)/n_0$. Thus, in the LDM, the binding energy $E_n(A)$ and energy density $\mathcal{E}_n/V$ of a nucleus has the form:

In the CLDM, one must specify the energy density $\mathcal{E}_n(n)$ as a function of arbitrary $n$ to allow for the density to change. We must thus introduce an additional parameter $K \approx 250$MeV called the compressibility. The conventional parametrization of the nuclear equation of state is expressed in terms of the energy per nucleon:

where $\delta$ expresses the deviation from the saturation density.

4.1 Thermodynamic Equilibrium

Now consider a nucleus of volume $V$ in a box (cell) of volume $V_c > V$. Let $w=V/V_c$ be the "volume fraction". To determine the thermodynamic equlibirum, we can minimize the energy as a function of particle number $A = Vn$ at fixed total volume $V_c$ with respect to the size of the nucleus expressed in terms of $w$:

where we have introduced $\mu_n = \mathcal{E}'(n)$ as the baryon chemical potential in the nucleus. Note that we have introduced a factor of $V_c^{-1}$ to show that this implies zero pressure $P_n(n) = \mu_n n - \mathcal{E}_n(n)$ in the nucleus is zero, matching the zero-pressure of the surrounding empty space (vaccuum).

Note that this implies that the energy per nucleon is minimized $\partial(\mathcal{E}_n/n)\partial n = 0$ in a saturating nucleus, not the energy density $\mathcal{E}_n(n)$. This is a consequence of the Maxwell construction of phase equilibrium between the vacuum and the nucleus. However, the total energy density of the box is minimized since this is proportional to the energy per baryon:

Moral: in a fixed total volume, the energy density of the total volume is minimized, but not the energy-density of a given phase when there may be phase-coexistance. The latter must be determined by allowing material to redistribute to satisfy chemical (we shall discuss this below) and mechanical equilibrium in terms of equal chemical potentials $\mu$ and pressures $P$.

We now do this numerically and make some plots:

Notice that at saturation density, the energy per baryon (right) is minimized, not the energy density. Now we perform the minimization numerically:

Here we see how the volume fraction of the nucleus increases as a function of the total baryon density in the box. Numerically we will have problems at the boundaries unless we give some constraints to the minimization since $0\leq w \leq 1$. Here is how we can provide boundaries:

Another approach is to enforce the constraint via Lagrange multipliers.

Here we see that for sufficiently negative chemical potentials, either the volume fraction is zero or the density is zero, meaning that the system has no nucleus: the negative chemical potential has esentially excluded all nucleons. (There is some strange numerical behaviour due to the fact that for negative chemical potentials, the system is singular. This should straighten itself out once we allow the system more freedom to have a uniuque solution).

4.2 Coulomb Interaction

The next step is to include the Coulomb interactions. The physical situation is that of a lattice of nuclei, so the full calculation needs to account for the structure of the lattice. However, once can obtain a reasonable approximation (typically accurate to a few percent) using a "Wigner-Seitz cell". This amounts to placing the nucleus at the center of a spherical shell with the same volume $V_c = \tfrac{4\pi}{3}r_c^3$ as the corresponding unit cell of the true lattice.

Here is the energy of a solid uniformly charged sphere of charge $Z$ and radius $r$: