\subsection{The Kempers approach}
This derivation is inspired by the approach proposed by Kempers in 1989 and 2001. The underlying assumptions are that the steady state of a two bulb system connected by a tube of negligible volume, where the bulbs are at equal pressure and uniform temperature is the state with a maximum in the canonical partition function ($Q$). The bulbs will be denoted with superscript $\alpha$ and $\beta$, symbols lacking a superscript will be used to refer to the system as a whole. It is further assumed that the Soret coefficient of a mixture of hard spheres (denoted with superscript $HS$) is accurately described by kinetic gas theory. To explain the thought process behind the new suggested approach, this section begins by describing the central results presented by Kempers.

\begin{table}[tb!]
    \centering
    \begin{tabular}{c l c}
    \toprule
    Symbol & Description & Unit \\
    \midrule
        \symbline{\varphi}{Arbitrary variable}{-}
        \symbline{f}{Arbitrary function}{-}
        \symbline{T}{Temperature}{\kelvin}
        \symbline{v_i}{Partial molar volume of $i$}{\cubic\meter\per\mol}
        \symbline{V}{Total volume of system}{\cubic\meter}
        \symbline{\mu_i}{Chemical potential of $i$}{\joule\per\mol}
        \symbline{A}{Helmholtz energy}{\joule}
        \symbline{k_B}{Boltzmann constant}{\joule\per\kelvin}
        \symbline{n_i}{Moles of $i$}{\mol}
        \symbline{\vn}{Vector of mole numbers}{\mol}
        \symbline{\lambda_i, \phi, \gamma}{Lagrange multipliers}{-}
        \symbline{Q}{Canonical partition function}{-}
    \midrule
        \multicolumn{2}{c}{Superscripts} & \\
    \midrule
        \symbline{\alpha, \beta}{Bulb indices}{-}
        \symbline{ig}{Ideal gas state}{-}
        \symbline{HS}{Hard sphere state}{-}
    \midrule
        \multicolumn{2}{c}{Operators} & \\
    \midrule
        \symbline{\nabla}{Spatial gradient}{\per\meter}
        \symbline{\DeltaAB}{Difference between bulbs}{-}
        \symbline{\DeltaIG}{Difference between real and ideal gas state}{-}
        \symbline{\DeltaHS}{Difference between real and hard sphere state}{-}
    \bottomrule
    \end{tabular}
    \caption{Symbols used in the text}
    \label{tab:my_label}
\end{table}

\subsection{The 1989-approach}
The two bulbs may be treated as independent subsystems, under the constraint of particle conservation. In that case, $Q = Q^{\alpha}Q^{\beta}$ and $n_i^{\alpha} + n_i^{\beta} = n_i$ $\forall$ $i$, where $n_i$ with no subscript refers to the total number of moles in the system. Additionally, the two subsystems are held at constant and equal volume, such that $\sum_i v_i^{\alpha}n_i^{\alpha} = \sum_i v_i^{\beta}n_i^{\beta} = \frac{V}{2}$. Under the current assumptions, the steady state compositions for given temperatures $T^{\alpha}$ and $T^{\beta}$ is then given by the constrained maximisation
\begin{equation}
    \begin{split}
        \max_{n_i^\alpha, n_i^{\beta}}\left\{Q^{\alpha}Q^{\beta}\right\} &= \max_{n_i^\alpha, n_i^{\beta}}\left\{\exp\lrp-\frac{A^{\alpha}}{k_B T^{\alpha}} - \frac{A^{\beta}}{k_B T^{\beta}}\rrp\right\} \\ &=\min_{n_i^\alpha, n_i^{\beta}}\left\{\frac{A^{\alpha}}{T^{\alpha}} + \frac{A^{\beta}}{T^{\beta}}\right\} \\
        n_i^{\alpha} + n_i^{\beta} &= n_i, \hspace{.5cm} \forall \hspace{.5cm} i \\
        \sum_i v_i^{\alpha}n_i^{\alpha} &= \sum_i v_i^{\beta}n_i^{\beta}
    \end{split}
    \label{eq:1989_problem}
\end{equation}

The subscripts on the $\min$ and $\max$ operators indicate the optimisation variables, and are to be understood such that for some function $f(\varphi_1, \varphi_2)$,
\begin{equation}
    \begin{split}
        \max_{\varphi_1, \varphi_2}\left\{f\right\} = \lrp\varphi_1^{extrema}, \varphi_2^{extrema}\rrp = \min_{\varphi_1, \varphi_2}\left\{-f\right\}\\
        f\lrp\varphi_1^{extrema}, \varphi_2^{extrema}\rrp = f_{max}
    \end{split}
