\section{Derivatives With Respect to Composition}
\label{Sec: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.  \MarginPar{Note that this
  is common practice and that our EOS also makes this approximation.}  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$^-{3}$] of isotope $i$ is can be formed from the mass 
density and the molar mass of that isotope as follows:
\begin{equation}\label{eq:number density}
  n_i = \frac{\rho_i N_\text{A}}{A_i},
\end{equation}
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:
\begin{equation}\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}
\end{equation}
where we have defined the mass fraction, $X_i = \frac{\rho_i}{\rho}$.  Note
\begin{equation}\label{eq:mass fraction sums to 1}
  \sum_i X_i = 1.
\end{equation}

We write the average molar mass and average proton number as:
\begin{align}\label{eq:abar}
  \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}\\
  \label{eq:zbar}
  \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}.
\end{align}
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
\begin{equation}\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}.
\end{equation}
From \eqref{eq:abar} and \eqref{eq:zbar} we have
\begin{align}\label{eq:abar_X_k}
  \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)\\
  \label{eq:zbar_X_k}
  \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),
\end{align}
where after differentiation we have used \eqref{eq:mass fraction sums to 1}
to write
\[
\left(\sum_i \frac{X_i}{A_i}\right)^{-1} = \bar{A}.
\]
We therefore have
\begin{equation}\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}}.
\end{equation}
Before it was brought to our attention by Frank Timmes, we were missing the
second term in \eqref{eq:abar_X_k}.  The only place where such terms 
appear in our algorithm is in a sum over all species, such as:
\begin{equation}\label{eq:sum over species}
  \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.
\end{equation}
\MarginPar{Note that if we are dealing with symmetric nuclei 
  (as we have been up through paper III) where $Z_k = N_k = \half A_k$ for all
  $k$, then the last term in \eqref{eq:sum over species} is identically zero
  as well.}
The second term in \eqref{eq:sum over species} is identically zero because
\[
\sum_k \dot{\omega}_k \equiv 0.
\]
This second term arises from what was added to \eqref{eq:abar_X_k} by 
Frank's correction.  Therefore, although important for individual derivatives
with respect to composition, this correction term has no effect on our 
solution.

\section{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 
\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
\begin{eqnarray}
 \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{eqnarray}
Since the parcel originates at r, $\rho_{parcel}(r) = \rho_{background}(r)$ and
so the stability criterion is
% \bigg [ \bigg(\frac{d\rho}{dr}\bigg)_{parcel} - \bigg(\frac{d\rho}{dr}\bigg)_{background} \bigg ] \Delta r &>& 0 \\ 
\begin{equation}
 \bigg(\frac{d\rho}{dr}\bigg)_{parcel} > \bigg(\frac{d\rho}{dr}\bigg)_{background}
\label{eqn:basicStability}
\end{equation}

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
\[
 \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 
\begin{equation}
  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}\
\label{eqn:lnEOS}
\end{equation}
For convenience we introduce 
%\begin{eqnarray}
\[
  \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}
\]
%\end{eqnarray}
Then we can rearrange \ref{eqn:lnEOS} to get
\begin{equation}
  \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}
\end{equation}
%
Then the general stability criterion is
\begin{equation}
  \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}
\label{eqn:genStability}
\end{equation}
%
Here's where various assumptions/simplifications get used.  
\begin{enumerate}
\item
If no assumptions are made, you can't get any further than equation 
(\ref{eqn:genStability}).  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 ture for ionization or the 
fluid is fully ionized. 
\item
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
\begin{equation}
  \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}
\end{equation}
The Ledoux stability criterion is obtained by assuming that the parcel moves
adiabatically.
\item
If we assume   
that the background is in chemical equilibrium and the parcel acheives 
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 
\begin{equation}
\bigg ( \frac{d \ln T}{d \ln P} \bigg )_{parcel} > 
  \bigg ( \frac{d \ln T}{d \ln P} \bigg )_{background}
\label{eqn:Schwarz}
\end{equation}
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 equation (\ref{eqn:lnEOS}) we have 
\begin{equation}
 d \rho = \frac{\partial \rho}{\partial P} \bigg |_{S} d P + \frac{\partial \rho}{\partial S} \bigg |_{P} dS
\end{equation}
We can subsitute this into equation (\ref{eqn:basicStability}) for stability, 
and assuming the parcel moves adiabatically, we get 
\begin{equation}
  \bigg ( \frac{\partial \rho}{\partial S} \bigg |_{P} \frac{dS}{dr} 
  \bigg )_{background}< 0
\end{equation}
One of Maxwell's relations is 
\begin{equation}
\frac{\partial \rho^{-1}}{\partial S} \bigg |_{P} = \frac{\partial T}{\partial P} \bigg |_{S}
\end{equation}
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
\begin{equation}
 \bigg ( \frac{d S}{d r} \bigg )_{background} > 0 
\label{eqn:stabilityEntr}
\end{equation}
%
\end{enumerate}

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 statisfied but the Schwarzchild is not, 
i.e.
\begin{equation}
  \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}
\end{equation}
(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 sythesize 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 comparision to it's surroundings.  
%(Assuming ideal gas, plus black body radiation pressure if wanted.)  
% I think this qualifcation might be necessary, but want to think about
% degeneracy pressure too.
If the parcel can retain
it's identitiy against diffusive mixing with the background long enough for 
significant heat exhange 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 tendancy 
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
\begin{equation}
\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}
\end{equation}
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 equation (\ref{eqn:stabilityEntr})
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.


\section{Adiabatic Excess}
\label{Sec: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:
\begin{eqnarray*}
  \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{eqnarray*} 
where the subscript ``ad'' means along an adiabat.  We can combine the
exponents to get the following relation
\begin{equation}\label{eq:Gamma relations}
  \Gamma_1 = \left(\frac{\Gamma_2}{\Gamma_2-1}\right)\left(\Gamma_3-1\right).
\end{equation}

The adiabatic excess is defined as
\begin{equation}\label{eq:adiabatic excess}
  \Delta\nabla = \nabla_\text{actual} - \nabla_\text{ad}
\end{equation}
where
\begin{equation}\label{eq:thermal gradient}
  \nabla \equiv \frac{d\ln T}{d\ln P}
\end{equation}
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 \eqref{eq:Gamma relations} 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
\begin{equation}\label{eq:dp}
  d\ln p = \chi_\rho d\ln\rho + \chi_T d\ln T
\end{equation}
where
\[
\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 \eqref{eq:dp} by $d\ln\rho$ and taking this along an adiabat
we have
\begin{equation}\label{eq:dp2}
  \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}.
\end{equation}
Using the $\Gamma$'s, we have
\begin{equation}\label{eq:Gamma1 relation with Gamma2}
  \Gamma_1 = \chi_\rho + \chi_T\left(\Gamma_3-1\right).
\end{equation}
Combining \eqref{eq:Gamma relations} and \eqref{eq:Gamma1 relation
  with Gamma2} to eliminate $\Gamma_3$, we have:
\begin{equation}\label{eq:nabla_ad}
  \nabla_\text{ad} = \frac{\Gamma_1 - \chi_\rho}{\chi_T\Gamma_1}
\end{equation}
which uses only terms which are easily returned from an EOS call.


