# Thermodynamic models

Now that our fundamental thermodynamics is refreshed and we understand that the minimisation needs a lot of thermodynamic derivatives, we're ready to develop consistent thermodynamic models to generate those derivatives.

## Generating functions and thermodynamic consistency

The thermodynamic potentials, when expressed as a function of their natural variables, provide a complete description of the state of a thermodynamic system. For example, consider the Gibbs free energy expressed in terms of its natural variables $G\left(T,\,p,\,\left\{N_i\right\}\right)$. If we know the values of the natural variables, we can calculate $G$. From the fundamental thermodynamics section we know the derivatives of $G$ are related as follows, \begin{align} {\rm d}G &= -S\,{\rm d}T + V\,{\rm d}p + \sum_i^{N_C}\mu_{i}\,{\rm d} N_{i}. \end{align} This implies that the derivatives of G in its natural variables are directly other thermodynamic state variables, \begin{align} \left(\frac{\partial G}{\partial T}\right)_{p,\left\{N_i\right\}} &= -S & \left(\frac{\partial G}{\partial p}\right)_{T,\left\{N_i\right\}} &= V & \left(\frac{\partial G}{\partial N_i}\right)_{T,p,\left\{N_{j\neq i}\right\}} &= \mu_i \end{align} From here we can derive all other thermodynamic potentials, for example, \begin{align} H &= G - T\,S & U &= H - p\,V & A &= G - p\,V. \end{align} We can take further derivatives in the natural variables to evaluate properties such as the heat capacity $C_p = \left(\partial H/\partial T\right)_{p,\left\{N_i\right\}}$. In this way, every interesting thermodynamic property can be created by only specifying the functional form of $G(T,\,p,\,\left\{N_i\right\})$ and then taking derivatives. This is known as a generating function approach, and this can be applied to any thermodynamic potential in its natural variables. Most commonly either the Gibbs or Helmholtz $A(T,\,V,\,\left\{N_i\right\})$ free energies are used as generating functions due to the convenience of their natural variables.

It is sometimes convenient to generate enthalpy rather than entropy from the Gibbs free energy. Consider a reversible process in a closed system, \begin{align*} {\rm d} \left(\frac{G}{T}\right) &= \frac{1}{T} {\rm d} G - \frac{G}{T^2}{\rm d}T\\ &= \frac{1}{T} {\rm d} G - \frac{H-TS}{T^2}{\rm d}T\\ &= \frac{1}{T} \left({\rm d} G + S{\rm d}T\right) - \frac{H}{T^2}{\rm d}T\\ &= \frac{V}{T} {\rm d}p - \frac{H}{T^2}{\rm d}T. \end{align*} This illustrates that the volume and enthalpy can also be directly obtained from the functional form of $G(T,p)/\,T$, by taking derivatives in its natural variables (while holding the other natural state variables, including $N_i$, constant).

Specifying the entire thermodynamic system using a single
thermodynamic potential guarantees *thermodynamic consistency*. This
is the requirement that all properties are correctly related via
their derivatives as required by the laws of thermodynamics. This
may be accidentally violated if the system is specified another way,
i.e., via any of the derivatives. The most common way thermodynamic
consistency was violated in the past was to combine polynomial
functions for heat capacity and density which cannot be integrated
into a consistent thermodynamic potential.

## Simplifying relationships for thermodynamic derivatives

A large number of thermodynamic derivatives are required to implement the minimisation. This section presents a number of useful expressions and approaches which allow us to interrelate various derivatives to reduce the number which must be implemented to a minimal amount.

First, generalised partial properties are introduced to demonstrate that all extensive properties can be expressed in terms of their derivatives in their extensive variables. The Bridgman tables provide a convenient method of expressing any derivative of the thermodynamic potentials or their variables in terms of just three "material derivatives". A "generalised" product rule is then introduced to demonstrate how derivatives may be interchanged, particularly for the triple product rule. Finally, a key relationship between the partial molar and partial "volar" properties is derived.

### Partial properties

As demonstrated in the section on solution for the internal energy, functions which are first-order homogeneous in their variables can be immediately expressed in terms of these derivatives. This is useful as it provided a relationship between the internal energy and the other thermodynamic potentials; however, a generalised rule would eliminate the need to generate expressions for thermodynamic properties IF their derivatives are available.

Unfortunately, the two variable sets considered here contain both
homogeneous first-order ($\left\{N_i\right\}$, $V$) and
inhomogeneous ($T$, $p$) variables. In this case, Euler's method
does not extend to expressions which are functions of both; However,
a similar solution can be derived for these expressions provided the
inhomogeneous variables are restricted to intensive
properties *and* the extensive variables together uniquely
specify the total size of the system.

The derivation presented here is a generalisation of the proof for
partial *molar* properties which you can find in any
thermodynamic text (e.g., Smith, van Ness, and Abbott, Sec.11.2, 7th
Ed). Consider an extensive thermodynamic property, $M$, which is a
function of extensive $\left\{X_i\right\}$ and intensive
$\left\{y_i\right\}$ variables. The total differential is as
follows,
\begin{align}
{\rm d} M =
\sum_i \left(\frac{\partial M}{\partial X_i}\right)_{\left\{X_{i\neq j}\right\},\left\{y_{k}\right\}} {\rm d}X_i
+\sum_i \left(\frac{\partial M}{\partial y_i}\right)_{\left\{X_{k}\right\},\left\{y_{j\neq i}\right\}} {\rm d}y_i.
\end{align}
The extensive property $M$ can be converted to a corresponding
intensive property, $m$, by dividing by the system amount, i.e.,
$M=m\,N$. If the extensive properties are held constant and *
they are sufficient to determine the total system amount*, $N$,
then the system size, $N$, is also constant and may be factored out
of $M$ in the intensive partial differential terms.
\begin{align}
{\rm d} M =
\sum_i \left(\frac{\partial M}{\partial X_i}\right)_{\left\{X_{i\neq j}\right\},\left\{y_{k}\right\}} {\rm d}X_i
+N \sum_i \left(\frac{\partial m}{\partial y_i}\right)_{\left\{X_{k}\right\},\left\{y_{j\neq i}\right\}} {\rm d}y_i.
\end{align}
In addition, ${\rm d} M={\rm d} (N\,m) = m\,{\rm d} N+N\,{\rm
d}m$ and ${\rm d} X_i={\rm d} (N\,x_i) = x_i\,{\rm d} N+N\,{\rm
d}x_i$. Inserting these and factoring out the terms in $N$ and
${\rm d}N$ yields,
\begin{align}\label{eq:pmolarintermediate}
\left(m - \sum_i \left(\frac{\partial M}{\partial X_i}\right)_{\left\{X_{i\neq j}\right\},\left\{y_{k}\right\}} x_i\right){\rm d} N
+ N\left({\rm d}m - \sum_i \left(\frac{\partial M}{\partial X_i}\right)_{\left\{X_{i\neq j}\right\},\left\{y_{k}\right\}} {\rm d}x_i - \sum_i \left(\frac{\partial m}{\partial y_i}\right)_{\left\{X_{k}\right\},\left\{y_{j\neq i}\right\}} {\rm d}y_i\right) =
0.
\end{align}
As ${\rm d}N$ and $N$ can vary independently this equation is only
satisfied if the terms in parenthesis are each zero. Multiplying the
first term in parenthesis by $N$ and setting it equal to zero yields
the required identity,
\begin{align}\label{eq:partialProperty}
M = \sum_i \bar{\bar{m}}_i\,X_i.
\end{align}
where $\bar{\bar{m}}_i=\left(\partial M/\partial
X_i\right)_{\left\{X_{i\neq j}\right\},\left\{y_{k}\right\}}$ is a
partial property.
This is a generalised solution for an extensive property in terms of
its partial properties which are its derivatives in its extensive
variables and is the primary objective of our derivation; however,
an additional important equation for the partial properties can be
easily obtained. The derivative of Eq.\eqref{eq:partialProperty}
(first divided by $N$) is,
\begin{align}
{\rm d}m = \sum_i x_i\,{\rm d}\bar{\bar{m}}_i + \sum_i \bar{\bar{m}}_i\,{\rm d}x_i.
\end{align}
This can be used to eliminate ${\rm d}m$ from the right term of
Eq.\eqref{eq:pmolarintermediate} which is also set equal to zero
(and multiplied by N) to give,
\begin{align}\label{eq:GibbsDuhem}
\sum_i X_i\,{\rm d}\bar{\bar{m}}_i
- \sum_i \left(\frac{\partial M}{\partial y_i}\right)_{\left\{X_{k}\right\},\left\{y_{j\neq i}\right\}} {\rm d}y_i = 0.
\end{align}
This is the generalised **Gibbs/Duhem** equation which
interrelates the intensive properties of the system to the partial
properties.

