notes-lecture-nuclear-burning.org

400A - Nuclear burning

Materials: Onno Pol’s lecture notes Chapter 6.1 and 6.2, Kippenhahn book Chapter 9 and 18, Cox & Giuli vol. I, Chapter 17.7, Clayton Chapter 4, Gamow 1928, Arnould & Takahashi 1999, Jose & Iliadis 2011

Microphysics of nuclear burning

Summary of where we are

We have four equations of the stellar structure assuming spherical symmetry, LTE, and hydrostatic (∂ t ≡ 0) equilibrium:

Mass conservation

Hydrostatic equilibrium

which follows from the momentum conservation equation.

Equation of state

Energy transport

where ∇ =∂ log(T)/∂ log(P) is the local temperature gradient, equal to the radiative gradient in stably stratified regions:

with κ = (1/κ_rad + 1/κ_cond)^-1 the combination “in parallel” of the radiative and conductive opacity (assumed to be known from atomic physics), and ∇ ≡ ∇_ad the adiabatic gradient (within ∼10^-7-8 precision) for convective regions. We also have a criterion (Schwarzschild or Ledoux) to determine which region is which.

Energy conservation

with ε_grav = T∂ s/∂ t the change in internal energy (typically ≅ 0 in gravothermal equilibrium, that is for phases when the evolution is much slower than τ_KH), and ε_ν>0.

The next step

With all the equations above we almost have a complete set of solvable differential equations. What is left to discuss is the specific energy generation per unit time ε_nuc due to nuclear burning (which will occupy us in this and the next lecture) and the specific energy loss ε_ν due to neutrino losses (which we will treat later).

In doing so, we will introduce another time-dependent equation (besides the implicit dependence in ε_grav): as nuclear burning proceeds, the chemical composition changes, on a timescale that we will see is extremely long compared to the thermal and dynamical timescales.

Energy generation as a self-regulating process

We have seen in the virial theorem lecture that one can derive from first principle that stars have a “negative heat capacity” and obtain a scaling between mass, radius and average temperature:

where C was a constant depending on the details of the mass distribution.

The star has a finite temperature, thus it radiates away thermal energy at the surface in the form of a luminosity L. This means the gravitational potential energy is decreasing (because of the virial theorem!), thus it must contract leading to an increase in T. For the Sun this would go on for ∼ 15 Myrs before it runs out of all of its internal energy, which is way too short a timescale to explain many Earthly observables.

We know this must stop either because the ideal gas EOS does not work (QM effects → degeneracy pressure → white dwarfs), or because some other energy source intervenes, balancing out the energy losses at the surface. Stars don’t shine because they burn, stars burn because they shine.

Initially, people considered chemical reactions and radioactivity as a possible energy sources, but it was shown very quickly that these are insufficient.

• Q: how would you show that these sources are insufficient? (Hint: consider the φ factor for the amount of mass released in the equations below)

It took major advances in quantum mechanics (QM) and their application to atomic and nuclear theory to work out in the late 1930s that the energy source is nuclear burning, a result obtained by Hans Bethe and collaborators.

In a sense, you can think of a star as an inertial confinement nuclear reactor where the confinement is provided by the self-gravity. This also implies that nuclear burning in a star is a self-regulating process: because nuclear reactions are a consequence of the energy losses, during long-lived equilibrium phases of the stellar life, they produce just enough energy to balance the losses! In other words, L_nuc = L ( + L_ν).

If a star were to not produce enough energy to verify that (i.e., L>L_nuc), then it will lose too much energy, meaning it will contract, and increase it’s temperature because of the virial theorem, which in turn regulates the energy generation by nuclear reactions (as we will see in a moment), until L_nuc reaches L. Viceversa, if L_nuc > L for some reason (e.g., there is a thermonuclear explosion in the star), than the extra energy release will cause an expansion of the star and by the virial theorem lower the mean temperature until L_nuc = L.

The nuclear timescale

Naturally, nuclear energy generation consumes nuclear fuel: to provide L_nuc the composition of the star slightly changes in time. This drives the evolution of star.