\end{equation}

By the Lagrange multiplier method, the target function corresponding to this optimisation is 
\begin{equation}
    \mathcal{L}\lrp \vn^\alpha, \vn^\beta, \mathbf{\lambda}, \gamma \rrp = \frac{A^{\alpha}}{T^{\alpha}} + \frac{A^{\beta}}{T^{\beta}} - \sum_i \lambda_i \lrp n_i^\alpha + n_i^{\beta} - n_i \rrp - \gamma \sum_i v_i^{\alpha}n_i^{\alpha} - v_i^{\beta}n_i^{\beta}
\end{equation}
differentiating with respect to the optimisation variables yields the set of equations
\begin{equation}
    \begin{split}
        \frac{\mu_i^{\alpha}}{T^{\alpha}} - \lambda_i - \gamma v_i^\alpha &= 0\\
        \frac{\mu_i^{\beta}}{T^{\beta}} - \lambda_i + \gamma v_i^\beta &= 0\\
    \end{split}
\end{equation}
Subtracting these two equations eliminates the $\lambda_i$, such that
\begin{equation}
    \begin{split}
        \frac{\mu_i^{\alpha}}{T^{\alpha}} - \frac{\mu_i^{\beta}}{T^{\beta}} - \gamma \lrp v_i^\alpha + v_i^\beta \rrp &= 0 \\
        \DeltaAB \pfrac{\mu_i}{T} - \gamma \lrp v_i^\alpha + v_i^\beta \rrp &= 0, \hspace{1cm} i = \{1, ..., N\}
    \end{split}
\end{equation}
where the notation $\DeltaAB\varphi \equiv \varphi^{\alpha} - \varphi^{\beta}$ has been introduced. Solving this equation for $\gamma$ with $i = k$ gives
\begin{equation}
    \gamma = \frac{1}{v_k^{\alpha} + v_k^{\beta}}\DeltaAB \pfrac{\mu_k}{T}
\end{equation}
Inserting this back into the remaining $N - 1$ equations gives
\begin{equation}
    \DeltaAB \pfrac{\mu_i}{T} - \pfrac{v_i^\alpha + v_i^\beta}{v_k^\alpha + v_k^\beta}\DeltaAB \pfrac{\mu_k}{T} = 0, \hspace{1cm} i = \{1, ..., N\} \setminus \{k\}.
\end{equation}
Dividing through by the length of the tube, $\Delta z$ and taking the limit $\Delta z \to 0$, $\DeltaAB T \to 0$ then gives the set of equations
\begin{equation}
    \nabla \pfrac{\mu_i}{T} - \pfrac{v_i}{v_k}\nabla \pfrac{\mu_k}{T} = 0, \hspace{1cm} i = \{1, ..., N\} \setminus \{k\}.
    \label{eq:kempers_89_result}
\end{equation}
So far, this derivation has exactly followed the approach proposed by Kempers in 1989. Keep this result in mind as we side-step to the 2001-approach.

\subsection{The 2001-approach}
In the 2001-approach, the function subject to maximisation is instead the ratio of the partition function to the partition function in the ideal gas state. That is, the objective is to find
\begin{equation}
    \max_{n_i^{\alpha}, n_i^{\beta}, n_i^{\alpha,ig}, n_i^{\beta,ig}} \lcp \frac{Q}{Q^{ig}}\rcp = \max_{n_i^{\alpha}, n_i^{\beta}, n_i^{\alpha,ig}, n_i^{\beta,ig}} \lcp \frac{Q^{\alpha}Q^{\beta}}{Q^{\alpha,ig} Q^{\beta,ig}}\rcp