The most well known applications of Eqs.\eqref{eq:partialProperty}
and \eqref{eq:GibbsDuhem} are for the partial *molar*
properties when $p$, $T$, and $\left\{N_i\right\}$ are the state
variables. In this case Eq.\eqref{eq:partialProperty} is
\begin{align}\label{eq:partialMolarProperty}
M = \sum_i \bar{m}_i\,X_i.
\end{align}
where $\bar{m}_i=\left(\partial M/\partial
N_i\right)_{p,T,\left\{N_{j\neq i}\right\}}$ is the partial molar
property. The most important partial molar property is the chemical
potential,
\begin{align*}
\mu_i = \bar{g}_i = \left(\frac{\partial G}{\partial N_i}\right)_{T,p,\left\{N_{j\neq i}\right\}},
\end{align*}
and it has already been proven that
Eq.\eqref{eq:partialMolarProperty} applies in this case, i.e.,
Eq.\eqref{eq:Grule} gives $G=\sum_iN_i\,\mu_i$. The corresponding
Gibbs-Duhem equation for the chemical potential in
$T,p,\left\{N_i\right\}$ is the most well-known form,
\begin{align}\label{eq:GibbsDuhemTpMu}\nonumber
\sum_i N_i\,{\rm d}\mu_i
- \left(\frac{\partial G}{\partial T}\right)_{p,\left\{N_i\right\}} {\rm d}T
- \left(\frac{\partial G}{\partial p}\right)_{T,\left\{N_i\right\}} {\rm d}p
\\
= \sum_i N_i\,{\rm d}\mu_i
+ S\,{\rm d}T
- V\,{\rm d}p
= 0.
\end{align}

As derivatives are distributive, and $T$ and $p$ are held constant, the partial molar properties for the thermodynamic potentials satisfy similar relationships as the original potentials (Eqs.\eqref{eq:Urule}-\eqref{eq:Grule}). \begin{align}\label{eq:partialmolarrelationstart} \mu_i&= \bar{h}_i - T\,\bar{s}_i\\ \bar{h}_i&= \bar{u}_i + p\,\bar{v}_i\\ \bar{a}_i&= \bar{u}_i - T\,\bar{s}_i\\ \bar{u}_i&= T\,\bar{s}_i - p\,\bar{v}_i + \mu_i \label{eq:partialmolarrelationend} \end{align} Thus, the partial molar properties not only provide a derivative in the molar amounts, but also completely describe all thermodynamic potentials and many extensive quantities (such as $V$). Simcem therefore only requires the implementation of expressions for three partial molar amounts ($\mu_i$, $\bar{v}_i$, $\bar{h}_i$) to specify all the thermodynamic potentials and their first derivative in the molar amounts for the $\left\{T,p,\left\{N_i\right\}\right\}$ variable set.

For the $\left\{T,V,\left\{N_i\right\}\right\}$ variable set, Eq.\eqref{eq:partialProperty} is, \begin{align}\label{eq:partialVolarProperty} M = \sum_i \breve{m}_i\, N_i + \left(\frac{\partial M}{\partial V}\right)_{T,\left\{N_{i}\right\}}\,V \end{align} where $\breve{m}_i = \left(\partial M/\partial N_i\right)_{T,V,\left\{N_{j\neq i}\right\}}$ is deemed here to be a partial "volar" quantity. A more useful form of this expression is Eq.\eqref{eq:molartovolarproperty} which is derived in the later section on the transformation between $\bar{m}_i$ and $\breve{m}_i$.

### Material properties (Bridgman tables)

There are a number of interesting properties which are derivatives of the thermodynamic potentials. For example, consider the isobaric heat capacity, \begin{align*} C_p = \left(\frac{\partial H}{\partial T}\right)_{p} = -T \left(\frac{\partial^2 G}{\partial T^2}\right)_{p} \end{align*}

Using partial molar properties allows us to illustrate a complication in the calculation of these material properties for multi-phases/sub-systems. For example, expressing the (extensive) heat capacity in terms of the partial molar (AKA specific) heat capacity using Eq.\eqref{eq:partialMolarProperty} yields the following expression, \begin{align}\label{eq:mixCp} C_p = \left(\frac{\partial H}{\partial T}\right)_{p}= \left(\frac{\partial}{\partial T} \sum_{i} N_{i}\, \bar{h}_i\right)_{p} = \sum_i\left[\bar{h}_i\left(\frac{\partial N_i}{\partial T}\right)_{p} +N_i\, \bar{c}_{p,i}\right] \end{align} where $\bar{c}_{p,i}=\left(\partial \bar{h}_i/\partial T\right)_{p}$ is the partial molar isobaric heat capacity. The term on the LHS of Eq.\eqref{eq:mixCp} arises from changes to the enthalpy caused by species transferring in and out of the system as the equilibrium state changes. This results in additional contributions to the apparent heat capacity above the partial molar isobaric heat capacity. For example, when a single-component fluid is at its boiling point, the apparent heat capacity of the overall system is infinite as $\partial N_i/\partial T$ is infinite due to the discontinuous change from liquid to vapour causing the instananeous transfer of molecules from one phase to another.

To complement the "equilibrium" thermodynamic $C_p$ above, it is convenient to define "frozen" thermodynamic derivatives where there are no molar fluxes, i.e., the "frozen" isobaric heat capacity, \begin{align}\label{eq:frozenCp} C_{p,\left\{N_j\right\}^{N_c}} = \sum_i^{N_C}N_i\,\bar{c}_{p,i} \end{align} The "frozen" properties are required while calculating the gradient of thermodynamic potentials during minimisation, and arise as all molar quantities are held constant while these derivatives are taken.

