*
Cococubed.com


Self-heating in reaction networks

Astronomy Research
   2024 Radiative Opacity
   2024 Neutrino Emission from Stars
   2023 White Dwarfs & 12C(α,γ)16O
   2023 MESA VI
   2022 Earendel, A Highly Magnified Star
   2022 Black Hole Mass Spectrum
   2021 Skye Equation of State
   2021 White Dwarf Pulsations & 22Ne
   Software Instruments

AAS Journals
   2024 AAS YouTube
   2024 AAS Peer Review Workshops

2024 ASU Energy in Everyday Life
2024 MESA Classroom
Outreach and Education Materials
    Solar Systems Astronomy
    Energy in Everyday Life
    Geometry of Art and Nature
    Calculus

Other Stuff:
   Bicycle Adventures
   Illustrations
   Presentations



Contact: F.X.Timmes
my one page vitae,
full vitae,
research statement, and
teaching statement.

$ \def\drvop#1{{\frac{d}{d{#1}}}} \def\drvf#1#2{{\frac{d{#1}}{d{#2}}}} \def\ddrvf#1#2{{\frac{d^2{#1}}{d{#2}^2}}} \def\partop#1{{\frac{\partial}{\partial {#1}}}} \def\ppartop#1{{\frac{\partial^2}{\partial {#1}^2}}} \def\partf#1#2{{\frac{\partial{#1}}{\partial{#2}}}} \def\ppartf#1#2{{\frac{\partial^2{#1}}{\partial{#2}^2}}} \def\mpartf#1#2#3{{\frac{\partial^2{#1}}{\partial{#2} \ {\partial{#3}}}}} $ A pdf of this note is avaliable.

Baryon number is an invariant. Define the abundance of species $Y_i$ by \begin{equation} Y_i = \frac{n_i}{n_B} = \frac{N_i}{N_B} \label{eq:y} \end{equation} where $N_i$ is the number of particles of isotope $i$, $N_B$ is the number of baryons, $n_i$ is the number density [cm$^{-3}$] of isotope $i$ and $n_B$ is baryon number density [cm$^{-3}$]. The number of baryons in isotope $i$ divided by the total number of baryons is the baryon fraction $X_i$, \begin{equation} X_i = Y_i \ A_i = \frac{n_i \ A_i}{n_B} \end{equation} where $A_i$ is the atomic mass number, the number of baryons in an isotope. Usually the baryon fraction is called the ``mass fraction''. Note \begin{equation} \sum X_i = \frac{n_B}{n_B} = 1 \label{eq:mconserv} \end{equation} is invariant under nuclear reactions. Define the baryon density, in atomic mass units, as \begin{equation} \rho_B = n_B \ m_u = \frac{n_B}{N_A} \hskip 0.2in {\rm g \ cm}^{-3} \end{equation} where $m_u$ is the atomic mass unit [g] and $N_A$ is the Avogadro number [g$^{-1}]$ in a system of units where the atomic mass unit is {\it defined} as 1/12 mass of an unbound atom of $^{12}$C is at rest and in its ground state.