\end{equation}
The constraints are expanded to include $n_i^{\alpha,ig} + n_i^{\beta,ig} = n_i$, where the superscript $ig$ indicates the ideal gas state. That is, the number of particles in the ideal gas state is treated as and independent variable, with the constraint that the total number of particles in the ideal gas state must equal that in the real state. Kempers does not place the constraint that $\sum_i v_i^{\alpha,ig}n_i^{\alpha,ig} = \sum_i v_i^{\beta,ig}n_i^{\beta,ig} = \frac{V}{2}$. The resulting optimisation problem is
\begin{equation}
    \begin{split}
        \min_{n_i^{\alpha}, n_i^{\beta}, n_i^{\alpha,ig}, n_i^{\beta,ig}} &\lcp \frac{A^{\alpha} - A^{\alpha,ig}}{T^{\alpha}} + \frac{A^{\beta} - A^{\beta,ig}}{T^{\beta}}\rcp \\
        n_i^{\alpha} + n_i^{\beta} &= n_i\\
        n_i^{\alpha,ig} + n_i^{\beta,ig} &= n_i\\
        \sum_i v_i^{\alpha}n_i^{\alpha} &= \sum_i v_i^{\beta}n_i^{\beta}
    \end{split}
\end{equation}
Because the mole numbers in the ideal gas state and the real state are treated as independent variables, and the Helmholtz energy in each state is only a function of the mole numbers in that state (i.e. $A^{\alpha} = A^{\alpha}(\vn^{\alpha}, T^{\alpha}, V^{\alpha})$, $A^{\alpha, ig} = A^{\alpha,ig}(\vn^{\alpha,ig}, T^{\alpha}, V^{\alpha})$ etc.), this is exactly equivalent to conducting the two independent optimisations
\begin{align}
        \min_{n_i^\alpha, n_i^{\beta}}&\left\{\frac{A^{\alpha}}{T^{\alpha}} + \frac{A^{\beta}}{T^{\beta}}\right\} & \min_{n_i^{\alpha,ig}, n_i^{\beta, ig}}&\left\{- \frac{A^{\alpha, ig}}{T^{\alpha}} - \frac{A^{\beta, ig}}{T^{\beta}}\right\} \\
        n_i^{\alpha} + n_i^{\beta} &= n_i & n_i^{\alpha, ig} + n_i^{\beta, ig} &= n_i \\
        \sum_i v_i^{\alpha}n_i^{\alpha} &= \sum_i v_i^{\beta}n_i^{\beta} & \\
\end{align}
The leftmost minimisation is exactly equivalent to the problem in Equation \eqref{eq:1989_problem}, resulting in the set of Equations \eqref{eq:kempers_89_result}. The rightmost minimisation yields the Lagrange target-function
\begin{equation}
    \mathcal{L}\lrp \vn^\alpha, \vn^\beta, \mathbf{\phi}\rrp = - \frac{A^{\alpha, ig}}{T^{\alpha}} - \frac{A^{\beta, ig}}{T^{\beta}} - \sum_i \phi_i (n_i^{\alpha, ig} + n_i^{\beta, ig} - n_i).
\end{equation}
Differentiation with respect to the target variables yields the set of equations
\begin{equation}
    \begin{split}
        - \frac{\mu_i^{\alpha,ig}}{T^{\alpha}} - \phi_i &= 0\\
        - \frac{\mu_i^{\beta,ig}}{T^{\beta}} - \phi_i &= 0
    \end{split}
\end{equation}
Subtracting the first equation from the second yields the set of equations
\begin{equation}
    \DeltaAB\pfrac{\mu_i^{ig}}{T} = 0.
\end{equation}
As before, dividing by the length of the tube $\Delta z$, and taking the limit $\Delta z \to 0$, $\DeltaAB T \to 0$ then gives
\begin{equation}
    \nabla \pfrac{\mu_i^{ig}}{T} = 0.
    \label{eq:kempers_2001_ig_result}
\end{equation}
To arrive at the set of equations proposed by Kempers in 2001, first subtract Equation \eqref{eq:kempers_2001_ig_result} for $i \neq k$ from equation \eqref{eq:kempers_89_result}, arriving at
\begin{equation}
    \begin{split}
        \nabla \pfrac{\mu_i}{T} - \nabla \pfrac{\mu_i^{ig}}{T} - \pfrac{v_i}{v_k}\nabla \pfrac{\mu_k}{T} &= 0, \hspace{1cm} i = \{1, ..., N\} \setminus \{k\}\\
        \DeltaIG \lcp \nabla \pfrac{\mu_i}{T} \rcp - \pfrac{v_i}{v_k}\nabla \pfrac{\mu_k}{T} &= 0
    \end{split}
    \label{eq:kempers_01_intermediate_1}