The $C_p$ is just one material property; however, there are many
other thermodynamic derivatives which may be
calculated. Fortunately,
the Bridgman
tables is a convenient method to express any material property
as a function of just three key material derivatives. The heat
capacity is one and the other two are the isothermal ($\beta_T$) and
isobaric ($\alpha$) expansivities,
\begin{align*}
\alpha\,V=\left(\frac{\partial V}{\partial T}\right)_p &= \sum_i\left[\bar{v}_i\left(\frac{\partial N_i}{\partial T}\right)_{p,\left\{N_{j\neq i}\right\}} + N_i\,\bar{\alpha}_i\,\bar{v}_i\right]
\\
\beta_T\,V=-\left(\frac{\partial V}{\partial p}\right)_T &= \sum_i\left[-\bar{v}_i\left(\frac{\partial N_i}{\partial p}\right)_{T,\left\{N_{j\neq i}\right\}} + N_i\,\bar{\beta}_{T,i}\,\bar{v}_i\right]
\end{align*}
where $\bar{\alpha}_i\,\bar{v}_i= \left(\partial
\bar{v}_{i}/\partial T\right)_{p,\left\{N_j\right\}}$ and
$\bar{\beta}_{T,i}\,\bar{v}_i=-\left(\partial v_i/\partial
p\right)_{T,\left\{N_j\right\}}$. Again, "frozen" material
derivatives are available,
\begin{align}\label{eq:frozenAlphaBeta}
\left(\alpha\,V\right)_{\left\{N_i\right\}}&= \sum_i^{N_C}N_i\,\bar{\alpha}_i\,\bar{v}_i & \left(\beta_T\,V\right)_{\left\{N_i\right\}}&= \sum_i^{N_C}N_i\,\bar{\beta}_i\,\bar{v}_i
\end{align}
The terms $\left(\partial N_i/\partial T\right)_{p,\left\{N_{j\neq
i}\right\}}$ and $\left(\partial N_i/\partial
p\right)_{T,\left\{N_{j\neq i}\right\}}$ which appear in the
material properties must be determined from the solution to the
thermodynamic minimisation problem. They quantify how
the *equilibrium* molar amounts change for a variation in
temperature and pressure and thus must account for the constraints
placed on the equilibrium state and the movement of the minimia.

The Bridgman table approach decomposes every thermodynamic derivative into a look-up table for the numerator and denominator expressed in terms of the three material derivatives, $C_p$, $\alpha$, and $\beta_T$. For example, consider the following "unnatural" derivative, \begin{align} \left(\frac{\partial G}{\partial V}\right)_{T} = \frac{\left(\partial G\right)_{T}}{\left(\partial V\right)_{T}} = \frac{-V}{\beta_T\,V} = -\beta_T^{-1} \end{align} Thus, to generate any derivative required for minimisation, only the three "material" derivatives and three partial molar properties are required.

### Product rules

This section is almost directly copied from this math.StackExchange.com post.

Consider any expression for a thermodynamic property, $M$, written in terms of the state variables, $M=M\left(\vec{X}\right)$. It is straightforward to transform this function into an implicit function $F=F\left(M,\,\vec{X}\right)=M(\vec{X}) - M=0$ which helps to illustrate that $M$ can be treated as a variable on an equal basis as $\vec{X}$. To allow a uniform representation, the arguments of $F$ are relabeled such that $F(x_1,x_2,\ldots,x_N)=0$. As the function $F$ remains at zero, holding all variables except two constant taking the total derivative and setting ${\rm d}F=0$ yields the following identity, \begin{align*} \left(\frac{\partial x_i}{\partial x_j}\right)_{x_k\forall k\neq i,j} &= - \frac{\left(\partial F/\partial x_j\right)_{x_k\forall k\neq j}}{\left(\partial F/\partial x_i\right)_{x_k\forall k\neq i}} & \text{for $i\neq j$}. \end{align*} This rule can be combined repeatedly and the terms on the RHS eliminated provided the variables loop back on themselves. For example, \begin{align*} \left(\frac{\partial x_1}{\partial x_2}\right)_{x_k\forall k\neq 1,2}\left(\frac{\partial x_2}{\partial x_3}\right)_{x_k\forall k\neq 2,3}\ldots\left(\frac{\partial x_n}{\partial x_1}\right)_{x_k\forall k\neq n,1} &= (-1)^n \cancelto{1}{\frac{\partial F/\partial x_2}{\partial F/\partial x_1}\frac{\partial F/\partial x_3}{\partial F/\partial x_2}\ldots\frac{\partial F/\partial x_1}{\partial F/\partial x_n}} & \text{for $n > 1$}. \end{align*} where the variables held constant were dropped for brevity on the right hand side. The first value $n=2$ yields the well known expression, \begin{align} \left(\frac{\partial x_1}{\partial x_2}\right)_{x_k\forall k\neq 1,2} \left(\frac{\partial x_2}{\partial x_1}\right)_{x_k\forall k\neq 1,2} &= 1, \end{align} or, more familiarly (and without the explicit variables held constant), \begin{align} \frac{\partial x_1}{\partial x_2} &= \left(\frac{\partial x_2}{\partial x_1}\right)^{-1}. \end{align} Finally, $n=3$ yields the triple product rule, \begin{align} \left(\frac{\partial x_1}{\partial x_2}\right)_{x_3}\left(\frac{\partial x_2}{\partial x_3}\right)_{x_1}\left(\frac{\partial x_3}{\partial x_1}\right)_{x_2} &= -1 \end{align} where it is implied that all other variables other than $x_1$, $x_2$, and $x_3$ are held constant. This rule has wide application but is particularly attractive when derivatives hold an "awkward" thermodynamic property constant. For example, consider the following case \begin{align} \left(\frac{\partial p}{\partial T}\right)_G = -\left(\frac{\partial G}{\partial p}\right)_T^{-1}\left(\frac{\partial G}{\partial T}\right)_p = \frac{S}{V} \end{align} The LHS is difficult to directly calculate in the state variables chosen here; however, the RHS arising from the triple product rule is in terms of natural derivatives of $G$ which are straightforward to calculate. The triple product rule is used extensively in the following sections to express complex expressions in terms of natural derivatives.

### Relation between $\bar{x}_i$ and $\breve{x}_i$

Consider a property, $M$, which is a function for four variables $x_1,x_2,x_3,x_4$. The total derivative is, \begin{align*} {\rm d}M = \sum_i \left(\frac{\partial M}{\partial x_i}\right)_{x_{i\neq j}} {\rm d}x_i \end{align*} Holding two variables constant and taking a derivative wrt a third yields, \begin{align*} \left(\frac{\partial M}{\partial x_1}\right)_{x_2,x_3} = \left(\frac{\partial M}{\partial x_1}\right)_{x_2,x_3,x_4} + \left(\frac{\partial M}{\partial x_4}\right)_{x_1,x_2,x_3} \left(\frac{\partial x_4}{\partial x_1}\right)_{x_2,x_3} \end{align*} Three relablings of this expression can be used to interrelate two derivatives which only differ in one variable which is held constant. \begin{align*} \left(\frac{\partial M}{\partial x_1}\right)_{x_2,x_3} = \left(\frac{\partial M}{\partial x_1}\right)_{x_2,x_4} + \left(\frac{\partial M}{\partial x_4}\right)_{x_1,x_2} \left(\frac{\partial x_4}{\partial x_1}\right)_{x_2,x_3} - \left(\frac{\partial M}{\partial x_3}\right)_{x_2,x_3,x_4} \cancelto{0}{\left(\left(\frac{\partial x_3}{\partial x_4}\right)_{x_1,x_2}\left(\frac{\partial x_4}{\partial x_1}\right)_{x_2,x_3}+\left(\frac{\partial x_3}{\partial x_1}\right)_{x_2,x_4}\right)} \end{align*} where the final term in parenthesis is cancelled to zero using the triple product rule. \begin{align} \left(\frac{\partial M}{\partial x_1}\right)_{x_2,x_3} = \left(\frac{\partial M}{\partial x_1}\right)_{x_2,x_4} + \left(\frac{\partial M}{\partial x_4}\right)_{x_1,x_2} \left(\frac{\partial x_4}{\partial x_1}\right)_{x_2,x_3} \end{align} This equation is particularly useful for changing between partial quantities while pressure or volume is held constant. For example, \begin{align} \left(\frac{\partial M}{\partial N_i}\right)_{T,p,\left\{N_{j\neq i}\right\}} = \left(\frac{\partial M}{\partial N_i}\right)_{T,V,\left\{N_{j\neq i}\right\}} + \left(\frac{\partial M}{\partial V}\right)_{T,\left\{N_{j}\right\}} \left(\frac{\partial V}{\partial N_i}\right)_{T,p,\left\{N_{j\neq i}\right\}} \end{align} or, in the notation used so far, \begin{align}\label{eq:molartovolarproperty} \bar{m}_i = \breve{m}_i + \bar{v}_i \left(\frac{\partial M}{\partial V}\right)_{T,\left\{N_{j}\right\}}. \end{align} This is a more useful form of Eq.\eqref{eq:partialVolarProperty}. The partial molar volume is inconvenient to derive directly when volume is an explicit variable; however, it may be expressed more conveniently using the triple product rule, \begin{align}\label{eq:TVpartialmolarvolume} \bar{v}_i = - \left(\frac{\partial p}{\partial n_i}\right)_{T,V,\left\{N_{j\neq i}\right\}} \left(\frac{\partial p}{\partial V}\right)_{T,\left\{N_{j}\right\}}^{-1}. \end{align} This allows us to obtain partial molar properties conveniently when working with volume as a state variable.

