Chapter 19

Nuclear Astrophysics

The Gamow Peak & Thermonuclear Reaction Rates

In stellar interiors, nuclear reactions occur between thermally distributed nuclei that must tunnel through the Coulomb barrier. The reaction rate per pair of particles is:

$$\langle\sigma v\rangle = \left(\frac{8}{\pi\mu}\right)^{1/2} \frac{1}{(kT)^{3/2}} \int_0^\infty S(E)\,\exp\left(-\frac{E}{kT} - \frac{b}{\sqrt{E}}\right) dE$$

where $S(E)$ is the astrophysical S-factor (which varies slowly with energy), and$b = \pi\alpha Z_1 Z_2\sqrt{2\mu c^2}$ parameterizes the Coulomb barrier. The Gamow energy characterizes the Coulomb barrier height:

$$E_G = \left(\pi\alpha Z_1 Z_2\right)^2 \cdot 2\mu c^2 = \frac{b^2}{4}$$

The integrand is a product of the Maxwell-Boltzmann distribution (falling exponential) and the tunneling probability (rising with energy). Their product peaks at the Gamow peak energy:

$$E_0 = \left(\frac{E_G\,(kT)^2}{4}\right)^{1/3} \qquad \text{(Gamow peak energy)}$$

The width of the Gamow peak is $\Delta = \frac{4}{\sqrt{3}}\sqrt{E_0 kT}$. For the p+p reaction at $T = 15$ MK: $E_0 \approx 6$ keV and $\Delta \approx 7$ keV. This is far below the Coulomb barrier (~500 keV), showing that quantum tunneling is essential for stellar nucleosynthesis.

The pp Chain

The proton-proton chain is the dominant energy source in stars with $M \lesssim 1.3\,M_\odot$(including our Sun). The net reaction is:

$$4p \to {^4\text{He}} + 2e^+ + 2\nu_e + 26.73\;\text{MeV}$$

The chain begins with the rate-limiting step (a weak interaction!):

$$p + p \to d + e^+ + \nu_e + 0.42\;\text{MeV} \qquad (\tau \sim 10^{10}\;\text{yr})$$

followed by rapid deuterium burning:

$$d + p \to {^3\text{He}} + \gamma + 5.49\;\text{MeV} \qquad (\tau \sim 1\;\text{s})$$

The chain then branches: pp-I ($^3$He+$^3$He, 86%), pp-II (via $^7$Be, 14%), and pp-III (via $^8$B, 0.02%). The pp-III branch produces high-energy neutrinos detected by SNO and Super-Kamiokande. The energy generation rate scales as$\epsilon \propto T^4$.

The CNO Cycle

In stars more massive than ~1.3 $M_\odot$, the CNO cycle dominates hydrogen burning. Carbon, nitrogen, and oxygen serve as catalysts:

$^{12}\text{C} + p \to {^{13}\text{N}} + \gamma$

$^{13}\text{N} \to {^{13}\text{C}} + e^+ + \nu_e$

$^{13}\text{C} + p \to {^{14}\text{N}} + \gamma$

$^{14}\text{N} + p \to {^{15}\text{O}} + \gamma$ (slowest step)

$^{15}\text{O} \to {^{15}\text{N}} + e^+ + \nu_e$

$^{15}\text{N} + p \to {^{12}\text{C}} + {^4\text{He}}$

The net result is identical to the pp chain: $4p \to {^4\text{He}} + 2e^+ + 2\nu_e$. The CNO rate scales very steeply with temperature ($\epsilon \propto T^{16}$), making it extremely sensitive to the core temperature. The bottleneck is the$^{14}\text{N}(p,\gamma)^{15}\text{O}$ reaction, which determines the overall CNO rate.

Triple-Alpha Process & Advanced Burning

After hydrogen exhaustion, helium burning proceeds via the triple-alpha process at$T \sim 10^8$ K. This is a two-step process:

$$^4\text{He} + {^4\text{He}} \rightleftharpoons {^8\text{Be}} \qquad (T_{1/2} = 6.7 \times 10^{-17}\;\text{s})$$