We can estimate the timescale for this assuming the star is in gravothermal equilibrium, so L=L_nuc. The nuclear timescale is the time it takes to lose the energy generated by nuclear reactions:

where φ is the fraction of rest mass of nuclei converted in energy by nuclear burning, f_burn is the fraction of the stellar mass M that is affected by burning (we need a stellar model to estimate that). For the Sun, f_burn ≅ 0.1 from detailed models (as you can verify with your MESA-web model!).

φ depends on the nuclear physics details. For hydrogen burning into helium (the two lightest elements), the proton mass is m_p = 1.0081 m_u (where the atomic mass unit m_u is defined in such a way that the mass of $¹²\mathrm{C}$ is exactly 12m_u - this is more convenient to measure experimentally to make a standard), and the mass of helium 4 is m(⁴He) = 4.0039m_u, so the fraction of rest mass of 4 protons turning into a helium nucleus is:

Note that to turn 4 protons into a helium, because of charge and leptonic number conservation, there needs to be 2 positrons and 2 neutrinos produced! Plugging in φ and M_o and L_o we get for the nuclear timescale:

So we now have estimates that allow us to see the complete ordering for the main global timescales for stellar evolution:

Moreover, Eq. \ref{eq:tau_nuc_scaling} shows that when accounting for nuclear energy generation (which we have implicitly assumed here to give a value of φ), the nuclear timescale matches the timescale found from geological evidence on Earth, and the ordering \ref{eq:timescale_ordering} tells us that the evolution of the Sun is on a very slow timescale compared to thermal and dynamical timescale, validating the assumptions we made so far to derive the equations for stellar structure and allowing for a quasi-static approximation to deal with the stellar evolution.

Let’s now dive into the details of the nuclear physics that allow for nuclear burning.

The energy reservoir: binding energy per nucleon

The energy reservoir that stars tap into is the nuclear binding energy: nuclear reactions rearrange nucleons in nuclei to create more bound configuration and extract binding energy as heat source. Therefore, it is useful to consider the nuclear binding energy per nucleon of all nuclei in the periodic table:

This plot shows empirical data, and there are several notable features coming from the internal structure of the nuclei that any model of nuclear (strong) interactions needs to explain:

• the binding energy per nucleon raises sharply with A among the light nuclei;
• the helium nucleus (α particle) has an exceptionally high binding energy per nucleon of ∼ 7.5MeV/nucleon;
• there are peaks for nuclei with A divisible by 4 and N=Z, that is nuclei that can be approximately thought as bound collections of α particles (for example $¹²\mathr{C}$ ∼ 3 α particles bound together, $¹⁶\mathrm{O}$ ∼4 α particles, etc..);
• the nuclear binding energy is roughly constant at about ∼ 8MeV/nucleon for most heavy nuclei;
• the slight drop in B/A after the maximum is caused by the Coulomb repulsion between the protons in heavy nuclei, and this is why those nuclei require a higher number of neutrons than protons to hold together, N>Z.
• there is a local maximum of the binding energy per nucleon at the iron group, around $⁵⁶\mathr{Fe}$/$⁶²\mathrm{Ni}$.

The presence of a maximum implies that there cannot be energy release from the fusion of nuclei heavier than iron: indeed those heavy nuclei are typically fuel for nuclear fission reactors, where energy is extracted breaking them apart. Stars, which do nuclear burning to release energy and balance the losses at the surface, have no interest in producing element heavier than iron during their stable/hydrostatic lifetime. The question of the formation of elements heavier than iron (which clearly exist!) is something that is actively worked on and requires out-of-equilibrium processes that can only occur in the presence of a neutron rich environment (e.g., AGB stellar winds, neutron stars interacting with something else, etc.). We may have a guest lecture on this later on, and there are projects related to these.

Conversely, moving leftward on this chart, by fusing together light nuclei, stars can release nuclear binding energy and sustain themselves. This is also convenient since stars are mostly made of H and He, so they have a lot of light elements available as energy sources.