## Gibbs Models ($p$ as a variable)

In this section, the calculation of the required properties for minimisation with phases specified by the following set of variables is considered, \begin{align*} \vec{X}_\alpha &= \left\{T_\alpha,\,p_\alpha,\,\left\{N_{i}\right\}^{N_{C,\alpha}}\right\} . \end{align*} In this particular variable set, the required constraint derivatives to implement the derivative of the Lagrangian are as follows, \begin{align} \left(\frac{\partial V_\alpha}{\partial T}\right)_{p,\left\{N_{i}\right\}} &= \left(\alpha\,V\right)_{\left\{N_{i}\right\}}, \\ \left(\frac{\partial V_\alpha}{\partial p}\right)_{T,\left\{N_{i}\right\}} &= -\left(\beta_T\,V\right)_{\left\{N_{i}\right\}}, \\ \left(\frac{\partial V_\alpha}{\partial N_i}\right)_{T,p,\left\{N_{j\neq i}\right\}} &= \bar{v}_i \end{align} The thermodynamic potential derivatives required to specify the derivatives of the Lagrangian (Eq.\eqref{eq:Fderiv}) as generated from Eqs.\eqref{eq:GradEqsStart}-\eqref{eq:GradEqsEnd} are \begin{align} \left(\frac{\partial f}{\partial T}\right)_{p,\left\{N_i\right\}},\, \left(\frac{\partial f}{\partial p}\right)_{T,\left\{N_i\right\}},\, \left(\frac{\partial f}{\partial N_i}\right)_{T,p,\left\{N_{i\neq j}\right\}}, \end{align} where $f=\left\{H,G,-S,U,A\right\}$. These derivatives are easily expressed using the Bridgman tables in terms of the three standard material derivatives and the results are given in the table below.

$Y_1$ | $Y_2$ | $f$ | $\left(\frac{\partial f}{\partial N_i}\right)_{T,p,\left\{N_{j\neq i}\right\}}$ | $\left(\frac{\partial f}{\partial T}\right)_{p,\left\{N_j\right\}}$ | $\left(\frac{\partial f}{\partial p}\right)_{T,\left\{N_j\right\}}$ |
---|---|---|---|---|---|

$p$ | $T$ | $G$ | $\mu_i$ | $-S$ | $V$ |

$V$ | $T$ | $A$ | $\bar{a}_i$ | $-(S+p\,\alpha\,V)$ | $p\left(\beta_{T}\,V\right)_{\left\{N_i\right\}}$ |

$p$ | $S$ | $H$ | $\bar{h}_i$ | $C_{p,\left\{N_j\right\}}$ | $V(1-T\,\alpha)$ |

$p$ $V$ | $H$ $U$ | $-S$ | $-\bar{s}_i$ | $-T^{-1}\,C_{p,\left\{N_j\right\}}$ | $\left(\alpha\,V\right)_{\left\{N_i\right\}}$ |

$V$ | $S$ | $U$ | $\bar{u}_i$ | $C_{p,\left\{N_j\right\}}-p\left(\alpha\,V\right)_{\left\{N_i\right\}}$ | $p\left(\beta_{T}\,V\right)_{\left\{N_i\right\}}-T\,\left(\alpha\,V\right)_{\left\{N_i\right\}}$ |

In summary, a model using this variable set must provide implementations of $\mu_i$, $\bar{v}_i$, $\bar{s}_i$, $\bar{\alpha}_i$, $\bar{\beta}_{T,i}$, and $\bar{c}_{p,i}$. These are all straightforward to obtain by performing derivatives of a Gibbs free energy function in its natural variables or integration and differentiation if a mechanical equation of state, $V\left(T,p,\left\{N_i\right\}\right)$, is available.

All other partial molar properties are obtained using Eqs.\eqref{eq:partialmolarrelationstart}-\eqref{eq:partialmolarrelationend} expressed in the following form. \begin{align} \bar{a}_i&= \mu_i - p\,\bar{v}_i\\ \bar{h}_i&= \mu_i + T\,\bar{s}_i\\ \bar{u}_i&= \mu_i + T\,\bar{s}_i - p\,\bar{v}_i \end{align} Most other relevant thermodynamic properties are calculated using the Bridgman tables.

### Ideal gas model

The ideal gas model is not only a good approximation for gases at low pressures, but it is also a key reference "state" for most complex equations. A thermodynamic path can be constructed to the ideal gas state at low pressures for most models, thus in this case a "ideal-gas" contribution can be factored out from these models.

The ideal gas chemical potential is defined as follows, \begin{align}\label{eq:idealgasmuTP} \mu^{(ig)}_{i}(T,\,p) &= \mu_{i}^{(ig)}(T) + R\,T\ln \left(\frac{p}{p_0}\right) + R\,T\ln \left(\frac{N_{i}}{N_\alpha}\right). \end{align} where $p_0$ is the reference pressure at which the temperature-dependent term, $\mu^{(ig)}_{i}(T)$, is measured and $N_\alpha$ is the total moles of the phase $\alpha$.

Using the chemical potential as a generating function, the minimial set of partial molar properties is derived below. \begin{align*} \bar{v}_{i}^{(ig)} &= R\,T/p & \bar{s}_{i}^{(ig)} &= -\frac{\partial \mu_{i}^{(ig)}(T)}{\partial T} - R \ln\left(\frac{p}{p_0}\right) - R \ln\left(\frac{N_{i}}{N_\alpha}\right) \\ \bar{c}_{p,i}^{(ig)} &= -T\frac{\partial^2 \mu^{(ig)}_{i}(T)}{\partial T^2} & \alpha^{(ig)} &= T^{-1} \\ \beta_{T}^{(ig)} &= p^{-1} & & \end{align*} where the partial molar heat capacity has been implicitly defined as $C_{p,\left\{N_i\right\}}^{(ig)} = \sum_i N_i\,\bar{c}_{p,i}^{(ig)}$.