The first law of thermodynamics is \begin{equation} dE = \delta Q -\delta W = \delta Q - PdV \hskip 0.2in {\rm erg} \label{eq:first} \end{equation} where $dE$ is the change in internal energy and $\delta Q$ is the heat added to the system. $\delta W = PdV$ is the mechanical work performed in expanding or contracting the system, where P [erg cm$^{-3}$] is the presure and $dV$ is the change in volume [cm$^3$]. If one chooses the entropy S [erg/K], volume $V$, and $N_i$ as the independent thermodynamic variables, then the total differential of the internal energy is \begin{align} dE & = \left ( \partf{E}{S}\right )_{V,N_i} dS + \left ( \partf{E}{V}\right )_{S,N_i} dV + \sum_i \left ( \partf{E}{N_i}\right )_{S,V} dN_i \notag \\[8pt] & = T dS - PdV + \sum_i \left ( \partf{E}{N_i}\right )_{S,V} dN_i \ . \label{eq:svn} \end{align} Comparing equations ($\ref{eq:first}$) and ($\ref{eq:svn}$), gives the familiar equality \begin{equation} \delta Q \equiv dE + PdV = T dS + \sum_i \left ( \partf{E}{N_i}\right )_{S,V} dN_i \label{eq:heat} \end{equation}
On the other hand, one may {\it choose} the temperature T, volume $V$, and $N_i$ as the independent thermodynamic variables. Then the total differential of the internal energy is \begin{align} dE & = \left ( \partf{E}{T}\right )_{V,N_i} dT + \left ( \partf{E}{V}\right )_{T,N_i} dV + \sum_i \left ( \partf{E}{N_i}\right )_{T,V} dN_i \notag \\ & = C_v dT - PdV + \sum_i \left ( \partf{E}{N_i}\right )_{T,V} dN_i \ , \label{eq:tvn} \end{align} where $C_v$ is the specific heat [erg/K] at constant volume. Comparing equations ($\ref{eq:first}$) and ($\ref{eq:tvn}$) gives the following equality \begin{equation} \delta Q \equiv dE + PdV = T dS + \sum_i \left ( \partf{E}{N_i}\right )_{S,V} dN_i = c_v dT + \sum_i \left ( \partf{E}{N_i}\right )_{T,V} dN_i \end{equation} Note the $\partial E / \partial N_i$ terms are sometimes called the chemical potential $\mu_i$. The heat added to the system from nuclear reactions is \begin{equation} \delta Q = E_{\rm nuc} - E_{\nu} \end{equation} which gives the equalities \begin{align} \delta Q & \equiv dE + PdV \notag \\ & = T dS + \sum_i \left ( \partf{E}{N_i}\right )_{S,V} dN_i \notag \\ & = C_v dT + \sum_i \left ( \partf{E}{N_i}\right )_{T,V} dN_i \notag \\ & = E_{\rm nuc} - E_{\nu} \end{align} Multiplying by the constant $N_A /N_B$ [1/g] gives the specific form of the first law: \begin{align} \delta q & \equiv de + P d \left(\frac{1}{\rho} \right) \notag \\ & = T ds + \sum_i \left ( \partf{e}{Y_i}\right )_{S,\rho} dY_i \notag \\ & = c_v dT + \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} dY_i \notag \\ & = \epsilon_{\rm nuc} - \epsilon_{\nu} \hskip 1.5in {\rm erg} \ {\rm g}^{-1} \end{align} Which could have been gotten by using the $(T,\rho,Y_i)$ Helmholtz free energy triplet as the independent variables. Operating with the time derivative $d/dt$ gives \begin{align} \drvf{e}{t} - \frac{P}{\rho^2} \drvf{\rho}{t} & = T \drvf{s}{t} + \sum_i \left ( \partf{e}{Y_i}\right )_{S,\rho} \drvf{Y_i}{t} \notag \\ & = c_v \drvf{T}{t} + \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} \notag \\ & = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} \hskip 1.5in {\rm erg} \ {\rm g}^{-1} \ {\rm s}^{-1} \label{eq:specific} \end{align} One may solve any combination of the expressions in equation ($\ref{eq:specific}$).