The change in binding energy per nucleon Δ(B/A) is what powers stars, and we can write the binding energy as the difference between the sum of the masses of the Z protons (m_p) and N=A-Z neutrons (m_n) minus the mass of the nucleus (measured experimentally):

which is a positive quantity precisely because the strong interaction between protons and neutrons keeps them bound. In the second step, we approximate the proton mass with the hydrogen atom mass, making a mistake of the order of ∼ 10 eV by neglecting the binding energy of the electron (which is acceptable since we are dealing with MeV = 10⁶eV as the relevant nuclear scale).

Because of the apparent peaks in the B/A vs. A curve, stars will tend to produce more of the elements with the most tight nuclei, whose production releases more energy! The nuclear structure determines the most abundant chemicals in the Universe.

Moreover, from the fact that there is a sharp jump form H to He, and then the curve rapidly flattens, we know that the rate at which nuclei are converted to achieve the equilibrium condition L_nuc = L will be slower when H is converted into He, and then it will need to rapidly speed up when He is converted in carbon and onwards (because the energy release per nucleon significantly drops lowering the factor φ in τ_nuc).

In the very late evolutionary phases (post helium burning), the required burning rate may become so high that τ_nuc becomes shorter than the global thermal timescale (but luckily not of the local thermal timescale, so all the equations we have derived so far still hold, since they are differential equations that describe local quantities). We will return on this later.

Q-value

For any nuclear reaction that may happen a+X → b+Y with a,b,X,Y generic particles, it is useful to define the so called Q-value as the mass difference between the reagents and the products:

which, if Q>0 is the energy released by the reaction which is thus exoenergetic - of the kind that stars need to generate energy and sustain themselves against gravity, or if Q<0 is the energy input needed to get the endoenergetic reaction going.

N.B.: In exoenergetic nuclear reactions, the energy release comes from the mass defect, caused by the increase in binding energy per nucleon: the total mass of the outgoing particles is lower than the ones incoming because their binding energy is higher or in other words their total energy is more negative.

The energy release by nuclear reaction per unit time and unit mass is just:

where the sum runs over all the possible reactions, r_i is the rate per unit time and volume of the reactions, and the division by the mass density ρ gives the right dimensions [ε_nuc] = [E]/([t][M]).

So, what is left to do is calculate the volumetric reaction rate r_i that can occur in a star.

Variety of possible nuclear reactions

A generic nuclear reaction X+a → Y+b is often written as X(a,b)Y to make it easy to express chains of reactions, e.g., X(a,b)Y(c,d)Z(e,f)A…

Depending on the nature of the incoming particle (X and a in our generic reaction), or in other words on the microphysics that determines the interaction, there can be of various kinds of reactions.

Charged-particles reactions

When X and a are charged nuclei, then the reaction can only occur if something allows them to overcome the Coulomb repulsion. These can be resonant or not (the distinction will come back later).

Example:

Reactions involving neutrons

In this case the force involved is the strong force, and there is no Coulomb repulsion to overcome. However, these require an environment that is neutron rich, which is astrophysically a rare occurrence, since the neutron is an unstable particle that decays in ∼15min to a proton n→ p+e^-+$\bar{ν_e}$. This half-life however can significantly change for neutrons bound in nuclei as opposed to free neutrons, that is the β^- decay time of a neutron rich nucleus can be much longer than the half-life of a free neutron.

Depending on the available flux of neutrons in the environment, we distinguish:

• r-process for rapid neutron captures (i.e., the rate of neutron captures is high w.r.t. the rate of neutron decays)
• s-process for slow neutron captures (i.e., each nucleus captures at best one neutron before decaying).

These processes are involved in the formation of elements heavier than iron, but they require particular astrophysical environment (e.g., the merger of two neutron stars or a neutron star with the core of another star, or the envelope of an AGB star).

Weak reactions

These can typically be spotted by the presence of a neutrino and/or the conversion of a nucleon from one eigenstate of isospin to another (in simpler words, the conversion of a proton into a neutron or viceversa).