The model is not completely specified until the pure function of temperature, $\mu_{i}^{(ig)}(T)=h_{i}^{(ig)}(T) - T\,s_{i}^{(ig)}(T)$, is specified, and this is discussed in the following section.

This is not a proof of how mixing entropy arises (which is discussed below), but a demonstration that mixing entropy is consistent with Dalton's law holding true. Consider a single component ideal gas, we have, \begin{align*} g_{i,pure} &= g_i^0(T) + R\,T\ln \left(\frac{p}{p_0}\right) \end{align*} If Dalton's law is true, the pressure of a single species is reduced when mixed as it only exerts a partial pressure, ($p_i=x_i\,p$). Assuming the total pressure and temperature remains constant during mixing, \begin{align*} g_{i,mixed} &= g^0_i(T) + R\,T\ln \left(\frac{p_i}{p_0}\right) = g^0_i(T) + R\,T\ln \left(\frac{p}{p_0}\right) + R\,T\ln x_i\\ g_{i,mixed} - g_{i,pure} &= R\,T\ln x_i. \end{align*} As $x\le1$, $g_{mixed} - g_{pure}\le0$, which is consistent that two ideal gases will be mixed at equilibrium under constant $T,p$. Thus, Dalton's law is consistent with mixing entropy

### Mixing entropy from lattice arguments

First, it is noted that entropy is extensive. The entropy of
two separate systems is the sum of each,
\begin{align*}
S = S_1 + S_2
\end{align*}
Fundamentally, entropy is a measure of the number of states of
a system; however, the number of states
has *multiplicative* scaling. For example, consider a
collection of light switches. A single light switch has two
states but a collection of $N$ light switches has $2^N$
possible states. The number of states is clearly
multiplicative and yet entropy is additive. The entropy must
therefore be the logarithm of the number of states as $\ln 2^N
= N \ln 2$. This line of arguing leads to what is known as
Boltzmann's entropy,
\begin{align*}
S = R \ln W
\end{align*}
where $W$ is the number of states and $R$ makes the connection
to the thermal scale.

Now consider a binary mixture where each location of the system is on a lattice. Each "site" of the lattice may be occupied by one of two end-members/molecules, $N_1$ or $N_2$. The total number of available sites is $N=N_1+N_2$ (the lattice is full). As each end-member/molecule is indistinguishable from other molecules of the same type, the number of ways we can arrange these molecules are, \begin{align*} W = N! / (N_1!\, N_2!) \end{align*} The entropy is then, \begin{align*} S &= R \ln \left(\frac{N!}{N_1!\,N_2!}\right) \\ &= R \left(\ln \left(N!\right) - \ln\left(N_1!\right) - \ln\left(N_2!\right)\right) \end{align*} Using Stirling's approximation $\lim_{n\to\infty}\ln(n!)=n\,\ln(n)-n$. \begin{align*} S &= R \left(N\,\ln \left(N\right) - N_1\,\ln\left(N_1\right) - N_2\ln\left(N_2\right)\right) \\ &= R\,N \left(\ln \left(N\right) - \frac{N_1}{N}\,\ln\left(N_1\right) - \frac{N_2}{N}\ln\left(N_2\right)\right) \\ &= R\,N \left(\ln \left(N\right) - \frac{N_1}{N}\,\ln\left(N_1\right) - \frac{N_2}{N}\ln\left(N_2\right)\right) \\ &= -R\,N \left(\frac{N_1}{N}\,\ln\left(\frac{N_1}{N}\right) + \frac{N_2}{N}\ln\left(\frac{N_2}{N}\right)\right) \end{align*} For a gas, consider that the number of sites is actually the number of molecules in the system and all we are doing is allowing the gases to mix (each species may substitute for another species provided the overall molar balance is maintained), then it is clear that: \begin{align*} S = -R\,N \left(x_1\ln x_1 + x_2\ln x_2\right) \end{align*} In comparison with the ideal gas properties, this is the entropy of mixing as derived from Dalton's law!

### Ideal-gas isobaric contribution $\mu_{i}^{(ig)}(T)$

The term $\mu_{i}^{(ig)}(T)$ is a function of temperature which is directly related to the thermodynamic properties of a single isolated molecule (molecules of ideal gases do not interact). There are many theoretical approaches to expressing these terms for solids, simple models (i.e., rotors), and quantum mechanics can even directly calculate values for real molecules. However these results are obtained, this term is typically parameterised using polynomial equations.

The most common parameterisation is to use a polynomial for the heat capacity. This is common to the NIST, NASA CEA, ChemKin, and many other thermodynamic databases. \begin{align}\label{eq:cppoly} \bar{c}_{p,i}^{(ig)}(T) &= \sum_{k} C_k\,T^k \end{align} A heat capacity polynomial can be related to the enthalpy and entropy through two additional constants, To demonstrate how this is correct, consider a closed system, \begin{align*} {\rm d}H = T\,{\rm d} S + V\,{\rm d}p. \end{align*} At constant pressure the following expression holds, \begin{align*} \left(\frac{\partial S}{\partial T}\right)_{p} = \frac{1}{T}\left(\frac{\partial H}{\partial T}\right)_p=\frac{C_p}{T}. \end{align*} Thus, \begin{align*} H^{(ig)}(T) - H^{(ig)}(T_{ref}) &= \int_{T_{ref}}^T C_p^{ig}(T)\,{\rm d}T & S^{(ig)}(T)-S^{(ig)}(T_{ref}) &= \int_{T_{ref}}^T\frac{C_p^{(ig)}(T)}{T}{\rm d}T \end{align*} where $T_{ref}$ is a reference Performing this integration for the polynomial of Eq.\eqref{eq:cppoly}, \begin{align*} \bar{h}_i^{(ig)}(T) &= \int \bar{c}_{p,i}^{(ig)}(T) {\rm d}T = C_{-1}\,\ln T + \sum_{k\neq-1} \frac{C_k}{k+1}\,T^{k+1} + \bar{h}_i^0 \\ \bar{s}_i^{(ig)}(T) &= \int \frac{\bar{c}_{p,i}^{(ig)}(T)}{T} {\rm d}T = C_{0}\,\ln T + \sum_{k\neq 0} \frac{C_k}{k}\,T^{k} + \bar{s}_i^0, \end{align*} where the two additional constants $\bar{h}_i^0$ and $\bar{s}_i^0$ must be determined through construction of a thermodynamic path to a known value of the entropy and enthalpy. This is usually through heats of reaction, solution, and equilibria such as vapour pressure, linking to the elements at standard conditions.

To allow this data to be expressed as a single function in Simcem, it is expressed directly as $\mu_i^{(ig)}(T)$ using the following transformation, \begin{align}\label{eq:Mu_shomate_form} \mu^{(ig)}_{i}(T) &= \tilde{C}_{lnT}\ln T + \tilde{C}_{TlnT}\,T\,\ln T + \sum_{k} \tilde{C}_k\,T^{k} \end{align} where, \begin{align*} \tilde{C}_k = \begin{cases} C_{-1} & \text{for $k=\ln T$.}\\ -C_{0} & \text{for $k=T\ln T$.}\\ C_{-1}+h_i^0 & \text{for $k=0$.}\\ C_{0}-s_i^0 & \text{for $k=1$.}\\ -\frac{C_{k-1}}{(k-1)k} & \text{for all other polynomial $C_p$ terms.} \end{cases} \end{align*}