$$^8\text{Be} + {^4\text{He}} \to {^{12}\text{C}}^* \to {^{12}\text{C}} + \gamma + 7.37\;\text{MeV}$$

This process is only possible because of the Hoyle state: an excited $0^+$ state in $^{12}$C at 7.654 MeV, just above the $^8$Be+$\alpha$ threshold (7.367 MeV). Fred Hoyle predicted this resonance from the observed cosmic abundance of carbon — a remarkable example of the anthropic principle in action.

Subsequent burning stages in massive stars proceed through C, Ne, O, and Si burning, each at progressively higher temperatures and shorter timescales, building up an onion-shell structure culminating in an iron core.

s-process and r-process Nucleosynthesis

Elements heavier than iron cannot be produced by charged-particle fusion (endothermic). Instead, they are built up by neutron capture:

$$(Z,A) + n \to (Z,A+1) + \gamma$$

s-process (slow neutron capture)

Neutron capture is slow compared to beta decay ($\tau_n \gg \tau_\beta$). The path follows the valley of stability. Occurs in AGB stars with neutron densities$\sim 10^7$ cm$^{-3}$. Produces about half of the elements heavier than iron, up to $^{209}$Bi. At steady-state, the s-process obeys:

$$\sigma_A \cdot N_A = \text{const} \qquad \text{(local approximation)}$$

r-process (rapid neutron capture)

Neutron capture is rapid compared to beta decay ($\tau_n \ll \tau_\beta$). The path runs along very neutron-rich nuclei far from stability. Requires extreme neutron densities $\sim 10^{24}$ cm$^{-3}$. Occurs in neutron star mergers (confirmed by GW170817/kilonova observation) and possibly in core-collapse supernovae. Produces the other half of heavy elements, including all thorium, uranium, and the r-process peaks at $A \approx 80, 130, 195$ (corresponding to neutron magic numbers N = 50, 82, 126 at the neutron drip line).

Neutron Stars & Nuclear Pasta

Neutron stars are the densest observable objects, with central densities reaching$5\text{--}10\,\rho_0$ (where $\rho_0 = 2.8 \times 10^{14}$ g/cm$^3$ is nuclear saturation density). The equation of state (EoS) relates pressure to density:

$$P = P(\rho, T, Y_p) \qquad \text{(EoS of dense nuclear matter)}$$

The TOV (Tolman-Oppenheimer-Volkoff) equation determines the structure:

$$\frac{dP}{dr} = -\frac{G\,\rho\,m}{r^2}\left(1 + \frac{P}{\rho c^2}\right)\left(1 + \frac{4\pi r^3 P}{mc^2}\right)\left(1 - \frac{2Gm}{rc^2}\right)^{-1}$$

In the inner crust ($0.1\text{--}0.5\,\rho_0$), competition between nuclear attraction and Coulomb repulsion produces exotic shapes known as "nuclear pasta":

Gnocchi: Spherical nuclei in a neutron gas ($\sim 0.1\,\rho_0$)
Spaghetti: Rod-like nuclear structures ($\sim 0.2\,\rho_0$)
Lasagna: Flat slab-like sheets ($\sim 0.3\,\rho_0$)
Anti-spaghetti: Cylindrical voids in nuclear matter ($\sim 0.4\,\rho_0$)
Anti-gnocchi (Swiss cheese): Spherical voids ($\sim 0.5\,\rho_0$)

Nuclear pasta may constitute up to ~50% of the crust mass and significantly affects neutrino transport, thermal conductivity, and possibly gravitational wave emission from rotating neutron stars.

Nuclear Cosmochronometry

Long-lived radioactive isotopes serve as cosmic clocks to date astrophysical events. The principle uses the decay equation:

$$N(t) = N_0\,e^{-\lambda t} \qquad \Rightarrow \qquad t = \frac{1}{\lambda}\ln\frac{N_0}{N(t)}$$

Key chronometer pairs include:

- $^{232}$Th/$^{238}$U ratio ($T_{1/2} = 14.05/4.47$ Gyr): Ages of metal-poor stars
- $^{235}$U/$^{238}$U ratio ($T_{1/2} = 0.704/4.47$ Gyr): Age of the Solar System
- $^{187}$Re/$^{187}$Os ($T_{1/2} = 41.6$ Gyr): Age of the Galaxy

From the observed $^{235}$U/$^{238}$U ratio of 0.00725 today and the estimated production ratio of ~1.4 from r-process models, the age of the r-process elements (and approximately the Galaxy) is estimated at $\sim 12\text{--}14$ Gyr, consistent with the CMB age of 13.8 Gyr.

Supernova Nucleosynthesis

Core-collapse supernovae ($M \gtrsim 8\,M_\odot$) are responsible for producing many intermediate-mass and heavy elements. When the iron core exceeds the Chandrasekhar mass:

$$M_{\text{Ch}} \approx 1.44\left(\frac{Y_e}{0.5}\right)^2 M_\odot$$

the core collapses in ~100 ms. The gravitational binding energy released is enormous:

$$E_{\text{grav}} \approx \frac{3GM_{\text{NS}}^2}{5R_{\text{NS}}} \approx 3 \times 10^{53}\;\text{erg} \approx 3 \times 10^{46}\;\text{J}$$

About 99% of this energy is carried away by neutrinos. The shock wave ejects the outer layers, producing elements through explosive nucleosynthesis (burning at$T > 5 \times 10^9$ K). Key products include $^{56}$Ni (which decays to$^{56}$Fe and powers the light curve), $^{44}$Ti, and various alpha-elements (O, Ne, Mg, Si, S, Ar, Ca). The neutrino-driven wind from the proto-neutron star may also contribute to the r-process or produce lighter heavy elements (the "weak r-process" or $\nu p$-process).

Python Simulation: Astrophysical Nuclear Processes

Gamow peak visualization for the p+p reaction, pp chain vs CNO cycle energy production rates, solar abundance pattern, and s-process/r-process nucleosynthesis paths.

Nuclear Astrophysics Visualizations

Python

Gamow peak, burning rates, abundances, and nucleosynthesis paths

script.py162 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ──────────────────────────────────────────────
# Nuclear Astrophysics: Gamow Peak &
# Nucleosynthesis Processes
# ──────────────────────────────────────────────

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.patch.set_facecolor('#0a0a1a')

for ax_row in axes:
    for ax in ax_row:
        ax.set_facecolor('#0a0a1a')
        for spine in ax.spines.values():
            spine.set_edgecolor('#334155')

# ─── Panel 1: Gamow Peak ───
# For p+p reaction at T = 1.5e7 K (solar core)
kB = 8.617e-5  # MeV/K
T = 1.5e7  # K
kT = kB * T  # MeV

E = np.linspace(0.001, 0.1, 1000)  # MeV

# Maxwell-Boltzmann distribution
MB = np.exp(-E / kT) / kT

# Gamow factor (Coulomb barrier penetration)
# E_G = (pi * alpha * Z1 * Z2)^2 * 2 * mu * c^2
# For p+p: Z1=Z2=1, mu = m_p/2
m_p = 938.272  # MeV/c^2
mu = m_p / 2
alpha = 1.0/137.036
E_G = (np.pi * alpha * 1 * 1)**2 * 2 * mu  # MeV
# E_G ~ 0.494 MeV for p+p

tunneling = np.exp(-np.sqrt(E_G / E))