Example:

Photodisintegrations

When one of the particles is a photon and the outgoing particles can be seen as “fragments” of the ingoing nucleus. These can occur when very energetic γ ray photons, because their energy needs to be comparable to the binding energy of nuclei, of the order of ∼8MeV × A, can encounter particles. This can occur for example at the very late moments of massive star evolution.

Example:

Nuclear reactions in stars

All of the types of reactions listed above (and more) can occur at some point in the evolution (and explosion!) of stars. For example, during hydrogen core burning (which we have used to estimate φ and thus τ_nuc) the star burns 4 protons into α particles:

where the positrons need to be there to conserve electric charge throughout the process, and the neutrinos need to be there for conserving the leptonic number (+1 for the leptons electron e^-, muon μ^-, tau τ^- and the corresponding neutrinos ν_e, ν_μ, ν_τ and -1 for their antiparticles positron e, positive muon μ, and positive τ⁺ and the corresponding antineutrinos).

N.B.: because of the phenomenon of neutrino oscillations (i.e., the conversion of ν_e ↔ ν_μ or ν_τ) the leptonic number conservation is not exact in nature, or in other words, the symmetry associated to this conservation law is not exact. While this was discovered through the ”missing solar neutrinos problem”, and it is thus related to stellar/solar physics, it requires the propagation of neutrinos over distances much longer than the size of nuclei, therefore, for the purpose of treating nuclear reactions we can assume conservation of the leptonic number.

From Eq. \ref{eq:overall_h_burn} we can see that:

• protons need to encounter each other. Statistically, 4 protons are unlikely to meet each other at a point in space for reaction Eq. \ref{eq:overall_h_burn} to occur. Eq. \ref{eq:overall_h_burn} is a compound reaction that “summarizes” the more complex burning process of H into He that we will detail later. Nevertheless, the process will necessarily involve charged particle reactions.
• ν_e appear ⇒ there will be weak reactions involved
• as already seen earlier, we can calculate the φ factor (cf. Eq. \ref{eq:mass_defect}), and thus the Q-value for the overall reaction Q_{H burn} ≅ 26.5 MeV. Note that the Q-value does not really depend on the details of the burning process.

N.B.: we discuss Eq. \ref{eq:overall_h_burn} because H is the most abundant element in the Universe, that most stars are made of, and thus this is (typically) the first process occurring in stars. It is also the one with the highest Q-value (cf. B/A vs. A plot!), thus the one that for a given luminosity L of the star can proceed at the slowest rate and last the longest.

Charged particle reactions

For the rest of this lecture we will focus mostly on charged particle reactions: as we just saw these are involved since H core burning, and they are the main reactions during the hydrostatic lifetime of stars. Moreover, breakthroughs in QM (by George Gamow) applied to the interaction of charged particles are what opened the way to the discovery that the energy source in stars are nuclear reactions.

Bohr’s approximation

To discuss them, we will implicitly use Bohr’s approximation, which is not a completely accurate physical picture, but it is intuitive and allows to describe the main processes occurring in nuclei. In this approximation, we treat the generic reaction between charged particles X(a,b)Y as if it went through an intermediate step of forming a compound nucleus C:

where the nucleus C is formed in an excited state C^* since in the reference frame of X the particle a arrives with its own kinetic energy and internal binding energy that (together with X’s internal energy) are generally not exactly the total energy for the compound nucleus C.

The second step is the decay of this fictional compound excited particle C^* in the products Y and b.

In this approximation, the second step of the decay does not depend on the first step (in other words, C^* loses the memory of how it formed) as long as the half-life of C^* is long compared to the light-crossing time of C^* itself. This is because we consider charged particles by hypothesis, so electromagnetic forces mediated by photons are what determines the interactions of the building blocks of C^*, and on a timescale long compared to the light-crossing time they will equilibrate and lose memory of how they came to be.

The light crossing time of a nucleus can be estimated starting from the experimental result on the size of nuclei (something that also needs to be explained by models of the strong force):

