\section{Evaluating partial molar properties required by the Kempers model}
\label{sec:eos_differentials}

The Kempers model requires several partial molar properties, most of which are directly available in the EoS-package ThermoPack.\cite{thermopack} Here, the quantity $\ppder{\mu_i}{x_j}_{T,p}$ is related to other quantities that are directly available.

Begin from the definition of the chemical potential as the derivative of Helmholtz energy and invert the order of differentiation,

\begin{equation}
    \ppder{\mu_i}{n_j}_{T,p} = \left(\pder{}{n_j}\ppder{A}{n_i}_{T,V} \right)_{T,p} = \left(\pder{}{n_i}\ppder{A}{n_j}_{T,p}\right)_{T,V}.
\label{eq:dmudn_TV_1}
\end{equation}

From the total differential of the Helmholtz energy, $dA = -SdT - pdV + \sum_i \mu_i dn_i$,

\begin{equation}
    \begin{split}
        \ppder{A}{n_j}_{T,p} &= -p \ppder{V}{n_j}_{T,p} + \mu_j\\
        \left(\pder{}{n_i}\ppder{A}{n_j}_{T,p}\right)_{T,V} &= \ppder{\mu_j}{n_i}_{T,V} - v_j \ppder{p}{n_i}_{T,V} - p \underbrace{\ppder{v_j}{n_i}_{T,V}}_{=0}\\
        \ppder{\mu_i}{n_j}_{T,p} &= \ppder{\mu_j}{n_i}_{T,V} - v_j \ppder{p}{n_i}_{T,V}.
    \end{split}
\end{equation}

The total differential of the chemical potential at constant pressure and temperature is 
\begin{equation}
    d\mu_i = \sum_j \ppder{\mu_i}{n_j}_{T,p,n_{k \neq j}} dn_j.
\end{equation}

Further inserting for the total differential of the mole number of each species yields
\begin{equation}
    \begin{split}
        d\mu_i &= \sum_j \ppder{\mu_i}{n_j}_{T,p,{n_{k \neq j}}}\left(x_j dn + n dx_j\right)\\
        d\mu_i &= \sum_j \ppder{\mu_i}{n_j}_{T,p,{n_{k \neq j}}} n dx_j + \left( \underbrace{\sum_j \ppder{\mu_i}{n_j}_{T,p,{n_{k \neq j}}} x_j}_{=0}\right) dn \\
        &= \sum_j \ppder{\mu_i}{x_j}_{T, p, n_{k \neq j}}dx_j
    \end{split}
\end{equation}

where the second term on the right hand side of the second line vanishes by the Gibbs-Duhem relation, and the partial derivative of interest may be readily identified as $\ppder{\mu_i}{x_j}_{T, p, n_{k \neq j}} = n\ppder{\mu_i}{n_j}_{T, p, n_{k \neq j}} $. This expression was validated numerically by evaluating the total derivative of $\mu_i$ at varying composition and total number of moles.