# Gamow peak = product
gamow_peak = MB * tunneling
gamow_peak /= gamow_peak.max()
MB_norm = MB / MB.max()
tunneling_norm = tunneling / tunneling[len(tunneling)//4]

# Gamow peak energy
E0 = (0.5 * E_G * kT**2)**(1.0/3.0)

axes[0,0].plot(E*1000, MB_norm, color='#3b82f6', linewidth=2, label='Maxwell-Boltzmann', linestyle='--')
axes[0,0].plot(E*1000, tunneling_norm, color='#22c55e', linewidth=2, label='Tunneling prob.', linestyle='--')
axes[0,0].fill_between(E*1000, gamow_peak, alpha=0.3, color='#ef4444')
axes[0,0].plot(E*1000, gamow_peak, color='#ef4444', linewidth=2.5, label='Gamow peak')
axes[0,0].axvline(E0*1000, color='#f59e0b', linewidth=1, linestyle=':', label=f'E0 = {E0*1000:.1f} keV')

axes[0,0].set_xlabel('Energy [keV]', color='#94a3b8', fontsize=11)
axes[0,0].set_ylabel('Relative probability', color='#94a3b8', fontsize=11)
axes[0,0].set_title('Gamow Peak (p+p, T=15 MK)', color='white', fontsize=13)
axes[0,0].legend(fontsize=8, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0,0].tick_params(colors='#94a3b8')
axes[0,0].set_xlim(0, 50)
axes[0,0].set_ylim(0, 1.2)

# ─── Panel 2: pp chain energy production ───
T_arr = np.linspace(5, 40, 200)  # Million K
# pp chain rate ~ T^4 (approximate)
rate_pp = (T_arr/15)**4
# CNO rate ~ T^16 (approximate)
rate_cno = 0.1 * (T_arr/15)**16
rate_cno = np.clip(rate_cno, 0, 1e6)

axes[0,1].semilogy(T_arr, rate_pp, color='#ef4444', linewidth=2.5, label='pp chain')
axes[0,1].semilogy(T_arr, rate_cno, color='#3b82f6', linewidth=2.5, label='CNO cycle')
axes[0,1].axvline(15, color='#f59e0b', linewidth=1, linestyle=':', alpha=0.5)
axes[0,1].text(16, 1e-3, 'Solar core', color='#f59e0b', fontsize=9)

# Crossover temperature
idx_cross = np.argmin(np.abs(rate_pp - rate_cno))
if idx_cross > 0 and idx_cross < len(T_arr)-1:
    axes[0,1].axvline(T_arr[idx_cross], color='#22c55e', linewidth=1, linestyle='--', alpha=0.5)
    axes[0,1].text(T_arr[idx_cross]+0.5, 1e-5, f'Crossover\n~{T_arr[idx_cross]:.0f} MK',
                   color='#22c55e', fontsize=8)

axes[0,1].set_xlabel('Temperature [Million K]', color='#94a3b8', fontsize=11)
axes[0,1].set_ylabel('Relative energy production rate', color='#94a3b8', fontsize=11)
axes[0,1].set_title('pp Chain vs CNO Cycle', color='white', fontsize=13)
axes[0,1].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0,1].tick_params(colors='#94a3b8')
axes[0,1].set_ylim(1e-6, 1e6)

# ─── Panel 3: Solar abundance pattern ───
# Simplified abundance data (log scale, Si=6)
elements = ['H', 'He', 'C', 'N', 'O', 'Ne', 'Mg', 'Si', 'S', 'Ca', 'Fe', 'Ni']
Z_vals = [1, 2, 6, 7, 8, 10, 12, 14, 16, 20, 26, 28]
abundances = [12.0, 10.93, 8.43, 7.83, 8.69, 7.93, 7.60, 7.51, 7.12, 6.34, 7.50, 6.22]

axes[1,0].bar(range(len(elements)), abundances, color=['#ef4444' if z <= 2 else
              '#f59e0b' if z <= 8 else '#22c55e' if z <= 14 else '#3b82f6'
              for z in Z_vals], alpha=0.8)
axes[1,0].set_xticks(range(len(elements)))
axes[1,0].set_xticklabels(elements, color='#94a3b8', fontsize=9)
axes[1,0].set_ylabel('log(abundance) [H=12]', color='#94a3b8', fontsize=11)
axes[1,0].set_title('Solar System Abundances', color='white', fontsize=13)
axes[1,0].tick_params(colors='#94a3b8')