The most comprehensive (and most importantly, free!) source of ideal gas heat capacities and constants is the NASA CEA database, but additional constants are widely available. When collecting $c_p$ data, a distinction must be made between the calculated ideal gas properties and the measured data. The measured data may include additional contributions from the interactions of the molecules and thus must be fitted with the full chemical potential model and not just the ideal gas model.

\begin{align*} \mu^{(ig)}_{i}(T) &= \bar{h}_i^{(ig)}(T) - T\,\bar{s}_i^{(ig)}(T) \\ &= C_{-1}\,\ln T + \sum_{k\neq-1} \frac{C_k}{k+1}\,T^{k+1} + h_i^0 - C_{0}T\ln T - \sum_{k\neq 0} \frac{C_k}{k}\,T^{k+1} - T\,s_i^0 \\ &= \left(C_{-1}- C_{0}\,T\right)\ln T + h_i^0 - T\,s_i^0 + \sum_{k\neq-1} \frac{C_k}{k+1}\,T^{k+1} - \sum_{k\neq 0} \frac{C_k}{k}\,T^{k+1} \\ &= \left(C_{-1}- C_{0}\,T\right)\ln T + h_i^0 +C_{-1} + T\left(C_0-\,s_i^0\right) + \sum_{k\neq\{-1,0\}} \frac{C_k}{k+1}\,T^{k+1} - \sum_{k\neq \{-1,0\}} \frac{C_k}{k}\,T^{k+1} \\ &= \left(C_{-1}- C_{0}\,T\right)\ln T + h_i^0 +C_{-1} + T\left(C_0-\,s_i^0\right) + \sum_{k\neq\{-1,0\}} \left(\frac{1}{k+1}-\frac{1}{k}\right)C_k\,T^{k+1} \\ &= \left(C_{-1}- C_{0}\,T\right)\ln T + h_i^0 +C_{-1} + T\left(C_0-\,s_i^0\right) - \sum_{k\neq\{-1,0\}} \frac{C_k}{k(k+1)}T^{k+1} \end{align*} Comparison with Eq.\eqref{eq:Mu_shomate_form} yields the definitions of the coefficients.

### Incompressible phase

Another useful reference state is the incompressible phase. This model is used to describe liquid and solid phases when an equation of state describing all phases at once is unavailable. The generating chemical potential is as follows, \begin{align*} \mu^{(ip)}_{i}(T,\,p) &= \mu^{(ip)}_{i}(T) + \bar{v}^0_{i}\left(p-p_0\right) + R\,T\ln \left(\frac{N_{i}}{N_\alpha}\right) \end{align*} where $\mu^{(ip)}_{i}(T)$ is again a temperature-dependent term and $\bar{v}^0_{i}$ is the (constant) partial molar volume of the incompressible species $i$. Using the expression for the chemical potential as a generating function all required properties are recovered, \begin{align*} \bar{v}_{i} &= \bar{v}^0_{i} & \bar{s}_{i} &= -\frac{\partial \mu^{(ip)}_{i}(T)}{\partial T} - R \ln\left(\frac{N_{i}}{N_{\alpha}}\right) \\ \bar{c}_{p,i} &= -T\frac{\partial^2 \mu^{(ip)}_{i}(T)}{\partial T^2} & \bar{\alpha}_{i} &= 0 \\ \bar{\beta}_{T,i} &= 0 & & \end{align*} As $\alpha$ and $\beta_T$ are zero, some thermodynamic properties are not well specified but can be determined by other means. For example, $\bar{c}_p^{(ip)}=\bar{c}_v^{(ip)}$. In Simcem, ideal solids still contain a mixing entropy term which is included to allow ideal solid solutions to be constructed and used as a base class for more complex models.

## Helmholtz Models ($V$ as a variable)

In this section, the calculation of the required properties for minimisation with phases specified by the following set of variables is considered, \begin{align*} \vec{X}_\alpha &= \left\{T_\alpha,\,V_\alpha,\,\left\{N_{i}\right\}^{N_{C,\alpha}}\right\} . \end{align*} In this variable set, the required non-zero derivatives to compute the constraint functions are as follows, \begin{align}\label{eq:VMatDeriv1} \left(\frac{\partial p_\alpha}{\partial T}\right)_{V,\left\{N_{i}\right\}} = -\frac{\left(\alpha\,V\right)_{\left\{N_{i}\right\}}}{\left(\beta_T\,V\right)_{\left\{N_{i}\right\}}}, \\ \label{eq:VMatDeriv2} \left(\frac{\partial p_\alpha}{\partial V}\right)_{T,\left\{N_{i}\right\}} = -\left(\beta_T\,V\right)_{\left\{N_{i}\right\}}^{-1}, \\ \left(\frac{\partial p_\alpha}{\partial N_i}\right)_{T,V,\left\{N_{j\neq i}\right\}} \end{align} These derivatives are convenient to determine in this variable set and also indirectly specify $\alpha$ and $\beta_T$ which are two of the three required material derivatives to allow use of the Bridgman tables. The third molar derivative is also convenient to determine and provides a direct path to the molar volume as given in Eq.\eqref{eq:TVpartialmolarvolume} and reproduced below, \begin{align} \bar{v}_i = - \left(\frac{\partial p}{\partial N_i}\right)_{T,V,\left\{N_{j\neq i}\right\}} \left(\frac{\partial p}{\partial V}\right)_{T,\left\{N_{j}\right\}}^{-1}. \end{align} Thus, these three derivatives are required to be implemented by all models and used to derive other properties. The derivatives of the potentials are as follows \begin{align} \left(\frac{\partial f}{\partial T}\right)_{V,\left\{N_i\right\}},\, \left(\frac{\partial f}{\partial V}\right)_{T,\left\{N_i\right\}},\, \left(\frac{\partial f}{\partial N_i}\right)_{T,V,\left\{N_{i\neq j}\right\}}, \end{align} where $f=\left\{H,G,-S,U,A\right\}$. The derivatives for each potential are specified below in terms of the most convenient properties/derivatives for this variable set.

$Y_1$ | $Y_2$ | $f$ | $\left(\frac{\partial f}{\partial N_i}\right)_{T,V,\left\{N_{j\neq i}\right\}}$ | $\left(\frac{\partial f}{\partial T}\right)_{V,\left\{N_j\right\}}$ | $\left(\frac{\partial f}{\partial V}\right)_{T,\left\{N_j\right\}}$ |
---|---|---|---|---|---|

$p$ | $T$ | $G$ | $\breve{g}_i$ | $V\,\left(\frac{\partial p}{\partial T}\right)_{V,\left\{N_i \right\}} -S$ | $V \left(\frac{\partial p}{\partial V}\right)_{T,\left\{N_i\right\}}$ |

$V$ | $T$ | $A$ | $\breve{a}_i$ | $-S$ | $-p$ |

$p$ | $S$ | $H$ | $\breve{h}_i$ | $C_{V,\left\{N_i\right\}}+V\,\left(\frac{\partial p}{\partial T}\right)_{V, \left\{N_i\right\}}$ | $V\left(\frac{\partial p}{\partial V}\right)_{T,\left\{N_i\right\}}+T\left(\frac{\partial p}{\partial T}\right)_{V,\left\{N_i\right\}}$ |

$p$ $V$ | $H$ $U$ | $-S$ | $-\breve{s}_i$ | $-\frac{C_{V,\left\{N_i\right\}}}{T}$ | $-\left(\frac{\partial p}{\partial T}\right)_{V,\left\{N_i\right\}}$ |

$V$ | $S$ | $U$ | $\breve{u}_i$ | $C_{V,\left\{N_i\right\}}$ | $T\,\left(\frac{\partial p}{\partial T}\right)_{V,\left\{N_i\right\}} - p$ |