\begin{comment}
    This section will summarise how the various differentials presented in Table \ref{tab:eos_diffs} may be related to other differentials implemented in ThermoPack. Lowercase letters will consistently denote molar, or partial molar quantities, such that for the volume, $V$, the molar volume is denoted $v$, and the partial molar volume is denoted $v_i$.
    
    \subsubsection{\mathsec{$\ppder{\mu_1}{x_1}_{\rho, T}$}}
    The Helmholtz energy differential in volume and mole numbers is related to the differential in density by the total differential of the density
    \begin{equation}
        \begin{split}
            \drho &= - \frac{n}{V^2}\dV + \frac{1}{V} \dn\\
            \dV &= - Vv\drho + v\dn
        \end{split}
        \label{eq:dVdrhodn}
    \end{equation}
    which gives
    \begin{equation}
        \begin{split}
            \dA &= - S \dT - p \dV + \sum_i \mu_i \dn_i\\
            &= -S\dT - p \left(-Vv\drho + v\dn\right) + \sum_i \mu_i \dn_i\\
            &= - S \dT + pVv\drho + \sum_i \left(\mu_i - pv\right)\dn_i\\
            &\implies \ppder{A}{n_i}_{T,\rho} = \mu_i - pv
        \end{split}
    \end{equation}
    The differential $\ppder{\mu_i}{n_j}_{T,\rho,n_k}$ can then be identified by reversing the order of differentiation
    \begin{equation}
        \begin{split}
            \ppder{\mu_i}{n_j}_{T,\rho,n_k} &= \left(\pder{}{n_j}\ppder{A}{n_i}_{T,V,n_k}\right)_{T, \rho, n_k}\\
            &= \left(\pder{}{n_j}\ppder{A}{n_i}_{T, \rho, n_k}\right)_{T,V,n_k}\\
            &= \left(\pder{}{n_j}(\mu_i - pv)\right)_{T,V,n_k}\\
            &= \ppder{\mu_i}{n_j} - p \underbrace{\ppder{v}{n_j}_{T,V,n_k}}_{=0} - v \ppder{p}{n_j}_{T,V,n_k}\\
            &= \ppder{\mu_j}{n_i}_{T,V} - v\ppder{p}{n_j}_{T,V}.
        \end{split}
        \label{eq:dmudn_Trho_1}
    \end{equation}
    
    From this, it is clear that at constant temperature and density,
    \begin{equation}
        \begin{split}
            \dmu_i &= \sum_j \ppder{\mu_i}{n_j}_{T,\rho,n_j}\dn_j\\
                &= \sum_j \ppder{\mu_i}{n_j}_{T,\rho}(n\dx_j + x_j\dn)\\
                &\implies \ppder{\mu_i}{x_j}_{T,\rho} = n \ppder{\mu_i}{n_j}.
        \end{split}
    \end{equation}
    This expression was tested numerically for a binary system by varying the mole number of both components at constant $T, \rho$ and examining the total differential of the chemical potential. That is, verifying that for each step,
    \begin{equation}
        \Delta \mu_i = \ppder{\mu_i}{x_1}_{T,\rho}\Delta x_1 + \ppder{\mu_i}{x_2}_{T,\rho}\Delta x_2.
    \end{equation}
    
    The intermediate relations for $\ppder{A}{n_i}_{T,\rho, n_k}$ and $\ppder{v}{n_j}_{T,V,n_k}$ were validated in the same manner.
    
    \subsubsection{\mathsec{$\ppder{\mu_i}{\rho}_{T,\vx}$}}
    The Maxwell-relations 
    \begin{align}
        \ppder{\mu_i}{V}_{T,\vn} &= - \ppder{p}{n_i}_{T,V} & \ppder{\mu_i}{T}_{V,\vn} &= - \ppder{S}{n_i}_{T,V}
    \end{align}
    are obtained from the Helmholtz energy, the total differential of the chemical potential can then be written as
    \begin{equation}
        \dmu_i = - \ppder{S}{n_i}_{T,V}\dT - \ppder{p}{n_i}_{T,V}\dV + \sum_j \ppder{\mu_i}{n_j}_{T,V} \dn_j
    \end{equation}
    Employing the relation between the volume, mole number and molar density differential in Eq. \eqref{eq:dVdrhodn} and the expression in \eqref{eq:dmudn_Trho_1} this can be rewritten as (note the swap of the $i,j$ indices in the derivative of the chemical potential with respect to mole numbers)
    \begin{equation}
        \begin{split}
            \dmu_i &= - \ppder{S}{n_i}_{T,V}\dT + \ppder{p}{n_i}_{T,V}Vv\drho - \ppder{p}{n_i}_{T,V}v\dn + \sum_j \left[\ppder{\mu_j}{n_i}_{T,\rho} + v\ppder{p}{n_i}_{T,V}\right]\dn_j\\
            &= - \ppder{S}{n_i}_{T,V}\dT + \ppder{p}{n_i}_{T,V}Vv\drho + \sum_j \ppder{\mu_j}{n_i}_{T,\rho} (x_j\dn + n \dx_j)\\
            &= - \ppder{S}{n_i}_{T,V}\dT + \ppder{p}{n_i}_{T,V}Vv\drho + \sum_j \ppder{\mu_j}{n_i}_{T,\rho} n \dx_j + \underbrace{\sum_j\ppder{\mu_j}{n_i}_{T,\rho}x_j}_{=0}\dn
        \end{split}
    \end{equation}
    The final term must be zero, as the chemical potential is an intensive property. The final term can be confirmed vanish by the Gibbs-Duhem relation at constant $T,p$ by reinserting the expression in Eq. \eqref{eq:dmudn_Trho_1}, further inserting the relation
    \begin{equation}
        \ppder{\mu_i}{n_j}_{T,V} - v_i \ppder{p}{n_j} = \ppder{\mu_j}{n_i}_{T,p}
    \end{equation}
    and applying the Gibbs-Duhem relation at constant temperature and pressure. One must also employ the symmetry of the second derivatives, specifically $\ppder{\mu_i}{n_j} = \ppder{\mu_j}{n_i}$.
    
    Now the derivative of interest can be identified as
    \begin{equation}
        \ppder{\mu_i}{\rho}_{T,\vx} = \ppder{p}{n_i}_{T,V}Vv.
    \end{equation}
    
    \subsubsection{Second derivatives of \mathsec{$\mu_i(T,\rho,\vx)$}}
    With the total differential of the chemical potential identified in the previous section, the second derivatives of interest can be efficiently identified
    \begin{equation}
        \spder{\mu_i}{x_j}{T}_{\rho} = - \left(\pder{}{x_j}\ppder{S}{n_i}_{T,V}\right)_{T,\rho}
    \end{equation}
    
    The third derivatives of the Helmholtz energy are not available through ThermoPack, these were therefore obtained by numerical differentiation of the second derivative properties. Numerical derivatives were computed as 
    \begin{equation}
        \begin{split}
            f(x + \epsilon) &= f(x) + \epsilon f'(x) + \frac{\epsilon^2}{2}f''(x) + \mathcal{O}(\epsilon^3)\\
            f(x - \epsilon) &= f(x) - \epsilon f'(x) + \frac{\epsilon^2}{2}f''(x) + \mathcal{O}(\epsilon^3)\\
            \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon} &= f'(x) + \mathcal{O}(\epsilon^3)
        \end{split}
    \end{equation}
    
    \subsubsection{\mathsec{$\ppder{\rho}{x_1}_{T,p}$ and $\ppder{\rho}{T}_{p,\vx}$}}
    The total differential of the density $\rho(n,V)$ is
    \begin{equation}
            \drho = \frac{1}{V}\dn -\frac{n}{V^2}\dV,
    \end{equation}
    where the volume differential may be expressed as 
    \begin{equation}
        \dV = \ppder{V}{T}_{p,\vn}\dT + \ppder{V}{p}_{T,\vn} \d p + \sum_i \ppder{V}{n_i}_{T,p,n_k}\dn_i
    \end{equation}
    and $n = \sum_i n_i$. Thereby,
    \begin{equation}
        \begin{split}
            \drho &= \frac{1}{V}\sum_i \dn_i - \frac{n}{V^2}\left(\ppder{V}{T}_{p,\vn}\dT + \ppder{V}{p}_{T,\vn} \d p + \sum_i \ppder{V}{n_i}_{T,p,n_k}\dn_i\right)\\
            &= - \frac{n}{V^2}\ppder{V}{T}_{p,\vn}\dT - \frac{n}{V^2}\ppder{V}{p}_{T,\vn} \d p + \sum_i \left(\frac{1}{V} - \frac{n}{V^2}v_i \right)\dn_i\\
            &= - \frac{n}{V^2}\ppder{V}{T}_{p,\vn}\dT - \frac{n}{V^2}\ppder{V}{p}_{T,\vn} \d p + \sum_i \left(\rho - \rho^2 v_i\right)\dx_i + \dn\underbrace{\sum_i\left(\frac{1}{V} - \frac{n}{V^2}v_i\right)x_i}_{=0}
        \end{split}
    \end{equation}
    where it is used that $\dn_i = n\dx_i + x_i \dn$, and the final sum vanishes as
    \begin{equation}
        \begin{split}
            \sum_i\left(\frac{1}{V} - \frac{n}{V^2}v_i\right)x_i &= \frac{1}{V}\left(\sum_i x_i - \frac{n}{V} \sum_i x_i v_i\right)\\
            &= \frac{1}{V} \left(1 - \frac{1}{v} v\right) = 0\\
        \end{split}
    \end{equation}
    Such that 
    \begin{equation*}
        \drho = - \frac{n}{V^2}\ppder{V}{T}_{p,\vn}\dT - \frac{n}{V^2}\ppder{V}{p}_{T,\vn} \d p + \sum_i \left(\rho - \rho^2 v_i\right)\dx_i
    \end{equation*}
    \begin{align}
        \ppder{\rho}{T}_{p,\vx} &= - \frac{n}{V^2}\ppder{V}{T}_{p,\vn} & \ppder{\rho}{p}_{T,\vx} &= - \frac{n}{V^2}\ppder{V}{p}_{T,\vn} & \ppder{\rho}{x_i}_{T,p} &= \rho - \rho^2 \ppder{V}{n_i}_{T,p}
    \end{align}
    
    Alternatively, the substitution $V = vn$ can be made to express the differential only in terms of intensive variables,
    \begin{align}
        \ppder{\rho}{T}_{p,\vx} &= - \rho^2 \ppder{v}{T}_{p,\vn} & \ppder{\rho}{p}_{T,\vx} &= - \rho^2\ppder{v}{p}_{T,\vn} & \ppder{\rho}{x_i}_{T,p} &= \rho - \rho^2 v_i
    \end{align}
    This can quickly be confirmed to hold for the case of an ideal gas, and has been tested numerically for non-ideal fluids.
\end{comment}