# Annotate Fe peak
axes[1,0].annotate('Fe peak', xy=(10, 7.5), xytext=(8, 9),
                   arrowprops=dict(arrowstyle='->', color='#f59e0b'),
                   color='#f59e0b', fontsize=9)

# ─── Panel 4: s-process and r-process paths (simplified) ───
# Show neutron number vs proton number
Z_stable = np.array([26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38])
N_stable = np.array([30, 32, 30, 34, 34, 38, 41, 42, 45, 44, 48, 48, 50])

# s-process path (near stability)
Z_s = np.array([26, 27, 28, 28, 29, 30, 30, 31, 32, 32, 33, 34, 34, 35, 36, 36, 37, 38])
N_s = np.array([30, 32, 30, 32, 34, 34, 36, 38, 38, 41, 42, 44, 45, 44, 46, 48, 48, 50])

# r-process path (neutron rich)
Z_r = np.array([26, 27, 28, 28, 29, 30, 30, 31, 32, 32, 33, 34, 34, 35, 36, 36, 37, 38])
N_r = N_s + np.array([4, 5, 6, 7, 6, 7, 8, 7, 8, 9, 8, 8, 9, 9, 9, 10, 10, 10])

axes[1,1].scatter(N_stable, Z_stable, color='#94a3b8', s=60, marker='s',
                  label='Stable nuclei', zorder=3)
axes[1,1].plot(N_s, Z_s, color='#f59e0b', linewidth=2, marker='o', markersize=4,
               label='s-process path', zorder=2)
axes[1,1].plot(N_r, Z_r, color='#3b82f6', linewidth=2, marker='^', markersize=4,
               label='r-process path', zorder=2)

axes[1,1].set_xlabel('Neutron number N', color='#94a3b8', fontsize=11)
axes[1,1].set_ylabel('Proton number Z', color='#94a3b8', fontsize=11)
axes[1,1].set_title('s-process and r-process Paths', color='white', fontsize=13)
axes[1,1].legend(fontsize=9, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white',
                 loc='upper left')
axes[1,1].tick_params(colors='#94a3b8')