It is interesting to consider equation ($\ref{eq:specific}$) for an idea gas. It has been shown previously that the specific internal energy \begin{equation} e = \frac{3}{2} \frac{P}{\rho} = \frac{3}{2} \frac{N_A kT}{\overline{{\rm A}}} = \frac{3}{2} \ N_A kT \ \sum_{i} Y_i \hskip 0.2 in {\rm erg \ g}^{-1} \end{equation} has the derivative \begin{equation} \left ( \partf{e}{Y_i} \right )_{\rho, N_{i \ne j}} = \frac{3}{2} N_A kT \ . \end{equation} Its also been previously shown in these notes that nuclear energy generation rate is \begin{equation} \dot{\epsilon}_{\rm nuc} = - N_A c^2 \sum_{i} \drvf{Y_i}{t} \ m_i \hskip 0.2 in {\rm erg \ g}^{-1} \ {\rm s}^{-1} \ . \end{equation} Neglecting thermal neutrino losses, the temperature evolution version of equation ($\ref{eq:specific}$) becomes \begin{align} c_v \drvf{T}{t} & = \dot{\epsilon}_{\rm nuc} - \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} \notag \\[8pt] & = - N_A c^2 \sum_{i} \drvf{Y_i}{t} \ m_i - \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} \notag \\[8pt] & = - N_A c^2 \sum_{i} \drvf{Y_i}{t} \ m_i - \frac{3}{2} N_A kT \sum_i \drvf{Y_i}{t} \ . \label{eq:evol1} \end{align} Physically equation ($\ref{eq:evol1}$) says the change in the energy is the energy release from nuclear reactions minus the thermal energy of the now heavier (and fewer) gas particles. For example, say there are 4 protons. Each has an energy of 3/2 $kT$. When the 4 protons combine to form 1 helium nucleus, there is an energy gain from the nuclear reaction but now there is only 1 particle with an energy of 3/2 $kT$. The system gains energy because \begin{equation} \frac{m_{4He} - 4 m_p}{4m_p} \cdot c^2 \simeq 0.007 \cdot c^2 \simeq 7\times 10^{18} \gg N_A k T \simeq 8 \times 10^{14} \left ( \frac{T}{10^7 {\rm K}} \right) \hskip 0.2in {\rm erg} \ {\rm g}^{-1} \end{equation} for typical hydrogen burning temperatures. At the other extreme, the system looses energy when the second term in equation ($\ref{eq:evol1}$) becomes larger than the first term. At temperatures of $\simeq$ 10$^{10}$ K, this occurs when the relative mass change from nuclear burning becomes smaller than $\eta = \Delta m/m = N_A k T / c^2 \simeq 8 \times 10^{16} / 10^{21} \simeq 10^{-4}$, about a factor of 100 smaller than hydrogen burning's 0.007. In general, when \begin{equation} \frac{N_A k T}{\eta c^2 } \ll 1 \end{equation} the second term in equation ($\ref{eq:specific}$) can be safely dropped.



Let's consider a few specific instantiations of equation ($\ref{eq:specific}$).

For self-heating at constant density, $d\rho = 0$, and equation($\ref{eq:specific}$) reduces to \begin{equation} \drvf{e}{t} = c_v \drvf{T}{t} + \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} \ . \end{equation} is a convenient form. If one evolves the specific energy, one must do a root find on the equation of state to get the temperature for a given $(e, \rho, Y_i)$ triplet. If one evolves the temperature directly, the $\partial e / \partial Y_i$ terms must be included: \begin{equation} \drvf{T}{t} = {1 \over c_v} \left [ \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} - \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} \right ] \hskip 0.5in {\rm K} \ {\rm s}^{-1} \ . \label{eq:temp} \end{equation}


For self-heating at constant pressure, one approach is to solve equation ($\ref{eq:temp}$) and do a root find on the equation of state to get the density for a given $(T, P, Y_i)$ triplet - a Gibbs free energy basis set. {\it This is the easiest approach}. Alternatively, the hard way, expanding the total specific energy differential of equation ($\ref{eq:specific}$) in the $(T, \rho, Y_i)$ Helmholtz free energy basis gives \begin{align} & \drvf{e}{t} - {P \over \rho^2} \drvf{\rho}{t} = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} \notag \\[8pt] & \left ( \partf{e}{T}\right )_{\rho,Y_i} \drvf{T}{t} + \left ( \partf{e}{\rho}\right )_{T,Y_i} \drvf{\rho}{t} + \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} - {P \over \rho^2} \drvf{\rho}{t} = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} \notag \\[8pt] & \left ( \partf{e}{T}\right )_{\rho,Y_i} \drvf{T}{t} + \left [ \left ( \partf{e}{\rho}\right )_{T,Y_i} - {P \over \rho^2} \right ] \drvf{\rho}{t} = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} - \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} \ . \label{eq:step1} \end{align}

Now, the first law is an exact differential, which requires \begin{equation} P = \rho^2 \left ( \partf{e}{\rho} \right )_{T,Y_i} + T \left ( \partf{P}{T} \right )_{\rho, Y_i} \ . \label{eq:max} \end{equation} Using equation ($\ref{eq:max}$) to eliminate $(\partial e / \partial \rho )_{T,Y_i}$ in equation ($\ref{eq:step1}$) gives \begin{equation} c_v \drvf{T}{t} - {T \over \rho^2} \left ( \partf{P}{T}\right )_{\rho,Y_i} \drvf{\rho}{t} = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} - \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} \label{eq:step2} \end{equation}