where $C_V = \left(\partial U/\partial T\right)_V$, and is related to the isobaric heat capacity using the following relationship, \begin{align}\label{eq:CpCv} C_{p,\left\{N_i\right\}} = C_{v,\left\{N_i\right\}} -T\,\left(\frac{\partial p}{\partial T}\right)_{V\,\left\{N_i\right\}}^2 \left(\frac{\partial p}{\partial V}\right)_{T,\left\{N_i\right\}}^{-1} \end{align} In summary, models should provide equations to calculate $p$, $\breve{a}_i$, $\breve{u}_i$, $C_{v,\left\{N_i\right\}}$, $\left(\partial p/\partial V\right)_{T,\left\{N_i\right\}}$, $\left(\partial p/\partial T\right)_{V,\left\{N_i\right\}}$, $\left(\partial p/\partial N_i\right)_{T,V,\left\{N_{j\neq i}\right\}}$. These derivatives are used to compute the frozen $\alpha$ and $\beta_T$ values using Eqs.\eqref{eq:VMatDeriv1} and \eqref{eq:VMatDeriv1}. The third material derivative, $C_{p,\left\{N_i\right\}}$, is then obtained from Eq.\eqref{eq:CpCv}. The partial molar volume is calculated from Eq.\eqref{eq:TVpartialmolarvolume}. The partial volar/molar quantities are related by Eq.\eqref{eq:molartovolarproperty}; however, the partial volar Helmholtz free energy is equal to the chemical potential, just like any molar derivative of a potential when its natural variables are held constant (see Eq.\eqref{eq:ChemPotDefinition}). Thus, the partial volar/molar Gibbs and Helmholtz free energies are closely related, \begin{align} \mu_i &= \breve{a}_i \\ \bar{a}_i &= \breve{a}_i - p\,\bar{v}_i \\ \breve{g}_i &= \mu_i - \bar{v}_i\,V\,\left(\frac{\partial p}{\partial V}\right)_{T,\left\{N_j\right\}} \end{align} All other required partial thermodynamic potential properties are derived using straightforward applications of Eq.\eqref{eq:molartovolarproperty} and the partial molar relations of Eqs.\eqref{eq:partialmolarrelationstart}-\eqref{eq:partialmolarrelationend}, the results of which are below, \begin{align} \bar{u}_i &= \breve{u}_i +\bar{v}_i\left(T\left(\frac{\partial p}{\partial T}\right)_{V,\left\{N_i\right\}} - p\right) \\ \bar{s}_i &= \frac{\bar{u}_i - \bar{a}_i}{T} \\ \breve{s}_i &= \bar{s}_i - \bar{v}_i \left(\frac{\partial p}{\partial T}\right)_{V\,\left\{N_j\right\}} \\ \bar{h}_i &= \bar{u}_i + p\,\bar{v}_i \\ \breve{h}_i &= \bar{h}_i - \bar{v}_i\left(V\left(\frac{\partial p}{\partial V}\right)_{T,\left\{N_i\right\}}+T\left(\frac{\partial p}{\partial T}\right)_{V,\left\{N_i\right\}}\right) \end{align}

### Ideal gas model

As experimental ideal-gas data is usually presented in terms of isobaric polynomials, the ideal gas functions here are based on the previous form of the ideal gas model given in Eq.\eqref{eq:idealgasmuTP}. Transforming the pressure variable to phase volume yields the following definition of the chemical potential and partial volar Helmholtz free energy, \begin{align} \breve{a}_i = \mu_i^{(ig)}(T,V_\alpha) &= \mu_{i}^{(ig)}(T) + R\,T\ln \left(\frac{N_\alpha\,R\,T}{V_\alpha\,p_0}\right) + R\,T\ln \left(\frac{N_{i}}{N_\alpha}\right). \end{align} The pressure derivatives are obtained from the well known relation $p=N\,R\,T/V$, \begin{align} \left(\partial p/\partial V\right)_{T,\left\{N_i\right\}} &= - N\,R\,T/V^2 \\ \left(\partial p/\partial T\right)_{V,\left\{N_i\right\}} &= N\,R/V \\ \left(\partial p/\partial N_i\right)_{T,V,\left\{N_{j\neq i}\right\}} &= R\,T/V \end{align} Finally, the internal energy and heat capacity are as follows, \begin{align} C_{v,\left\{N_i\right\}}^{(ig)}&=-T\sum_i N_i\,\frac{\partial^2 \mu_i(T)}{\partial T^2} - N\,R \\ \breve{u}_i^{(ig)} = \bar{u}_i^{(ig)} &= \mu_i^{(ig)}(T) - R\,T-T\frac{\partial \mu_i^{(ig)}(T)}{\partial T} \end{align}

## Other models

### Excess phases

Assume a phase/model is infinite in quantity, therefore its extensive state variables must be excluded from the system under study. The infinite phase can exchange mass with the system under study, thus it has a chemical potential. The change in the Gibb's free energy of the (infinite) system can be defined in terms of the changes in the species within it: \begin{align} \Delta G = \sum_i \mu_i\,\Delta N_i \end{align} where $\mu_i$ is the constant chemical potential. We note then that the entropy change and volume change of this system due to variation in the system temperature/pressure must be zero.

## Automatic differentiation

The best introduction to automatic differentiation I found are here from the 2010 Sintef winter school and these notes are largely based on those notes.

Our goal is to evaluate a function and its derivatives at once. By considering how derivatives are propagated through operators/functions we see that the $k$th derivative of a function $f(g)$ can be expressed in terms of the first $k$ derivatives of $g$. To prove this, consider the Taylor expansion of a general function $f$: \begin{align*} f(a+x) &= f(a) + \frac{1}{1!}\left(\frac{\partial f}{\partial x}\right)(a)(x-a) + \frac{1}{2!}\left(\frac{\partial^2 f}{\partial x^2}\right)(a)(x-a)^2 + \ldots \\ &= \sum_{k=0}^\infty \frac{1}{k!}\left(\frac{\partial^k f}{\partial x^k}\right)(a)(x-a)^k \\ &= \sum_{k=0}^\infty f_k(a)(x-a)^k \end{align*} Each Taylor term $f_k$ corresponds to the $k$th derivative (divided by $k!$). Thus if we can determine how the Taylor coefficients of the result of an operation are related to the Taylor coefficients of its arguments, we can propogate derivative information through the operation.

### Basic arithmetic