plt.tight_layout()
plt.savefig('nuclear_astrophysics.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())

print("=== Nuclear Astrophysics Summary ===")
print()
print(f"Solar core temperature: T = {T/1e6:.0f} MK")
print(f"kT = {kT*1000:.2f} keV")
print(f"Gamow energy (p+p): E_G = {E_G*1000:.1f} keV")
print(f"Gamow peak energy: E_0 = {E0*1000:.1f} keV")
Delta_E = 4 * np.sqrt(E0 * kT / 3)
print(f"Gamow peak width: Delta = {Delta_E*1000:.1f} keV")
print()
print("Hydrogen burning energy release:")
print("  pp chain: 4p -> He-4 + 2e+ + 2nu_e + 26.73 MeV")
print("  CNO cycle: Same net reaction, C,N,O as catalysts")
print(f"  pp dominates below ~{T_arr[idx_cross]:.0f} MK, CNO above")
print()
print("Key burning stages:")
print("  H burning:  T ~ 15 MK,    tau ~ 10 Gyr (main sequence)")
print("  He burning: T ~ 100 MK,   tau ~ 1 Gyr (red giant)")
print("  C burning:  T ~ 600 MK,   tau ~ 1000 yr")
print("  Si burning: T ~ 3000 MK,  tau ~ days")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation: Gamow Peak Calculator

Computes Gamow peak parameters for various thermonuclear reactions, tabulates stellar burning stages, and compares s-process and r-process characteristics.

Nuclear Astrophysics Calculator

Fortran

Gamow energies, burning stages, and nucleosynthesis process comparison

program.f90127 lines

program nuclear_astrophysics
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: pi = 3.14159265358979323846_dp
  real(dp), parameter :: kB = 8.617d-5  ! MeV/K
  real(dp), parameter :: alpha_fs = 1.0_dp/137.036_dp
  real(dp), parameter :: m_p = 938.272_dp  ! MeV/c^2
  real(dp), parameter :: m_n = 939.565_dp  ! MeV/c^2
  real(dp), parameter :: hbar_c = 197.327_dp  ! MeV fm

real(dp) :: T, kT, E_G, E0, Delta, mu_red
  real(dp) :: Z1, Z2, A1, A2, tau
  real(dp) :: rate_factor
  integer :: i

write(*,'(A)') '=================================================='
  write(*,'(A)') ' Nuclear Astrophysics Calculations'
  write(*,'(A)') '=================================================='
  write(*,'(A)') ''

! ─── Gamow Peak Calculations ───
  write(*,'(A)') '--- Gamow Peak Parameters ---'
  write(*,'(A)') ''

T = 1.5d7  ! K (solar core)
  kT = kB * T

write(*,'(A,ES10.2,A)') '  Solar core T = ', T, ' K'
  write(*,'(A,F8.4,A)') '  kT = ', kT*1000.0_dp, ' keV'
  write(*,'(A)') ''

write(*,'(A20,A10,A10,A12,A12,A12)') 'Reaction', 'Z1*Z2', &
       'mu(u)', 'E_G(keV)', 'E_0(keV)', 'Delta(keV)'
  write(*,'(A)') '──────────────────────────────────────────────────' // &
       '─────────────────────────'

! p + p
  Z1 = 1; Z2 = 1; A1 = 1; A2 = 1
  mu_red = A1*A2/(A1+A2) * 931.494_dp  ! MeV/c^2
  E_G = (pi * alpha_fs * Z1 * Z2)**2 * 2.0_dp * mu_red
  E0 = (0.5_dp * E_G * kT**2)**(1.0_dp/3.0_dp)
  Delta = 4.0_dp * sqrt(E0 * kT / 3.0_dp)
  write(*,'(A20,F10.0,F10.3,F12.1,F12.2,F12.2)') 'p + p', Z1*Z2, &
       A1*A2/(A1+A2), E_G*1000, E0*1000, Delta*1000

! p + C-12 (CNO)
  Z1 = 1; Z2 = 6; A1 = 1; A2 = 12
  mu_red = A1*A2/(A1+A2) * 931.494_dp
  E_G = (pi * alpha_fs * Z1 * Z2)**2 * 2.0_dp * mu_red
  E0 = (0.5_dp * E_G * kT**2)**(1.0_dp/3.0_dp)
  Delta = 4.0_dp * sqrt(E0 * kT / 3.0_dp)
  write(*,'(A20,F10.0,F10.3,F12.1,F12.2,F12.2)') 'p + C-12', Z1*Z2, &
       A1*A2/(A1+A2), E_G*1000, E0*1000, Delta*1000

! He-4 + He-4 (triple alpha, first step)
  Z1 = 2; Z2 = 2; A1 = 4; A2 = 4
  mu_red = A1*A2/(A1+A2) * 931.494_dp
  E_G = (pi * alpha_fs * Z1 * Z2)**2 * 2.0_dp * mu_red
  T = 1.0d8  ! He burning temperature
  kT = kB * T
  E0 = (0.5_dp * E_G * kT**2)**(1.0_dp/3.0_dp)
  Delta = 4.0_dp * sqrt(E0 * kT / 3.0_dp)
  write(*,'(A20,F10.0,F10.3,F12.1,F12.2,F12.2)') 'He4+He4 (T=100MK)', Z1*Z2, &
       A1*A2/(A1+A2), E_G*1000, E0*1000, Delta*1000

! C-12 + C-12 (C burning)
  Z1 = 6; Z2 = 6; A1 = 12; A2 = 12
  mu_red = A1*A2/(A1+A2) * 931.494_dp
  E_G = (pi * alpha_fs * Z1 * Z2)**2 * 2.0_dp * mu_red
  T = 6.0d8
  kT = kB * T
  E0 = (0.5_dp * E_G * kT**2)**(1.0_dp/3.0_dp)
  Delta = 4.0_dp * sqrt(E0 * kT / 3.0_dp)
  write(*,'(A20,F10.0,F10.3,F12.1,F12.2,F12.2)') 'C12+C12 (T=600MK)', Z1*Z2, &
       A1*A2/(A1+A2), E_G*1000, E0*1000, Delta*1000

write(*,'(A)') ''

! ─── Stellar Burning Stages ───
  write(*,'(A)') '--- Stellar Burning Stages (25 M_sun star) ---'
  write(*,'(A)') ''
  write(*,'(A15,A12,A15,A15,A20)') 'Stage', 'T(MK)', 'Duration', &
       'Q(MeV/nuc)', 'Products'
  write(*,'(A)') '──────────────────────────────────────────────────' // &
       '─────────────────────────'
  write(*,'(A15,A12,A15,F15.2,A20)') 'H burning', '15', '10 Myr', &
       6.68_dp, 'He-4'
  write(*,'(A15,A12,A15,F15.2,A20)') 'He burning', '100', '1 Myr', &
       1.60_dp, 'C-12, O-16'
  write(*,'(A15,A12,A15,F15.2,A20)') 'C burning', '600', '1000 yr', &
       0.54_dp, 'Ne, Na, Mg'
  write(*,'(A15,A12,A15,F15.2,A20)') 'Ne burning', '1200', '0.7 yr', &
       0.34_dp, 'O, Mg'
  write(*,'(A15,A12,A15,F15.2,A20)') 'O burning', '1500', '2 yr', &
       0.32_dp, 'Si, S'
  write(*,'(A15,A12,A15,F15.2,A20)') 'Si burning', '3000', 'days', &
       0.18_dp, 'Fe, Ni'
  write(*,'(A)') ''

! ─── Neutron Capture Processes ───
  write(*,'(A)') '--- s-process vs r-process ---'
  write(*,'(A)') ''
  write(*,'(A20,A25,A25)') '', 's-process', 'r-process'
  write(*,'(A)') '──────────────────────────────────────────────────' // &
       '─────────────────────'
  write(*,'(A20,A25,A25)') 'Site', 'AGB stars', 'NS mergers/SN'
  write(*,'(A20,A25,A25)') 'Neutron density', '~10^7 cm^-3', '~10^24 cm^-3'
  write(*,'(A20,A25,A25)') 'Timescale', '~years', '~seconds'
  write(*,'(A20,A25,A25)') 'Path', 'Along stability', 'Far from stability'
  write(*,'(A20,A25,A25)') 'End point', 'Bi-209', 'U,Th,Cf'
  write(*,'(A20,A25,A25)') 'Key isotopes', 'Ba,Sr,Pb', 'Eu,Pt,Au,U'
  write(*,'(A)') ''

! ─── Nuclear Pasta Phases ───
  write(*,'(A)') '--- Neutron Star Crust: Nuclear Pasta ---'
  write(*,'(A)') ''
  write(*,'(A20,A15,A25)') 'Phase', 'Density', 'Structure'
  write(*,'(A)') '────────────────────────────────────────────────────────'
  write(*,'(A20,A15,A25)') 'Gnocchi', '~0.1 rho_0', 'Spherical nuclei'
  write(*,'(A20,A15,A25)') 'Spaghetti', '~0.2 rho_0', 'Rod-like nuclei'
  write(*,'(A20,A15,A25)') 'Lasagna', '~0.3 rho_0', 'Slab-like nuclei'
  write(*,'(A20,A15,A25)') 'Anti-spaghetti', '~0.4 rho_0', 'Cylindrical holes'
  write(*,'(A20,A15,A25)') 'Anti-gnocchi', '~0.5 rho_0', 'Spherical holes'
  write(*,'(A20,A15,A25)') 'Uniform', '> 0.5 rho_0', 'Uniform nuclear matter'
  write(*,'(A,A)') '  (rho_0 = 2.8e14 g/cm^3', ' = nuclear saturation density)'

end program nuclear_astrophysics

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Radiation Detection Course Overview →