In addition, for constant pressure, the total differential is $dP = 0$ \begin{equation} dP = \left ( \partf{P}{T}\right )_{\rho,Y_i} dT + \left ( \partf{P}{\rho}\right )_{T,Y_i} d\rho + \sum_i \left ( \partf{P}{Y_i}\right )_{T,\rho} dY_i = 0 \ . \end{equation} Solving for $d\rho$, \begin{equation} d\rho = - {1 \over \left ( \partf{P}{\rho}\right )_{T,Y_i} } \left [ \left ( \partf{P}{T}\right )_{\rho,Y_i} dT + \sum_i \left ( \partf{P}{Y_i}\right )_{T,\rho} dY_i \right ] \ . \label{eq:drho} \end{equation} Substituting equation ($\ref{eq:drho}$) into equation ($\ref{eq:step2}$) \begin{equation} c_v \drvf{T}{t} + {T \over \rho^2} \left ( \partf{P}{T}\right )_{\rho,Y_i} {\left ( \partf{\rho}{P}\right )_{T,Y_i} } \left [ \left ( \partf{P}{T}\right )_{\rho,Y_i} \drvf{T}{t} + \sum_i \left ( \partf{P}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} \right ] \notag \end{equation} \begin{equation} = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} - \sum_i \left ( \partf{e}{Y_i}\right )_{T,\rho} \drvf{Y_i}{t} \end{equation} Rearranging, \begin{equation} \left [ c_v + {T \over \rho^2} \left ( \partf{P}{T}\right )^2_{\rho,Y_i} \left ( \partf{\rho}{P}\right )_{T,Y_i} \right ] \drvf{T}{t} \notag \end{equation} \begin{equation} = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} - \sum_i \left [ \left ( \partf{e}{Y_i}\right )_{T,\rho} + {T \over \rho^2} \left ( \partf{P}{T}\right )_{\rho,Y_i} \left ( \partf{\rho}{P}\right )_{T,Y_i} \left ( \partf{P}{Y_i} \right )_{T,\rho} \right ]\drvf{Y_i}{t} \end{equation} using the so-called thermal and density exponents \begin{equation} \chi_t = {T \over P} \partf{P}{T} \hskip 1.0in \chi_{\rho} = {\rho \over P} \partf{P}{\rho} \end{equation} one has \begin{equation} \left [ c_v + {P \over \rho T} {\chi_T \over \chi_{\rho} } \right ] \drvf{T}{t} = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} - \sum_i \left [ \left ( \partf{e}{Y_i}\right )_{T,\rho} + \left ( \partf{\rho}{P}\right )_{T,Y_i} \left ( \partf{P}{Y_i} \right )_{T,\rho} \right ]\drvf{Y_i}{t} \end{equation} or \begin{equation} c_p \drvf{T}{t} = \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} - \sum_i \left [ \left ( \partf{e}{Y_i}\right )_{T,\rho} + {T \over \rho^2} \left ( \partf{P}{T}\right )_{\rho,Y_i} \left ( \partf{\rho}{P}\right )_{T,Y_i} \left ( \partf{P}{Y_i} \right )_{T,\rho} \right ]\drvf{Y_i}{t} \end{equation} where $c_p$ [erg g$^{-1}$ K$^{-1}$] is the specific heat at constant pressure. Hence, the final temperature evolution equation for a fixed pressure and changing composition: \begin{equation} \drvf{T}{t} = {1 \over c_p} \left [ \dot{\epsilon}_{\rm nuc} - \dot{\epsilon}_{\nu} - \sum_i \left [ \left ( \partf{e}{Y_i}\right )_{T,\rho} + {T \over \rho^2} \left ( \partf{P}{T}\right )_{\rho,Y_i} \left ( \partf{\rho}{P}\right )_{T,Y_i} \left ( \partf{P}{Y_i} \right )_{T,\rho} \right ]\drvf{Y_i}{t} \right ] \end{equation} Whew. Equation ($\ref{eq:temp}$) with a root find is simpler.



For self-heating at entropy, a model for a burning convective region, one approach is to solve equation ($\ref{eq:temp}$) and do a root find on the equation of state to get the density for a given $(T, S, Y_i)$ triplet. Alternatively, the hard way, is to follow the example above to get an explicit evolution equation for the entropy.