Any compound nucleus C^* with lifetime longer than this allows us to use this two step approximation to treat the problem, where the two incoming particles X and a come “into contact” (i.e., within their r_n), form an excited compound nucleus C^*, which then de-excites in the final products Y and b independently of how it formed.

N.B.: the nuclear radius dependence on A can be flipped around to infer that the average density of nuclei is constant as A increases: ρ_n ≅ Am_u / (4π/3 × r_n³) ≅ 2×10¹⁴ A g cm^-3.

N.B.: this is necessarily an oversimplified picture, since nuclei are described by QM and don’t have a “hard” boundary at r_n, but rather their constituent nucleons have decaying wave-functions that solve the Schrodinger equation with the nuclear potential for their interactions. In reality, each particle can be described as a wave with De Broglie wavelength λ = h/p, and a physically more accurate picture should treat all the particles involved accounting for their wave nature.

To understand how two charged nuclei, both with positive charge, can “come into contact” within r_n, we need to consider the potential governing their interaction:

where the last term is the centrifugal potential in the rest-frame of the target nucleus X which depends on the reduced mass m_aX = m_am_X/(m_a+m_X) and quantum number \ell which determines the order of the wave-function of the system a+X. For simplicity, we can limit ourselves to consider \ell=0: we already have a repulsive Coulomb term to win over, and any extra repulsive term such a centrifugal barrier is only going to lower the reaction rate. The most important reactions are going to have \ell = 0, that is head-on collisions between a and X.

Electromagnetic potential term

For the electromagnetic term we can write:

where the first part is >0 and describes the Coulomb repulsion between the two nuclei of chage Z_ae and Z_Xe (both positive), and the electron screening term reduces the repulsion of the nuclei: in the stellar plasma we expect each nucleus to be statistically surrounded by a “cloud” of electrons of radius of the order of the Debye length of the plasma:

where N is the total number of particles in absence of screening (nuclei/ions+electrons), and χ = ∑_i Z_i² (N_i/N) + N_e/N with N_i and N_e number of ions and electrons in absence of screening.

For distances between a and X larger than r_Debye the electron screening reduces the Coulomb repulsion between the nuclei.

Nuclear potential term

Finally, for the nuclear potential, there isn’t a well known functional form from first principles, and it is typically derived experimentally. This is because the interactions between nucleons (=protons and neutrons) cannot be treated in a perturbative theory of the strong force. For more details than necessary to understand thermonuclear reactions, see also this appendix and nuclear physics textbooks such as ”Introductory nuclear physics” by K. Krane.

Besides the nuclear physics details here, the important point is that the nuclear potential is going to be attractive at short range (V_nuc(r≤ r_n)<0), but it has a repulsive core (that is there is a certain r_nuc,core below which V_nuc becomes very large and positive), otherwise the nuclei would not have a finite approximately constant density, and goes to zero at large distances (the strong force has a short range).

Combining electromagnetic and nuclear potential

Putting things together we can sketch the following graph for the potential felt by particle a and generated by the strong and electromagnetic force by particle X:

• at distances r \gg r_Debye electron screening nullifies the repulsive Coulomb potential
• for a relative energy at infinity of E, there is a distance of classical minimum approach r_c
• just outside r_n≡ r_n there is the maximum height of the Coulomb barrier E_C = Z_a Z_X e²/r_n.
• inside r_n (the nuclear radius of particle X), the potential is attractive, and allows for bound states with quantized energy levels.
• if the two nuclei get too close to each other, there is a repulsive core of the nuclear force that dominates over any electromagnetic effect.
• For 0<E<E_c, there are metastable energy level possible (represented in the figure by the gray bands). What makes them metastable is also what allows nuclear burning: quantum tunnelling through the Coulomb barrier.

Impossibility of nuclear reactions without QM

Without QM, for a nuclear reaction to happen (assuming Bohr’s approximation), the two charged particles would need a relative energy at infinity higher than the maximum of the Coulomb barrier, so that r_c \leq r_n. Assuming the energy is just coming from the thermal energy of the gas:

where we assume Z_a = Z_x = A = 1 to minimize the Coulomb barier, so r_n ≡ r₀ = 1.3 × 10^-13 cm. The temperature threshold we have derived is much larger than the mean temperature in the Sun as estimated with the Virial theorem. It is also much larger than the central temperature of the Sun which one could estimate assuming \langle T \rangle - T_surf,o ∼ T_center,o - \langle T \rangle. The Sun is not hot enough to have nuclear energy generation without QM.

N.B.: In a nuclear fusion experiment energies E\geq E_C are reached, however, the beams are not in thermal equilibrium. The key point is laboratory nuclear physics experiment occur at much higher energies than nuclear reactions in stars.

Charged particles nuclear reaction tunneling through the Coulomb barrier

The piece of puzzle that allows for charged particles nuclear reactions in stars is the QM tunnel effect, which was studied in the context of the C^* → Y+b decay in Gamow 1928.

N.B.: the tunnel effect is a purely wave mechanics phenomenon that has to do with constructive interference of waves. The QM element is that the particles involved here are waves!

From QM, we know that in the “classically forbidden region”, where E<V (so the classical kinetic energy term in E=K+V would be K<0), the wave function can still be non-zero. To calculate the wave-function Ψ of particle a in the potential of X, we can make the following ansatz:

• in the classically allowed region r\geq r_c the wave function will have the form of a propagating wave with phase dependent on (E-V(r))^1/2;
• in the classically allowed region with r\leq r_n, we assume the same functional form;
• in the classically forbidden region, we will have a superposition of evanescent waves with exponentially decaying amplitude

N.B.: The math simplifies significantly using the WKB approximation, that is writing Ψ = exp(Φ) and solving for Φ instead of Ψ.

By imposing that the wave function Ψ and its derivative ∂_rΨ are continuous at the boundaries r_c and r_n one can find solutions that have a non-zero amplitude inside r_n, that is tunneling solutions! The continuity of Ψ depends on the phase at the boundaries, for specific values of E it is possible to obtain solutions, these typically correspond to metastable energy levels of the potential, with energy E_metastable ± Δ E and Δ E ∼ h/τ the “width” in energy depending on the τ the half-life of the metastable state. N.B.: the width Δ E of the levels is set by the uncertainty principle.

These states are metastable because a bound system between a and X (i.e., C^*) can result in the tunneling of a out of the potential of X (the decay C^* → a+X).

Therefore, the nuclear reaction rates are going to be extremely sensitive to the relative energy of a and X: if this relative energy E hits a resonance of the compound state C^*, that is if E∼ E_metastable, the wave function of the system has non-zero amplitude inside the nucleus, and since the probability of finding a within r≤ r_n is ∝ ||Ψ||², this means there will be a non-zero probability of forming the compound nucleus.

In proximity of a resonance the cross section takes the typical shape of a Lorentzian with width determined by the energy width of the metastable state C^*{}:

where Γ is the “width” in energy of the metastable level, determined by the uncertainty principle: Γ\tau ∼ h with τ half-life of the metastable state and h Planck’s constant.

Conversely, for non-resonant reactions, since the thermal energies are of order of 100 of eV while the nuclear scale is at ∼10MeV we can neglect the energy dependence of σ_nuc.

Finally, actually carrying out the calculation give a probability of tunneling

with b = 2π Z_a Z_X e² (m_aX/2)^1/2 / \hbar.

N.B.: Nuclear resonances allow stars to work, and for example the existence of a specific resonance in the nucleus of $¹²\mathrm{C}$ (Hoyle state) is what allows helium to burn into carbon and ultimately allows for life (as we will see in more detail in the next lecture)!

However, because nuclear physics experiments operate at so much higher energy than stars, extrapolating to lower relative energies is complicated and dangerous: it is easy to miss an unknown resonance that would change significantly the rate.

Thermonuclear reaction rates

Now that we have discussed how a charged particle nuclear reaction is possible through quantum tunneling through the Coulomb barrier, to complete our quest for ε_nuc we need to obtain the rate per unit volume and time of each possible nuclear reaction.

