Notes on Thermodynamics

Derivatives With Respect to Composition

In the following we assume that the molar mass of species \(i\) is given by its atomic mass number, \(A_i = N_i + Z_i\) where \(N_i\) is the number of neutrons and \(Z_i\) is the number of protons in the isotope. This is a slight approximation that ignores the mass difference between protons and neutrons as well as some minor binding energy terms.

The number density [cm:math:^-{3}] of isotope \(i\) is can be formed from the mass density and the molar mass of that isotope as follows:

(324)\[\label{eq:number density} n_i = \frac{\rho_i N_\text{A}}{A_i},\]

where \(N_\text{A}\) is Avogadro’s number [# / mole]. The molar abundance, \(Y_i\), is a measure of the number of moles of species \(i\) per gram in the system:

(325)\[\label{eq:molar abundance} Y_i = \frac{n_i}{\rho N_\text{A}} = \frac{\rho_i}{\rho}\frac{1}{A_i} \equiv \frac{X_i}{A_i}\]

where we have defined the mass fraction, \(X_i = \frac{\rho_i}{\rho}\). Note

(326)\[\sum_i X_i = 1.\]

We write the average molar mass and average proton number as:

(327)\[\bar{A} = \frac{\sum_i A_i n_i}{\sum_i n_i} = \left(\sum_i X_i\right) \left(\sum_i \frac{X_i}{A_i}\right)^{-1}\]

(328)\[\bar{Z} = \frac{\sum_i Z_i n_i}{\sum_i n_i} = \left(\sum_i Z_i \frac{X_i}{A_i}\right)\left(\sum_i \frac{X_i}{A_i}\right)^{-1}.\]

Our algorithm requires terms involving the derivative thermodynamic variables (\(p\) or \(e\), e.g.) with respect to composition. Our EOS does not return such derivatives but instead returns derivatives of these variables with respect to \(\bar{A}\) and \(\bar{Z}\). Using the chain rule, we have

(329)\[\label{eq:p_xk} \frac{\ptl p}{\ptl X_k} = p_{X_k} = \frac{\ptl p}{\ptl \bar{A}}\frac{\ptl \bar{A}}{\ptl X_k} + \frac{\ptl p}{\ptl \bar{Z}}\frac{\ptl \bar{Z}}{\ptl X_k}.\]

From Eq.327 and Eq.328 we have

(330)\[\frac{\ptl \bar{A}}{\ptl X_k} = \left(\sum_i \frac{X_i}{A_i}\right)^{-1} - \frac{\bar{A}^2}{A_k} = -\frac{\bar{A}}{A_k}\left(\bar{A} - A_k\right)\]

(331)\[\frac{\ptl \bar{Z}}{\ptl X_k} = \left(\frac{Z_k}{A_k}\right)\left(\sum_i \frac{X_i}{A_i}\right)^{-1} - \frac{\bar{Z}}{A_k}\left(\sum_i \frac{X_i}{A_i}\right)^{-1} = -\frac{\bar{A}}{A_k}\left(\bar{Z} - Z_k\right),\]

where after differentiation we have used Eq.326 to write

(332)\[\left(\sum_i \frac{X_i}{A_i}\right)^{-1} = \bar{A}.\]

We therefore have

(333)\[\label{eq:p_Xk_full} p_{X_k} = -\frac{\bar{A}}{A_k}\left(\bar{A} - A_k\right) \frac{\ptl p}{\ptl\bar{A}} - \frac{\bar{A}}{A_k}\left(\bar{Z} - Z_k\right) \frac{\ptl p}{\ptl\bar{Z}}.\]

Before it was brought to our attention by Frank Timmes, we were missing the second term in Eq.330. The only place where such terms appear in our algorithm is in a sum over all species, such as:

(334)\[\sum_i p_{X_i}\dot{\omega}_i = -\bar{A}^2\frac{\ptl p}{\ptl\bar{A}}\sum_i \frac{\dot{\omega}_i}{A_i} +\bar{A}\frac{\ptl p}{\ptl\bar{A}}\sum_i \dot{\omega}_i -\bar{A}\bar{Z}\frac{\ptl p}{\ptl\bar{Z}}\sum_i \frac{\dot{\omega}_i}{A_i} +\bar{A}\frac{\ptl p}{\ptl \bar{Z}}\sum_i\frac{Z_i}{A_i}\dot{\omega}_i.\]

The second term in Eq.334 is identically zero because

(335)\[\sum_k \dot{\omega}_k \equiv 0.\]

This second term arises from what was added to Eq.330 by Frank’s correction. Therefore, although important for individual derivatives with respect to composition, this correction term has no effect on our solution.

Convective stability criterion

Here we look at the criterion for convective stability in the case of non-uniform chemical composition. This section follows Cox & Giuli :raw-latex:`\cite{cg-ed2}` closely (see chapter 13).

Consider a fluid parcel that gets displaced upwards (against gravity) from a radial location \(r\) to \(r + \Delta r\). The parcel is stable to convection if the displaced parcel’s density is greater than the surrounding fluid and gravity pushes the parcel back towards where it came from. Then the criterion for stability should be

(336)\[\begin{split}\begin{aligned} \rho_{parcel}(r+\Delta r) - \rho_{background}(r + \Delta r) &>& 0 \\ \bigg[\rho_{parcel}(r) + \bigg(\frac{d\rho}{dr}\bigg)_{parcel}\Delta r\bigg] - \bigg[\rho_{background}(r) + \bigg(\frac{d\rho}{dr}\bigg)_{background}\Delta r\bigg] &>& 0 \end{aligned}\end{split}\]

Since the parcel originates at r, \(\rho_{parcel}(r) = \rho_{background}(r)\) and so the stability criterion is

(337)\[\bigg(\frac{d\rho}{dr}\bigg)_{parcel} > \bigg(\frac{d\rho}{dr}\bigg)_{background}\]

Since the total pressure, \(P\), always increases inward in a star in hydrostatic equilibrium, we can use \(P\) instead of \(r\) as the independent radial variable. Then condition for stability can be written as

(338)\[\bigg( \frac{d \ln \rho}{d \ln P}\bigg )_{parcel} < \bigg(\frac{d \ln \rho}{d \ln P}\bigg)_{background}\]

Using the equation of state \(P = P( \rho, T, \bar{\mu})\), where \(\bar{\mu}\) is the average mass per molecule, we can write

(339)\[d \ln P = \frac{\partial \ln P}{\partial \ln \rho} \bigg |_{T, \bar{\mu}}d \ln \rho + \frac{\partial \ln P}{\partial \ln T} \bigg |_{\rho, \bar{\mu}} d \ln T + \frac{\partial \ln P}{\partial \ln \bar{\mu}}\bigg |_{\rho, T} d \ln \bar{\mu}\\]

For convenience we introduce

(340)\[\chi_{\rho} = \frac{\partial \ln P}{\partial \ln \rho}\bigg |_{T,\bar{\mu}} \qquad \chi_T = \frac{\partial \ln P}{\partial \ln T} \bigg |_{\rho,\bar{\mu}} \qquad \chi_{\bar{\mu}} = \frac{\partial \ln P}{\partial \ln \bar{\mu}} \bigg |_{\rho, T}\]

Then we can rearrange Eq.339 to get

(341)\[\frac{d \ln \rho}{\partial \ln P} = \frac{1}{\chi_\rho} - \frac{\chi_T}{\chi_\rho} \frac{d \ln T}{d \ln P}- \frac{\chi_{\bar{\mu}}}{\chi_\rho} \frac{d \ln \bar{\mu}}{d \ln P}\]

Then the general stability criterion is

(342)\[\bigg ( \frac{1}{\chi_\rho} - \frac{\chi_T}{\chi_\rho} \frac{d \ln T}{d \ln P}- \frac{\chi_{\bar{\mu}}}{\chi_\rho} \frac{d \ln \bar{\mu}}{d \ln P} \bigg )_{parcel} < \bigg ( \frac{1}{\chi_\rho} - \frac{\chi_T}{\chi_\rho} \frac{d \ln T}{d \ln P}- \frac{\chi_{\bar{\mu}}}{\chi_\rho} \frac{d \ln \bar{\mu}}{d \ln P} \bigg )_{background}\]

Here’s where various assumptions/simplifications get used.

If no assumptions are made, you can’t get any further than Eq.342. Even in view of an infinitesimally small initial perturbation, you can’t, in general, assume the \(\chi\)’s in parcel are the same as the \(\chi\)’s in the background. This applies in the case where nuclear reactions and/or ionization change the composition of the parcel. This case tends not to be of much interest for two reasons. Either composition effects get incorporated implicitly through assuming chemical equilibrium. Or both of these terms can be neglected in the rising parcel. This would be justified if the timescale for reactions is long compared to the convective timescale, and either the same is true for ionization or the fluid is fully ionized.
If we assume that \(\bar{\mu}\) remains constant in the parcel, then \(\frac{d \ln \bar{\mu}}{d \ln P}\) drops out for the parcel. In this case, we can assume, in view of the arbitrarily small initial perturbation of the parcel, that \(\chi_\rho\) and \(\chi_T\) to have the same values in the parcel as in the background. Then the stability criterion becomes

(343)\[\bigg ( \frac{d \ln T}{d \ln P} \bigg )_{parcel} > \bigg ( \frac{d \ln T}{d \ln P} + \frac{\chi_{\bar{\mu}}}{\chi_T} \frac{d \ln \bar{\mu}}{d \ln P} \bigg )_{background} \label{eqn:Ledoux}\]

The Ledoux stability criterion is obtained by assuming that the parcel moves adiabatically.
If we assume that the background is in chemical equilibrium and the parcel achieves instantaneous chemical equilibrium, then \(\bar{\mu} = \bar{\mu}(\rho,T)\) for the background and the parcel. (Note that we aren’t requiring constant composition in the parcel here.) The effect of variable composition are then absorbed into \(\chi_\rho\) and \(\chi_T\). Again, we can take \(\chi_\rho\) and \(\chi_T\) to have the same values in the parcel as in the background. The criterion then is

(344)\[\bigg ( \frac{d \ln T}{d \ln P} \bigg )_{parcel} > \bigg ( \frac{d \ln T}{d \ln P} \bigg )_{background} \label{eqn:Schwarz}\]

We obtain the Schwarzchild criterion for stability if we also assume the parcel moves adiabatically.

The Scharwzchild criterion can be recast in terms of entropy if the EOS is taken as \(P(\rho, S)\) instead of \(P(\rho, T)\). Then, in place of Eq.339 we have

(345)\[d \rho = \frac{\partial \rho}{\partial P} \bigg |_{S} d P + \frac{\partial \rho}{\partial S} \bigg |_{P} dS\]

We can substitute this into Eq.337 for stability, and assuming the parcel moves adiabatically, we get

(346)\[\bigg ( \frac{\partial \rho}{\partial S} \bigg |_{P} \frac{dS}{dr} \bigg )_{background}< 0\]

One of Maxwell’s relations is

(347)\[\frac{\partial \rho^{-1}}{\partial S} \bigg |_{P} = \frac{\partial T}{\partial P} \bigg |_{S}\]

All thermodynamically stable substances have temperatures that increase upon adiabatic compression, i.e. \(\frac{\partial T}{\partial P} \big |_{S} > 0\). So Maxwell’s relation implies that \(\frac{\partial \rho}{\partial S} \big |_{P} < 0\). The stability criterion then becomes

(348)\[\bigg ( \frac{d S}{d r} \bigg )_{background} > 0\]

Determining which stability criterion we want to enforce in creating the initial model is complicated by the phenomenon of semiconvection, which occurs when the Ledoux criterion is satisfied but the Schwarzchild is not, i.e.

(349)\[\bigg ( \frac{d \ln T}{d \ln P} \bigg )_{parcel} < \bigg ( \frac{d \ln T}{d \ln P} \bigg )_{background} < \bigg ( \frac{d \ln T}{d \ln P} \bigg )_{parcel} - \bigg ( \frac{\chi_{\bar{\mu}}}{\chi_T} \frac{d \ln \bar{\mu}}{d \ln P} \bigg )_{background}\]

(Note that \(\chi_{\bar{\mu}}\) is negative, as pressure is inversely proportional to mass per particle, and \(\frac{d \ln \bar{\mu}}{d \ln P}\) is positive, since nuclear reactions synthesize more massive particles in the center of the star.) In this case, when a rising parcel eventually reaches neutral buoyancy, it will have a temperature excess in comparison to it’s surroundings. If the parcel can retain it’s identity against diffusive mixing with the background long enough for significant heat exchange to occur, then the parcel’s temperature will drop, it will contract increasing it’s density, and the parcel will move inwards. The time scale of semiconvection is much longer than the timescale of traditional convection.

When we set up an initial model, we want to minimize any initial tendency towards convective motions, as we want these to be driven by the heating due to nuclear reactions, not the initial configuration we supply. Thus I think we want to guard against semiconvection as well as “traditional” convection by using the stability criterion

(350)\[\bigg ( \frac{d \ln T}{d \ln P} \bigg )_{parcel} = \frac{d \ln T}{d \ln P} \bigg |_{S,\bar{\mu}} > \bigg ( \frac{d \ln T}{d \ln P} \bigg )_{background}\]

Although this looks like the Schwarschild criterion (and, because I’m not entirely sure on vocabulary, it might even be called the Schwarzchild criterion), this does not simplify to Eq.348 because we need to keep the explicit \(\bar{\mu}\) dependence in the EOS.

The question of whether we’re in chemical equilibrium or not might be a moot point since our EOS (or any other part of the code) doesn’t enforce chemical equilibrium. Thus, even in the case of chemical equilbrium, we can’t in general drop the explicit \(\bar{\mu}\) dependence from our equations. If we wanted to do that, then we would need \(\bar{\mu}(\rho,T)\) to be substituted for \(\bar{\mu}\) inside the EOS.

Adiabatic Excess

The adiabatic excess, \(\Delta\nabla\), is a quantity used to determine if a system is stable (\(\Delta\nabla < 0\)) or unstable (\(\Delta\nabla > 0\)) to convection under the Schwarzschild criterion (i.e. neglecting compositional gradients). Cox and Giuli (see chapter 9) define three different “adiabatic exponents” that we will use:

(351)\[\begin{split}\begin{aligned} \Gamma_1 &\equiv& \left(\frac{d\ln p}{d\ln\rho}\right)_\text{ad} \\ \frac{\Gamma_2}{\Gamma_2-1} &\equiv& \left(\frac{d\ln p}{d\ln T}\right)_\text{ad} \\ \Gamma_3 - 1 &\equiv& \left(\frac{d\ln T}{d\ln\rho}\right)_\text{ad},\end{aligned}\end{split}\]

where the subscript “ad” means along an adiabat. We can combine the exponents to get the following relation

(352)\[\Gamma_1 = \left(\frac{\Gamma_2}{\Gamma_2-1}\right)\left(\Gamma_3-1\right).\]

The adiabatic excess is defined as

(353)\[\label{eq:adiabatic excess} \Delta\nabla = \nabla_\text{actual} - \nabla_\text{ad}\]

where

(354)\[\label{eq:thermal gradient} \nabla \equiv \frac{d\ln T}{d\ln P}\]

is the thermal gradient. It is important to note that these thermal gradients are only along the radial direction. The “actual” gradient can be found from finite differencing the data whereas the adiabatic term, \(\nabla_\text{ad} = \left(\Gamma_2-1\right) / \Gamma_2\), will need to be calculated at each point using thermodynamic relations. Our EOS only returns \(\Gamma_1\) so we need find another relation to use with Eq.352 to solve for the adiabatic excess.

The Schwarzschild criterion does not care about changes in composition and we therefore write \(p = p(\rho,T)\) and

(355)\[d\ln p = \chi_\rho d\ln\rho + \chi_T d\ln T\]

where

(356)\[\chi_\rho = \left(\frac{d\ln p}{d\ln\rho}\right)_T,\qquad \chi_T = \left(\frac{d\ln p}{d\ln T}\right)_\rho.\]

Dividing Eq.355 by \(d\ln\rho\) and taking this along an adiabat we have

(357)\[\left(\frac{d\ln p}{d\ln\rho}\right)_\text{ad} = \chi_\rho + \chi_T \left(\frac{d\ln T}{d\ln\rho}\right)_\text{ad}.\]

Using the \(\Gamma\)’s, we have

(358)\[\Gamma_1 = \chi_\rho + \chi_T\left(\Gamma_3-1\right).\]

Combining Eq.352 and Eq.358 to eliminate \(\Gamma_3\), we have:

(359)\[\label{eq:nabla_ad} \nabla_\text{ad} = \frac{\Gamma_1 - \chi_\rho}{\chi_T\Gamma_1}\]

which uses only terms which are easily returned from an EOS call.