\end{equation}
where $\DeltaIG \varphi \equiv \varphi - \varphi^{ig}$. Finally, multiply Equation \eqref{eq:kempers_2001_ig_result} for $i = k$ with $\frac{v_i}{v_k}$ and add the result to Equation \eqref{eq:kempers_01_intermediate_1}, getting
\begin{equation}
    \DeltaIG \lcp \nabla \pfrac{\mu_i}{T} \rcp - \pfrac{v_i}{v_k} \DeltaIG \lcp \nabla \pfrac{\mu_k}{T} \rcp = 0.
    \label{eq:kempers_01_intermediate_2}
\end{equation}
To solve this equation for the Soret coefficient, begin by expanding the gradients
\begin{equation}
    \begin{split}
        \nabla \pfrac{\mu_i}{T} &= \lrp \pder{}{T} \frac{\mu_i}{T} \rrp_{p,x} \nabla T + \sum_j \lrp \pder{}{x_j} \frac{\mu_i}{T} \rrp_{T,p} \nabla x_j \\
        &= -\frac{h_i}{T^2} \nabla T + \frac{1}{T} \sum_j \ppder{\mu_i}{x_j}_{T,p} \nabla x_j
    \end{split}
    \label{eq:kempers_01_intermediate_3}
\end{equation}
Where the relationship
\begin{equation}
    \begin{split}
        \d H &= T \d S + V \d p + \sum_j \mu_j \dn_j\\
        \ppder{H}{n_i}_{T,p} \equiv h_i &= T \ppder{S}{n_i}_{T,p} + \mu_i\\
        &= - T \ppder{\mu_i}{T}_{\vn,p} + \mu_i\\
        &= - T^2 \lrp \pder{\pfrac{\mu_i}{T}}{T} \rrp_{\vn, p}
    \end{split}
\end{equation}
has been used. Insert this expression int Eq. \eqref{eq:kempers_01_intermediate_2} and divide through to obtain a set of $N - 1$ equations for the Soret coefficients,
\begin{equation}
    \DeltaIG \lcp -\frac{h_i}{T^2} - \frac{1}{T} \sum_j \ppder{\mu_i}{x_j} x_j (1 - x_j) S_{T,j} \rcp - \frac{v_i}{v_k} \DeltaIG \lcp -\frac{h_k}{T^2} - \frac{1}{T} \sum_j \ppder{\mu_k}{x_j} x_j (1 - x_j) S_{T,j} \rcp = 0
\end{equation}

\subsection{The new approach}

When one views the 2001-approach as two independent optimisations, it is not immediately clear why the real state is assumed to obey Equation \eqref{eq:kempers_89_result}, while the ideal gas state does not. The difference in the equations arises from the lack of a constraint on the bulb volumes in the ideal gas state.

If instead, one assumes that all systems obey Equation \eqref{eq:kempers_89_result}, and that the Soret coefficient in a mixture of hard spheres is accurately described by kinetic gas theory, one can set up Equation \eqref{eq:kempers_89_result} for both the hard sphere state and the real state,
\begin{equation}
    \begin{split}
        \nabla \pfrac{\mu_i}{T} - \frac{v_i}{v_k}\nabla \pfrac{\mu_k}{T} &= 0\\
        \nabla \pfrac{\mu^{HS}_i}{T} - \frac{v^{HS}_i}{v^{HS}_k}\nabla \pfrac{\mu^{HS}_k}{T} &= 0
    \end{split}
\end{equation}
Now, subtracting the two, as in the 2001-approach, and inserting the expansion of $\nabla \pfrac{\mu_i}{T}$ from Eq. \eqref{eq:kempers_01_intermediate_3}, one arrives at
\begin{equation}
    \DeltaHS \lcp -\frac{h_i}{T^2} - \frac{1}{T} \sum_j \ppder{\mu_i}{x_j} x_j (1 - x_j) S_{T,j} \rcp - \DeltaHS \lcp \frac{v_i}{v_k} \lrp -\frac{h_k}{T^2} - \frac{1}{T} \sum_j \ppder{\mu_k}{x_j} x_j (1 - x_j) S_{T,j} \rrp \rcp = 0,
\end{equation}
where $\DeltaHS \varphi \equiv \varphi - \varphi^{HS}$. Together with the condition that $\sum_j x_j(1 - x_j) S_{T,j} = 0$, and with $S_{T,i}^{HS}$ supplied from kinetic gas theory, this set of equations can be solved for the Soret coefficients in the real state.