To demonstrate this we start from two basic definitions: one for the Taylor coefficients of a constant, $c$: \begin{align*} (c)_k = \begin{cases} c & \text{for $k=0$}\\ 0 & \text{for $k>0$} \end{cases}, \end{align*} and another for the Taylor coefficients of a variable, $x$: \begin{align*} (x)_k = \begin{cases} x & \text{for $k=0$}\\ 1 & \text{for $k=1$}\\ 0 & \text{for $k>1$} \end{cases}. \end{align*} Our first operations are addition and subtraction, where it is easy to see the Taylor series of the result is trivially calcluated from the Taylor series of the arguments: \begin{align*} \left(f+g\right)_k &= f_k + g_k\\ \left(f-g\right)_k &= f_k - g_k \end{align*} For multiplication, consider Taylor expansions of both the result and the two functions: \begin{align*} \sum_{k=0}^\infty \left(f\,g\right)_k (x-a)^k &= \sum_{l=0}^\infty f_l (x-a)^l \sum_{m=0}^\infty g_m (x-a)^m \\ &= \sum_{l=0}^\infty\sum_{m=0}^\infty f_l\,g_m (x-a)^{l+m} \end{align*} We can then match terms in equal powers of $(x-a)$ on either side of the equation to yield the final result: \begin{align*} \left(f\,g\right)_k = \sum_{i=0}^k f_{i} g_{k-i} \end{align*} Division is slightly more challenging, instead consider that $f = (f\div g) g$, and insert the taylor expansions: \begin{align*} \sum_{k=0}^\infty f_k (x-a)^k &= \sum_{l=0}^\infty(f\div g)_l (x-a)^l \sum_{m=0}^\infty g_m (x-a)^m \\ &= \sum_{l=0}^\infty \sum_{m=0}^\infty (f\div g)_l g_m (x-a)^{l+m} \end{align*} Again equating terms with equal powers exactly as with multiplication: \begin{align*} f_k &= \sum_{l=0}^k(f\div g)_l g_{k-l} \\ &= \sum_{l=0}^{k-1} (f\div g)_l g_{k-l} + (f\div g)_k g_0 \\ (f\div g)_k &= \frac{1}{g_0}\left(f_k - \sum_{l=0}^{k-1} (f\div g)_l g_{k-l}\right) \end{align*} Where we removed our target term from the sum in the second line, then rearranged to make it the subject in the third line.

### Special functions

Now we focus on special functions. To solve these, we need the general derivative of a taylor series: \begin{align*} \frac{\partial f}{\partial x} = \sum_{k=1}^\infty k\,f_k (x-a)^{k-1} \end{align*} Considering $\ln$ first, we have the following differential relationship: \begin{align*} \frac{\partial}{\partial x} \ln f &= \frac{\partial f}{\partial x} \frac{1}{f} \\ f \frac{\partial}{\partial x} \ln f &= \frac{\partial f}{\partial x}, \end{align*} where we have multiplied by $f$ to avoid polynomial division when the Taylor series are inserted. Doing this now, \begin{align*} \sum_{i=0}^\infty f_i (x-a)^i \sum_{j=1}^\infty j\left(\ln f\right)_j (x-a)^{j-1} &= \sum_{k=1}^\infty k\, f_k (x-a)^{k-1} \end{align*} Multiplying both sides by $(x-a)$, factoring common terms as with multiplication: \begin{align*} \sum_{i=0}^\infty f_i \sum_{j=1}^\infty j\left(\ln f\right)_j (x-a)^{i+j} &= \sum_{k=1}^\infty k\,f_k (x-a)^k \\ \sum_{i=1}^k f_{k-i}\,i\left(\ln f\right)_i &= k\,f_k & \text{for $k>0$} \\ f_0\,k\,\left(\ln f\right)_k + \sum_{i=1}^{k-1} f_{k-i}\,i\left(\ln f\right)_i &= k\,f_k & \text{for $k>0$} \\ \left(\ln f\right)_k &= \frac{1}{f_0}\left(f_k - \frac{1}{k}\sum_{i=1}^{k-1} i\,f_{k-i}\,\left(\ln f\right)_i\right) & \text{for $k>0$} \end{align*} Where this expression only applies for $k>0$ as the first line has no constant terms within it; however, we have the trivial identity $(\ln f)_0 = \ln f_0$. Sine and Cosine have to be calcluated at the same time (regardless of which one is required). Starting from the general differentiation rules \begin{align*} \frac{\partial}{\partial x} \sin f &= \frac{\partial f}{\partial x} \cos f \\ \frac{\partial}{\partial x} \cos f &= -\frac{\partial f}{\partial x} \sin f \end{align*} Inserting the Taylor expansions and grouping terms with identical coefficients the result is found: \begin{align*} \left(\sin f\right)_k &= \frac{1}{k}\sum_{i=1}^k i\,f_i \left(\cos f\right)_{k-i} & \text{for $k>0$} \\ \left(\cos f\right)_k &= -\frac{1}{k}\sum_{i=1}^k i\,f_i \left(\sin f\right)_{k-i} & \text{for $k>0$}, \end{align*} again we have $\left(\sin f\right)_0 = \sin f_0$ and in general $\left(f(g)\right)_0 = f(g_0)$.

### Generalised Power law

Finally, we calculate the generalized power rule: \begin{align*} \frac{\partial}{\partial x} f^g &= f^g \left(\frac{\partial f}{\partial x} \frac{g}{f} + \frac{\partial g}{\partial x} \ln f\right) \\ f \frac{\partial}{\partial x} f^g &= f^g \frac{\partial f}{\partial x} g + f^g \frac{\partial g}{\partial x} f\,\ln f \end{align*} Inserting Taylor expansions: \begin{align*} \sum_{i=0}^\infty \sum_{j=1}^\infty j\,f_i\left(f^g\right)_j \left(x-a\right)^{i+j-1} &= \sum_{i=0}^\infty \left(f^g\right)_i\sum_{j=1}^\infty j\,f_j \sum_{k=0}^\infty g_k \left(x-a\right)^{i+j-1+k} \\ &\qquad + \sum_{i=0}^\infty \left(f^g\right)_i\sum_{j=1}^\infty j\,g_j \sum_{k=0}^\infty f_k \sum_{l=0}^\infty \left(\ln f\right)_l\left(x-a\right)^{i+j-1+k+l} \end{align*} Multiplying both sides by $\left(x-a\right)$ and grouping common indices/terms: \begin{align*} \sum_{i=0}^\infty \sum_{j=1}^\infty j\,f_i\left(f^g\right)_j \left(x-a\right)^{i+j} &= \sum_{i=0}^\infty \sum_{j=1}^\infty \sum_{k=0}^\infty j\left(f^g\right)_i\left[f_j g_k \left(x-a\right)^{i+j+k} + g_j f_k \sum_{l=0}^\infty \left(\ln f\right)_l\left(x-a\right)^{i+j+k+l}\right] \end{align*} Selecting all terms with the $m$th power: \begin{align*} \sum_{j=1}^m j\,f_{m-j}\left(f^g\right)_j &= \sum_{j=1}^{m} \sum_{i=0}^{m-j} j\left(f^g\right)_i\left[f_j g_{m-i-j} + g_j \sum_{k=0}^{m-i-j} f_k \left(\ln f\right)_{m-i-j-k}\right] \end{align*} Factoring the $m$th term from the left-hand side: \begin{align*} \left(f^g\right)_m &= \frac{1}{m\,f_0}\left(\sum_{j=1}^{m} \sum_{i=0}^{m-j} j\left(f^g\right)_i\left[f_j g_{m-i-j} + g_j \sum_{k=0}^{m-i-j} f_k \left(\ln f\right)_{m-i-j-k}\right] - \sum_{j=1}^{m-1} j\,f_{m-j}\left(f^g\right)_j\right) \end{align*} If the power is a constant, i.e., $g=a$, then the following simplified expression is obtained: \begin{align*} \left(f^g\right)_m &= \frac{1}{m\,f_0}\left(\sum_{j=1}^{m} j\left(f^g\right)_{m-j} f_j\,a - \sum_{j=1}^{m-1} j\,f_{m-j}\left(f^g\right)_j\right) \\ &= \frac{1}{m\,f_0}\sum_{j=1}^{m}\left(f^g\right)_{m-j} f_j\left(j\,a - m+j\right) \\ &= \frac{1}{f_0}\sum_{j=1}^{m}\left(\frac{(a + 1)j}{m} - 1\right) f_j \left(f^g\right)_{m-j}. \end{align*}