What brings together the generic nuclei a and X, giving them the energy E such as the probability of tunneling through the Coulomb barrier is not negligible is the thermal motion of the gas they compose: this is why we talk about thermonuclear reactions in a star (and why the energy scale for these reactions is much lower than the energy scale in laboratory experiments).

Once again, to put together an expression for the nuclear reaction rate we can start from dimensional analysis trying to combine the available pieces:

• number density of reactants [n_a] = [n_X] = [L]^-3
• their relative velocity [v] = [L]/[t] (which is related to their relative energy at infinity E which is of the order of the thermal energy, v∼(2k_BT/m_aX)^1/2)
• the cross section for the reaction [σ] = [L²] (whose calculation will depend on the details of the QM problem outlined above and we know will depend on the relative energy of the particles, and thus ultimately their relative - thermal - velocity at infinity): σ ≡ σ(v)

With these ingredients we can make a rate of the number of reactions X(a,b)Y per unit time and volume with:

where we implicitly assumes that a ≠ X. However, in a star sometimes there are reactions among identical particles (for example the weak reaction p+p→ D+e⁺ +ν_e). In such cases we should make sure to not double count particle pairs, so we can write more generally:

where δ_aX = 1 ⇔ a=X. This would be the expression if all particles had the same relative velocity v. In reality, we know the stellar gas is pretty close to LTE and thus the distribution in energy of particles is given by the Maxwell-Boltzmann distribution, and to get the effective rate of reactions we need to integrate over that.

The relative velocity is v=|v_a - v_X| (the velocities of the two species can be different if they have different masses), and substituting to the number density the integral over the velocities of the phase space densities n_i → ∫ dn_i(v) (N.B.: we have already done this many times when deriving the EOS, just in momentum instead of velocity), our expression for the rate becomes:

where dn_i(v_i) are Maxwell-Boltzmann distributions, unless we are considering thermo-nuclear burning in a (partially) degenerate environment.

• Q: can you think of stellar situations where there is burning in a (partially) degenerate environment?

We can explicit the Maxwell-Boltzmann distribution assuming that the nuclei are non-relativistic, since their thermal kinetic energy is of the order of k_bT \ll GeV ∼ m_u/c², and you can analytically verify that the product of two Maxwell-Boltzmann distribution keeps the same functional form by changing variables to express things in the center of mass frame of the a+X system, yielding:

where again m_aX = m_a m_X / (m_a + m_X) is the reduced mass between a and X, v their relative velocity, N_a and N_X are the total number of particles, and we get a factor of v² from assuming isotropic motion and using spherical-polar coordinates in velocity space, so d³v = 4π v²dv.

The term in the integral is the average over the distribution of velocities of \langle\sigma(v)v\rangle, which has the dimension of [L²]×[L]/[t] = [L³]/[t].

N.B.: this does not depend on the density ρ!

Often, this average quantity is approximated as a powerlaw from the known value at a certain temperature T₀.

this is convenient because powerlaw dependencies are intuitive, and often people will quote the exponent β in arguments - be aware it is a big oversimplification: in reality β\equivβ(T) itself, but since σ(v) is non-zero only for a small range of T, taking β\simeq constant is not that problematic.

N.B.: thermonuclear reaction rates, because of the probability of tunneling through the Coulomb barrier are extremely sensitive to T. This makes nuclear physics equations very stiff numerically and can be a problem when computing stellar models.

Since velocity is not a great quantity to use in QM problems, and we have seen above that to get nuclear reactions we do need to account for QM effects, we can rewrite the rate above using that for r→ +∞ the relative energy between a and X is purely kinetic (the potential goes to zero faster than r^-2 because of electron screening!). Thus v=(2E/m_aX)^1/2, and we can also use that the cross section σ is proportional to the tunneling probability so σ(E)\propto exp(-b/E^1/2):

We can graph the part within the integral to understand where the rate is going to peak:

The product of an exponentially decreasing Maxwell-Boltzmann distribution \propto exp(-E) times the exponentially growing tunneling probability \propto exp(-1/E^1/2) will result in a very peaked integrand, the so called Gamow peak: even without nuclear resonances making the cross section σ peak (because E is within the width of a metastable energy level), the nuclear reaction rate is still very peaked around a specific energy!

Nuclear physics in stellar evolution codes

When trying to model the structure and evolution of a star, we cannot carry out all the integrals we wrote down here on the fly. Instead, we rely on tabulated nuclear reaction rates as a function of T and ρ.

This is a topic of active research, with certain reactions being particularly uncertain (e.g.,¹²C(α,γ)¹⁶O which determines the C/O ratio in the Universe, or even the 3α reaction that determines the formation of carbon in the first place). See for example Shen et al. 2023.

Composition changes and stellar evolution term

We already wrote Eq. \ref{eq:eps_nuc} for the energy generation term ε_nuc entering in the stellar structure equation describing energy conservation (dL/dm = ε_nuc - ε_ν + ε_grav). With a theory (and experimental data) to calculate the nuclear reaction rates r_i for each possible reaction, we have now a complete set of equations for the structure of a star at any given point in time under the assumption of spherical symmetry (which in nature can and is broken occasionally!).

However, because of nuclear reactions, the nuclei in the stellar plasma progressively change (on a timescale ∼ τ_nuc, which is long compared to all other timescales, but short compared to the age of the Universe!). This introduces the equations that drive the evolution of stars.

The change in number density of each nuclear species i per unit time is determined by its production rate through all possible nuclear reactions j+k→ i+… that have i as an outgoing particle, minus the destruction rate i+j→ something:

where the δ_i,j expresses that if i=j (reaction of two nuclei of the same species), then two nuclei are destroyed.

Using that n_i = X_iρ/(A_im_u) and what we have learned to determine the rates r_ij and r_kl we can rewrite this as a function of variables that already appear in our stellar structure equations:

and we have one such equation for each species i to consider. In case there is also mixing (for example due to convection), then we need to add to each of these equation a mixing term (advective or diffusive).

This effectively completes the set of equations we need to study not only the structure of stars but also their long-term evolution.

Appendix: Nuclear liquid drop model

In the 1930s Gamow, and later Bethe, Weiszacker and collaborators developed a model of the nuclear interactions based on the analogy with a liquid drop, which explains most of the observed features available at the time (e.g., the constant density of the nuclei, the preference for N=Z=A/2, the existence of particularly bound nuclei with certain even N and Z numbers).

N.B.: this is still a rough approximation that has been updated by other models since then, but paints a physically intuitive picture that has a wide range of applicability. Still today, most nuclear potential models are either specialized to a narrow energy range (so called local potentials), or involve a large number of experimentally determined parameters like this model.

In the liquid drop model, the nuclear binding energy can be expressed as:

where the coefficients a_i are experimentally derived. Each term has a specific interpretation, keeping in mind the empirical relation for the nuclear radius r_n ≡ r_n(A) \propto A^1/3:

• a_vol is the coefficient for a volume term, that is negative and expresses the fact that on short range nucleons attract (F_n=- ∇ V_n) each other and nuclei hold together.
• a_surf is a surface correction on the previous term, and expresses the fact that nucleons at the edge of the nucleus feel the interaction with fewer other nucleons.
• a_Coulomb expresses the electromagnetic repulsion between protons
• a_sym expresses the fact that stable nuclei prefer to have N=Z (unless N>Z is necessary to increase A without increasing the Coulomb term for heavy nuclei)
• a_coupling \propto ± A^-3/4 which is >0 if both N and Z are even, <0 if both are odd, and zero otherwise. This term expresses the fact that stable nuclei tend to prefer filling the energy levels for protons and neutrons (think by analogy with atoms wanting to fill their electron levels to be stable), so they tend to prefer having an even number for each, one spin up and one spin down in each energy level. (N.B.: this means that there are some number of nucleons for which nuclei are particularly stable, these are the “magic numbers” that will come back in the next lecture)