# Consolidated theory of fluid thermodiffusion

Alexander Van-Brunt  Patrick E. Farrell  Charles W. Monroe
<EMAIL_ADDRESS>Mathematical Institute, University of Oxford,
Woodstock Road, Oxford, OX2 6GG, UK  Department of Engineering Science,
University of Oxford, Oxford, OX1 3PJ, United Kingdom  The Faraday
Institution, Harwell Campus, Didcot, OX11 0RA, United Kingdom

###### Abstract

We present the Onsager–Stefan–Maxwell thermodiffusion equations, which account
for the Soret and Dufour effects in multicomponent fluids by treating heat as
a pseudo-component. Unlike transport laws derived from kinetic theory, this
framework preserves the structure of the isothermal Stefan–Maxwell equations,
separating the thermodynamic forces that drive diffusion from the force that
drives heat flow. The Onsager–Stefan–Maxwell transport-coefficient matrix is
symmetric, and the second law of thermodynamics imbues it with simple spectral
characteristics. This new approach proves equivalent to both the intuitive
extension of Fick’s law and the generalized Stefan–Maxwell equations
popularized by Bird, Stewart, and Lightfoot. A general inversion process
allows the unique formulation of flux-explicit transport equations relative to
any choice of convective reference velocity. Stefan–Maxwell diffusivities and
thermal diffusion factors are tabulated for gaseous mixtures containing
helium, argon, neon, krypton, and xenon. The framework is deployed to perform
numerical simulations of steady three-dimensional thermodiffusion in a ternary
gas.

##  1 Introduction

Thermodiffusion generally describes situations where heat transfer is coupled
to mass transfer, a class of phenomena first recognized by Ludwig [1].
Primarily it refers to the process by which a temperature gradient drives
diffusion, known as the Soret effect [2], and the converse process, by which a
concentration gradient drives irreversible heat flux, called the Dufour effect
[3]. Rahman and Saghir provide a concise historical overview of
thermodiffusion research over the last century [4]. Interest continues today:
recently Platten discussed a host of new experimental techniques for probing
thermodiffusion [5]; Köhler and Morozov reviewed measurements for organic
solutions [6].

Chemical engineers have pursued separation processes based on thermodiffusion
for decades [7]. Chemical physicists, however, utilized the phenomenon first,
exploiting the Soret effect for isotope separation [8]. Thermodiffusion finds
contemporary applications in semiconductor fabrication [9, 10]. In
electrochemistry, thermoelectric phenomena such as the Peltier effect have
been tied to Soret and Dufour fluxes [11]. Metallurgists have investigated
thermodiffusion in liquid alloys [12]. Interest has also surged among
biochemical researchers, who study thermophoresis of biomolecules [13].

Gas-phase thermodiffusion has a sound theoretical foundation, laid by Chapman
and Enskog’s kinetic theory [14]. Hirschfelder, Curtiss, and Bird implemented
the kinetic theory for multicomponent gases, obtaining transport laws that
account for coupled heat flux and mass flux [15]. These and an alternative,
flux-explicit formulation are covered in the textbook of Bird, Stewart, and
Lightfoot [16]. A broader framework for coupled heat and mass transport was
created by Onsager [17], whose theory of irreversible thermodynamics applies
without modification to nonideal and condensed phases, as well as perfect
gases.

Structural differences among the available theoretical frameworks pose a
barrier to experimentalists, whose work has been constrained almost entirely
to thermodiffusion in binary systems. Importantly, properties measured in one
embodiment of the theory do not readily translate into another. The numerical
complexity of multicomponent simulations has also impeded progress.

Here we set out to unify the major continuum-scale theories of
thermodiffusion, and to deploy the elementary principles of irreversible
thermodynamics to produce more accurate and robust numerical methods. As a
practical matter, the unification effort reveals how to interconvert property
sets measured under different paradigms.

We propose a new constitutive scheme—the Onsager–Stefan–Maxwell
thermodiffusion equations—that simplifies analysis and numerics by accounting
for heat as a pseudo-species. This provides a natural route by which binary-
system properties can be extrapolated to simulate multicomponent systems. For
10 noble-gas pairs, we tabulate the properties needed to simulate binary
thermodiffusion with any of the popular theories and in addition detail how to
extrapolate these properties to simulate $n$-ary systems. We demonstrate the
applicability of the Onsager–Stefan–Maxwell framework by simulating steady,
three dimensional thermodiffusion in a ternary noble-gas mixture.

##  2 Approaches to thermodiffusion

The typical routes taken when creating thermodiffusion models are: derivation
by methods of irreversible thermodynamics; derivation from kinetic theory; and
intuitive formulation by extending Fick’s law and Fourier’s law. Each of these
choices may have particular benefits and disadvantages with regard to physical
rigour, ease of parameterization, and computational simplicity. Our task of
reconciling these different perspectives begins with summarizing the available
theories.

###  2.1 Irreversible thermodynamics

Transport laws based on irreversible thermodynamics derive from the local
entropy balance. Following the derivation of Jaumann [18], one begins with the
entropy continuity equation

$\frac{\partial\rho\hat{S}}{\partial
t}=-\vec{\nabla}\cdot\vec{N}_{S}+\dot{s},$ (2.1)

where $\rho$ is the local mass density of the phase, $\hat{S}$ the specific
entropy, $\vec{N}_{S}$ the total entropy flux, and $\dot{s}$ the volumetric
rate of entropy production. Through a lengthy thermodynamic analysis involving
the first law, the volumetric equation of state, Euler’s extensivity theorem,
and the Gibbs–Duhem equation, as well as various Legendre transformations and
Maxwell relations [15, 19], entropy continuity transforms into an expression
that describes the instantaneous energy dissipation. For an isothermal
viscous-fluid diffusion system at a null rate of strain, the dissipation rate
$T\dot{s}$ has units of power per volume,
$\mathrm{J}\,\mathrm{m}^{-3}\,\mathrm{s}^{-1}$, and takes the form

$T\dot{s}=-\vec{\nabla}\ln T\cdot\vec{q}\hskip
1.5pt^{\prime}+\sum_{i=1}^{n}\vec{d}_{i}\cdot\vec{v}_{i},$ (2.2)

in which $T$ is the absolute temperature. Here $\vec{d}_{i}$ denotes the
thermodynamic force that drives diffusion of species $i$, with units
$\mathrm{N}\,\mathrm{m}^{-3}$ (or $\mathrm{J}\,\mathrm{m}^{-4}$), and
$\vec{v}_{i}$ is its species velocity; $\vec{q}\hskip 1.5pt^{\prime}$ denotes
the irreversible heat flux, with units
$\mathrm{J}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}$.

To ensure consistency of units among the fluxes in the dissipation functional,
we define the _thermal velocity_ $\vec{v}_{0}$ as

$\vec{v}_{0}=\frac{\vec{q}\hskip 1.5pt^{\prime}}{\rho\hat{C}_{p}T},$ (2.3)

where $\hat{C}_{p}$ is the phase’s specific constant-pressure heat capacity,
in $\mathrm{J}\,\mathrm{k}\mathrm{g}^{-1}\,\mathrm{K}^{-1}$ (or
$\mathrm{m}^{2}\,\mathrm{s}^{-2}\,\mathrm{K}^{-1}$). This necessitates that
the conjugate _thermal driving force $\vec{d}_{0}$_ be

$\vec{d}_{0}=-\rho\hat{C}_{p}\vec{\nabla}T,$ (2.4)

simplified with the identity $T\vec{\nabla}\ln T=\vec{\nabla}T$. By
construction, $\vec{v}_{0}$ and $\vec{d}_{0}$ have units that match the
quantities $\vec{v}_{i}$ and $\vec{d}_{i}$ defined earlier. Henceforth we call
members of the collection $\\{\vec{d}_{i}\\}_{i=0}^{n}$ thermodiffusion
driving forces, and those of $\\{\vec{v}_{i}\\}_{i=0}^{n}$ thermodiffusive
velocities; we reserve the terms _diffusion driving forces_ and _species
velocities_ to refer to the smaller collections $\\{\vec{d}_{i}\\}_{i=1}^{n}$
and $\\{\vec{v}_{i}\\}_{i=1}^{n}$, respectively. A phase containing $n$
species affords $n+1$ thermodiffusive driving forces and $n+1$ thermodiffusive
velocities.

The principles of irreversible thermodynamics require that proper constitutive
laws for transport be formulated with the conjugate pairs of fluxes and forces
that appear in the dissipation functional [20]. In equation (2.2), for
instance, the thermodiffusive velocity $\vec{v}_{i}$ is conjugate to the
respective thermodiffusion driving force $\vec{d}_{i}$. Onsager postulated
that, within nonuniform systems not too far from equilibrium, global
thermodynamic state functions can be applied to describe the local states [17,
21]. Under this condition of ‘microscopic equilibrium’, Onsager posited
further that a linear transformation maps the thermodiffusive velocities into
their conjugate driving forces.

More concretely, thermodiffusion constitutive laws from irreversible
thermodynamics take the force-explicit forms

$\displaystyle\vec{d}_{0}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}\vec{v}_{0}+\sum\limits_{j=1}^{n}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0j}\vec{v}_{j}$
(2.5)
$\displaystyle\vec{d}_{i}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{i0}\vec{v}_{0}+\sum\limits_{j=1}^{n}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{ij}\vec{v}_{j},$
(2.6)

where the second equation holds for each $i\in\left\\{1,...,n\right\\}$.111For
simplicity we have assumed when writing equations (2.5) and (2.6) that the
material being modelled is isotropic, in which case the coefficients
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{ij}$
are scalars. In general, anisotropic systems can be handled by replacing every
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{ij}$
with a symmetric tensor [22], and replacing all scalar multiplications on the
right with tensor-vector products. These relationships, which we name the
_Onsager–Stefan–Maxwell thermodiffusion equations_ , introduce the Onsager
drag coefficients
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{ij}$,
a set of material properties with units
$\mathrm{J}\,\mathrm{s}\,\mathrm{m}^{-5}$ (i.e., force per volume per
velocity) that parameterizes how the thermodiffusive velocities map into the
thermodiffusive driving forces. These equations generalize the isothermal
Onsager–Stefan–Maxwell equations

$\vec{d}_{i}=\sum\limits_{j=1}^{n}{M}_{ij}\vec{v}_{j},$ (2.7)

for $i\in\\{1,\dots,n\\}$, expressed in terms of the isothermal drag
coefficients $M_{ij}$ [23]. Note that equations (2.5) and (2.6) have not been
written as compactly as possible because in later discussion it will sometimes
be convenient to consider the equation for the thermal force $\vec{d}_{0}$
separately, and to keep the coefficients whose second index equals $0$ out of
the sums.

It is often useful to think of the set of all Onsager drag coefficients as
being arrayed within an $(n+1)\times(n+1)$ matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
whose row and column indices range across the integers from $0$ to $n$, which
we will call the _non-isothermal drag matrix_. Onsager asserts that with a
proper choice of velocities and driving forces based in the energy
dissipation, such as equation (2.2), the reciprocal relation among the
transport coefficients is expressed by symmetry of
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
[17]. This symmetry has practical implications for the spectral structure of
the thermodiffusion problem.

Onsager’s formulation is not as abstract as it first appears. Hirschfelder et
al. substantiated the general concepts, manipulating thermodynamic laws and
hydrodynamic equations to cast the dissipation functional in terms of standard
quantities from thermodynamics [15, 19]. The general diffusion driving forces
within a non-isobaric, non-isothermal, multicomponent fluid are [15, 19]

$\vec{d}_{i}=-c_{i}\left(\vec{\nabla}\mu_{i}+\overline{S}_{i}\vec{\nabla}T-\frac{\overline{m}_{i}}{\rho}\vec{\nabla}p\right).$
(2.8)

In this definition, $p$ denotes the external pressure; for each species $i$,
$c_{i}$ is the molar concentration, $\overline{m}_{i}$ the molar mass,
$\mu_{i}$ the chemical potential, and $\overline{S}_{i}$ the partial molar
entropy. Thus practical implementation of laws (2.5) and (2.6) is possible if
the thermodynamic state variables are known functions of
$\left\\{c_{i}\right\\}_{i=1}^{n}$, $T$, and $p$. Determining these state
variables is sufficiently involved, however, that it is uncommon for
experimentalists to perform measurements sufficient to specify the drag matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
for condensed phases.

###  2.2 The kinetic theory of gases

Chapman–Enskog theory [14] leads to a set of thermodiffusion transport laws
whose structure differs from laws (2.5)-(2.6). Through the Chapman–Enskog
perturbation of the Boltzmann equation for a multicomponent gas, Hirschfelder
et al. [15] derived the so-called _generalized Stefan–Maxwell equations_ ,
which they wrote as

$\vec{d}_{i}-\sum_{j=1}^{n}\frac{RTc_{i}c_{j}}{\mathscr{D}_{ij}c_{\text{T}}}\left(\frac{D_{i}^{T}}{\rho_{i}}-\frac{D_{j}^{T}}{\rho_{j}}\right)\vec{\nabla}\ln
T=\sum_{j=1}^{n}\frac{RTc_{i}c_{j}}{\mathscr{D}_{ij}c_{\text{T}}}\left(\vec{v}_{i}-\vec{v}_{j}\right)$
(2.9)

for $i\in\left\\{1,...,n\right\\}$, where $R$ is the gas constant. Here
$c_{\textrm{T}}=\sum_{i=1}^{n}c_{i}$ is the total molar concentration (inverse
molar volume) and $\rho_{i}$ denotes the mass density of species $i$;
$\rho_{i}=c_{i}\overline{m}_{i}$. Equation (2.9) introduces another set of
transport properties, namely, the Stefan–Maxwell diffusion coefficients
$\mathscr{D}_{ij}$ and the coefficients of thermal diffusion $D_{i}^{T}$.

As pointed out by Newman [11], the parameters $D_{i}^{T}$ have units of
viscosity (_i.e._ $\mathrm{P}\mathrm{a}\,\mathrm{s}$, or
$\mathrm{J}\,\mathrm{s}\,\mathrm{m}^{-3}$, or
$\mathrm{k}\mathrm{g}\,\mathrm{m}^{-1}\,\mathrm{s}^{-1}$), which is
unconventional for a diffusion coefficient. To retain dimensional consistency
and simplify notation in the sequel, we therefore define

$\mathscr{D}_{i}^{T}=\frac{D_{i}^{T}}{\rho_{i}},$ (2.10)

which we call the _Soret diffusivity_ of species $i$. The Soret diffusivity
$\mathscr{D}_{i}^{T}$ has typical units, $\mathrm{m}^{2}\,\mathrm{s}^{-1}$.

Some attempts have been made to bring constitutive formulation (2.9) into
harmony with (2.5)-(2.6). For isothermal multicomponent diffusion, Lightfoot
et al. [23] pointed out that the Stefan–Maxwell equations can be understood in
terms of Onsager’s transport equations by identifying the entries $M_{ij}$ of
the isothermal transport matrix M as

$M_{ij}=\begin{cases}-\frac{RTc_{i}c_{j}}{\mathscr{D}_{ij}c_{\text{T}}}\>\>\>\text{if}\>\>i\neq
j\\\ \sum_{k\neq
i}^{n}\frac{RTc_{i}c_{k}}{\mathscr{D}_{ik}c_{\text{T}}}\>\>\>\text{if}\>\>i=j.\end{cases}$
(2.11)

Onsager reciprocity therefore implies symmetry of the coefficients
$\mathscr{D}_{ij}$.

Identity (2.11) fails to carry over to the non-isothermal case because the
generalized Stefan–Maxwell equations (2.9) mix driving forces and fluxes,
which may also violate the conditions required to ensure a symmetric Onsager
reciprocal relation [24, 25]. Specifically, equation (2.9) contains no thermal
velocity on the right; on the left, it also mixes the thermal driving force
with a diffusion driving force. Therefore
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}\neq\textbf{M}$,
even if comparison is restricted to the submatrix with $1\leq i,j\leq n$, as
we shall see in section 3.

###  2.3 Multicomponent Fick’s law

Many transport analyses begin in a more _ad hoc_ fashion, by postulating flux-
explicit forms of the transport laws with the structure of Fick’s law [16, 26,
20]. For a general multicomponent mixture this results in expressions for the
excess molar flux of species $i$, $\vec{J}_{i}$, of the form

$\vec{J}_{i}=-\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{D}}_{i0}\vec{\nabla}T-\sum_{j=1}^{n}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{D}}_{ij}\vec{\nabla}c_{j}$
(2.12)

for $i\in\left\\{1,...,n\right\\}$. To close this thermodiffusion model,
Fickian laws from equation (2.12) must be augmented by a complementary
extension of Fourier’s law of heat conduction, written as

$\vec{q}\hskip
1.5pt^{\prime}=-\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{D}}_{00}\vec{\nabla}T-\sum_{j=1}^{n}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{D}}_{0j}\vec{\nabla}c_{j},$
(2.13)

which expresses the irreversible heat flux.

Clearly the entries of the $(n+1)\times(n+1)$ Fickian thermodiffusion matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{D}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{D}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{D}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{D}}}$
have inconsistent units. In the multicomponent generalization of Fick’s law,
equation (2.12), the transport coefficients
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{D}}_{ij}$
that accompany concentration gradients can be seen as traditional
diffusivities, in $\mathrm{m}^{2}\,\mathrm{s}^{-1}$; the coefficients of
thermal diffusion
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{D}}_{i0}$
have units
$\mathrm{m}\mathrm{o}\mathrm{l}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}\,\mathrm{K}^{-1}$.
In the generalization of Fourier’s law from equation (2.13),
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{D}}_{00}$
is an effective thermal conductivity, in
$\mathrm{J}\,\mathrm{m}^{-1}\,\mathrm{s}^{-1}\,\mathrm{K}^{-1}$, but the other
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{D}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{D}}_{0j}$
are Dufour coefficients, in
$\mathrm{J}\,\mathrm{m}^{2}\,\mathrm{s}^{-1}\,{mol}^{-1}$. These discrepancies
emphasize that the reciprocal relation here is ambiguous; it cannot simply be
a symmetry of the
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{D}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{D}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{D}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{D}}}$
matrix.

A few paths have been proposed to connect the Onsager–Stefan–Maxwell equations
(2.5)-(2.6) with relations in the Fickian forms of equations (2.12)-(2.13).
Despite their linearity, transformation between force-explicit and flux-
explicit constitutive formulations requires more than a simple inversion of
the transport matrix. As pointed out by Helfand [27], inversion is thwarted by
the fact that a Gibbs–Duhem relation constrains the diffusion driving forces.
When $\vec{d}_{i}$ is defined with equation (2.8), the Gibbs–Duhem equation
requires that

$\sum_{i=1}^{n}\vec{d}_{i}=0,$ (2.14)

which further implies singularity of the isothermal drag matrix M. A pseudo-
inversion process is enabled, however, by considering kinematic constraints
implied by the invariance of diffusion phenomena with respect to the velocity
that drives convection.

Flux-law inversions implicitly based on the mass-average velocity were
addressed by Curtiss and Bird [7] and Fong et al. [28]. Inversions with
respect to other reference velocities have been implemented by Helfand [27]
and Newman [29]. Even when applied in isothermal cases, however, these
derivation procedures are somewhat impromptu. The Curtiss–Bird approach to
inversion, which produces the most well known flux-explicit model for
thermodiffusion, was only taken up to quaternary mixtures. No general method
so far has been put forward to derive a flux-explicit form of the
Onsager–Stefan–Maxwell equations relative to an arbitrary convective velocity.

##  3 Properties of the non-isothermal drag matrix

First we set out to reconcile the generalized Stefan–Maxwell equations (2.9)
with the Onsager–Stefan–Maxwell thermodiffusion equations (2.5)-(2.6). We then
discuss how the second law of thermodynamics and the principle of convection
invariance impose further structure beyond simple symmetry of the non-
isothermal drag matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$.

###  3.1 Construction and connection to kinetic theory

One can move toward the form of the generalized Stefan–Maxwell equations (2.9)
by isolating $\vec{v}_{0}$ from equation (2.5) and substituting the result
into equation (2.6). Regrouping terms and rearranging yields

$\vec{d}_{i}+\frac{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0i}\rho\hat{C}_{p}T}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\vec{\nabla}\ln
T=\sum_{j=1}^{n}\left(\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{ij}-\frac{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0i}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0j}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\right)\vec{v}_{j}.$
(3.1)

Agreement of the left side with (2.9) requires that

$\frac{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0i}\rho\hat{C}_{p}T}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}=-\sum_{j=1}^{n}\frac{RTc_{i}c_{j}}{c_{\text{T}}\mathscr{D}_{ij}}\left(\mathscr{D}_{i}^{T}-\mathscr{D}_{j}^{T}\right).$
(3.2)

Insertion of (2.11), followed by the simplification $\sum_{j=1}^{n}M_{ij}=0$
implied by Gibbs–Duhem equation (2.14) and the Onsager reciprocal relation,
produces

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0i}=-\frac{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}{\rho\hat{C}_{p}T}\sum_{j=1}^{n}M_{ij}\mathscr{D}_{j}^{T}.$
(3.3)

Similarly, agreement of the right side of equation (3.1) with (2.9) demands
for all $1\leq i,j\leq n$ that

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{ij}=M_{ij}+\frac{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0i}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0j}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}},$
(3.4)

where $M_{ij}$ relates to Stefan–Maxwell diffusivities through equation
(2.11). These conclusions can be summarized by writing the non-isothermal drag
matrix in block form, as

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}=\begin{bmatrix}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}&\widetilde{\textbf{m}}_{0}^{\top}\\\
\widetilde{\textbf{m}}_{0}&\textbf{M}+\dfrac{\widetilde{\textbf{m}}_{0}\widetilde{\textbf{m}}_{0}^{\top}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\end{bmatrix},$
(3.5)

in which the column matrix
$\widetilde{\textbf{m}}_{0}=\left[\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0i}\right]^{i=1}_{n}$
contains entries
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0i}$
defined by equation (3.3), superscript $\top$ indicates the matrix transpose,
and the first row’s entries follow from the reciprocal relation
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}^{\top}$.
This construction affirms that the generalized Stefan–Maxwell equations (2.9)
derived from kinetic theory are compatible with the Onsager–Stefan–Maxwell
thermodiffusion laws (2.5)-(2.6) from irreversible thermodynamics.

Despite the mixing of fluxes with driving forces in the generalized
Stefan–Maxwell equations, one can show that their structure is consistent with
the reciprocal relation
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}=\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}^{\top}$.
Substituting matrix (3.2) into equation (3.5), using flux law (2.5), then
returning to the original thermal flux and driving force with definitions
(2.3)-(2.4), yields

$\vec{q}\hskip
1.5pt^{\prime}=-\frac{\rho^{2}\hat{C}_{\rho}^{2}T}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\nabla
T-\sum_{i,j=1}^{n}\frac{RTc_{i}c_{j}}{c_{\text{T}}\mathscr{D}_{ij}}\left(\mathscr{D}_{i}^{T}-\mathscr{D}_{j}^{T}\right)\vec{v}_{j}.$
(3.6)

This precisely matches the expression for irreversible heat flux provided by
Hirschfelder, Curtiss, and Bird [15], although they define the coefficient of
$-\vec{\nabla}T$ as the thermal conductivity. For now we retain the symbol
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}$,
leaving a detailed discussion of heat conduction to section 4.5.

###  3.2 Entropy generation and spectral structure.

A benefit of casting the generalized Stefan–Maxwell equations into the
Onsager–Stefan–Maxwell form is the structure inherited by the transport matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$.
Given constitutive laws (2.5) and (2.6), the entropy production due to
thermodiffusion is equivalent to

$T\dot{s}=\sum_{i,j=0}^{n}\vec{v_{i}}\cdot\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{ij}\vec{v}_{j}.$
(3.7)

The second law of thermodynamics requires that $\dot{s}\geq 0$ for every
collection of thermodiffusive velocities, so the matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
must be positive semidefinite, as well as symmetric.

Direct calculation shows a connection between the spectral structures of the
isothermal and non-isothermal drag matrices, M and
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$,
respectively. Insertion of matrix (3.5) into equation (3.7), algebraic
rearrangement of terms involving the symmetry of
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$,
and elimination of $\vec{v}_{0}$ with equation (2.5) show that

$T\dot{s}=\frac{1}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\big{\|}\vec{d}_{0}\big{\|}^{2}+\sum_{i,j=1}^{n}\vec{v}_{i}\cdot
M_{ij}\vec{v}_{j},$ (3.8)

in which $\|\vec{u}\|^{2}$ = $\vec{u}\cdot\vec{u}$ for any vector $\vec{u}$.
Positive semidefiniteness of the non-isothermal matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
is thus implied by positive semidefiniteness of isothermal matrix M, so long
as the Onsager drag
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}$
is positive. It is noteworthy that the sum on the right of equation (3.8) can
be expressed in terms of Stefan–Maxwell coefficients as [30]

$\sum_{i,j=1}^{n}\vec{v}_{i}\cdot
M_{ij}\vec{v}_{j}=\frac{RT}{2c_{\text{T}}}\sum_{i,j=1}^{n}\frac{c_{i}c_{j}}{\mathscr{D}_{ij}}\left\|\vec{v}_{i}-\vec{v}_{j}\right\|^{2}.$
(3.9)

Thus, if all of the Stefan–Maxwell diffusivities in a material are constant,
they must also all be positive. One should be wary that for materials where
$\mathscr{D}_{ij}$ varies with composition, negative values are possible, and
have been observed.

The distinction between convection and diffusion requires that no dissipation
occurs unless there is relative motion of the species [21]. This feature
naturally arises from the generalized Stefan–Maxwell equations, as may be seen
from equation (3.9) by noting that
$\vec{v}_{i}-\vec{v}_{j}=\left(\vec{v}_{i}+\vec{u}\right)-\left(\vec{v}_{j}+\vec{u}\right)$
for any vector $\vec{u}$. In the isothermal case, invariance of the
dissipation with respect to convection can equivalently be stated as a
condition that M has a single null eigenvalue corresponding to the eigenvector
$\textbf{1}=[1,1,...,1]^{\top}\in\mathbb{R}^{n}$, as noted by Helfand [27].
This spectral structure is also necessary to enforce Gibbs–Duhem equation
(2.14). Indeed, in the isothermal case, equation (2.6) reduces to

$\vec{d}_{i}=\sum_{j=1}^{n}M_{ij}\vec{v}_{j},$ (3.10)

and the facts that $\textbf{M}=\textbf{M}^{\top}$ and 1 is a null eigenvector
of M imply equality (2.14).

A similar conclusion about spectral structure can be drawn for the non-
isothermal drag matrix. One can use equation (3.3) to show that the sum over
rows $1$ through $n$ of the first ($j=0$) column of
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
satisfies

$\sum_{i=1}^{n}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{i0}=\sum_{i,j=1}^{n}M_{ij}\mathscr{D}_{j}^{T}=0.$
(3.11)

Put another way, this says the column matrix 1 is orthogonal to the previously
defined column matrix $\widetilde{\textbf{m}}_{0}$, so that
$\textbf{1}^{\top}\widetilde{\textbf{m}}_{0}=0$. Consequently, the block
matrix from equation (3.5) must satisfy the matrix equation

$\begin{bmatrix}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}&\widetilde{\textbf{m}}_{0}^{\top}\\\
\widetilde{\textbf{m}}_{0}&\textbf{M}+\dfrac{\widetilde{\textbf{m}}_{0}\widetilde{\textbf{m}}_{0}^{\top}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\end{bmatrix}\begin{bmatrix}0\\\
\textbf{1}\end{bmatrix}=0\cdot\begin{bmatrix}0\\\ \textbf{1}\end{bmatrix},$
(3.12)

demonstrating that the the column matrix
$[\begin{array}[]{cc}0&\textbf{1}^{\top}\end{array}]^{\top}$ is a null
eigenvector of the non-isothermal drag matrix. Like the isothermal case, this
spectral structure of
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
ensures the enforcement of Gibbs–Duhem relation (2.14).

##  4 Inverting the constitutive laws

Many practitioners find it desirable to implement a flux-explicit (or
velocity-explicit) transport formulation. As mentioned before, the
transformation of force-explicit laws to a flux-explicit form is not entirely
straightforward: since the isothermal drag matrix satisfies
$\textbf{M}\textbf{1}=0\textbf{1}=\textbf{o}$, where o indicates a column of
zeroes, it is singular, and consequently equation (3.10) cannot be inverted
directly. To resolve this issue it is necessary to consider how the flux-
explicit laws behave with regard to the reference velocity selected to
describe convection. Once the process of choosing a reference velocity is
formalized, a kinematic relation derived from that choice can be exploited to
develop a unique inverted formulation of the force-explicit transport laws.

###  4.1 Convective reference velocities

Convection is fundamentally anchored in species fluxes: a valid choice of
convective velocity must always be a linear combination of the species
velocities. Newman made the simplest choice, employing the $k^{\textrm{th}}$
species’ velocity $\vec{v}_{k}$ as a reference for convection [29]. Other
common reference velocities include the mole-average velocity ${\vec{v}}\hskip
1.0pt^{c}$ employed by Waldmann and Van Der Valk [31, 32],

$\vec{v}\hskip
1.0pt^{c}=\frac{1}{c_{\textrm{T}}}\sum_{j=1}^{n}c_{j}\vec{v}_{j},$ (4.1)

and the mass-average velocity ${\vec{v}}\hskip 1.0pt^{\rho}$, defined as

$\vec{v}\hskip 1.0pt^{\rho}=\frac{1}{\rho}\sum_{j=1}^{n}\rho_{j}\vec{v}_{j},$
(4.2)

which was used by both Hirschfelder et al. [15] and Curtiss and Bird [7].

In Newman’s inversion process [29], after species $k$ is selected as the
reference for convection, it is shown that the isothermal drag matrix M can be
truncated by deleting its $k^{\textrm{th}}$ row and column without loss of
information. This first minor matrix is nonsingular and can be inverted
directly. Helfand considered isothermal, isobaric multicomponent diffusion
[27]; when writing the Onsager–Stefan–Maxwell laws, he used chemical-potential
driving forces $\vec{\nabla}\mu_{i}$ in place of $\vec{d}_{i}$ and molar
fluxes $c_{i}\vec{v}_{i}$ in place of $\vec{v}_{i}$. These choices make the
flux-explicit coefficient matrix share the same null eigenvector as the force-
explicit matrix, enabling a Moore–Penrose pseudo-inversion. Hirschfelder et
al. [15] and Curtiss and Bird [7] extensively manipulate the generalized
Stefan–Maxwell equations to derive velocity-explicit thermodiffusion laws; the
rather convoluted process they lay out is summarized in Bird, Stewart, and
Lightfoot’s book [16]. More recently Bothe et al. [33] compared the isothermal
Stefan–Maxwell equations with a generalized form of Fick’s law and showed
equivalent properties such as positivity preservation between the two
approaches.

Depending on the context, different reference velocities may be desirable in
flux-explicit constitutive laws [16, 34]. Here we write a general reference
velocity as $\vec{v}\hskip 1.0pt^{\psi}$, following Goyal and Monroe [19]. To
create a convective velocity, one selects a collection of thermodynamic
extensive quantities per unit volume $\left\\{\psi_{i}\right\\}_{i=1}^{n}$,
whose total volumetric amount $\psi_{\textrm{T}}$ is defined as

$\psi_{\textrm{T}}=\sum_{i=1}^{n}\psi_{i}=\textbf{1}^{\top}\bm{\psi},$ (4.3)

in which the rightmost expression introduces the column matrix
$\bm{\psi}=\left[\psi_{i}\right]^{i=1}_{n}$. Note that the quantities
$\psi_{i}$ are not limited to choices with these precise physical
interpretations; the only necessary restrictions on $\bm{\psi}$ are that it is
real and that $\psi_{\textrm{T}}\neq 0$ to avoid degeneracy.

The $\psi$-average velocity, $\vec{v}\hskip 1.0pt^{\psi}$, is defined as

$\vec{v}\hskip
1.0pt^{\psi}=\frac{1}{\psi_{T}}\sum_{j=1}^{n}\psi_{j}\vec{v}_{j}.$ (4.4)

To further define the excess velocity of species $i$ relative to the
$\psi$-average velocity, $\vec{v}_{i}^{\hskip 1.0pt\psi}$, take

$\vec{v}_{i}^{\hskip 1.0pt\psi}=\vec{v}_{i}-\vec{v}\hskip 1.0pt^{\psi}.$ (4.5)

By design, the collection of these excess velocities satisfies a kinematic
relation

$\sum_{j=1}^{n}\psi_{j}\vec{v}\hskip 1.0pt^{\psi}_{j}=0.$ (4.6)

In the language of linear algebra, this says that the space of excess species
velocities is orthogonal to $\bm{\psi}$. Thus we may also refer to the column
$\bm{\psi}$ itself as the kinematic relation associated with convective
velocity $\vec{v}\hskip 1.0pt^{\psi}$, because it defines the orthogonal
subspace of the species velocities in which the set of excess species
velocities must reside.

###  4.2 Inversion of the isothermal drag matrix

Invariance of the isothermal drag matrix M with respect to the choice of
convective velocity requires that

$\vec{d}_{i}=\sum_{j=1}^{n}M_{ij}\vec{v}_{j}=\sum_{j=1}^{n}M_{ij}\vec{v}_{j}^{\hskip
1.0pt\psi}$ (4.7)

for any valid kinematic relation $\bm{\psi}$. Bearing in mind equation (4.6),
one can augment the transport matrix and write equation (4.7) equivalently as

$\vec{d}_{i}=\sum_{j=1}^{n}M_{ij}\vec{v}_{j}^{\hskip
1.0pt\psi}+\gamma\psi_{i}\sum_{j=1}^{n}\psi_{j}\vec{v}_{j}^{\hskip
1.0pt\psi}=\sum_{j=1}^{n}M^{\psi}_{ij}\left(\gamma\right)\vec{v}_{j}^{\hskip
1.0pt\psi}$ (4.8)

for any $\gamma>0$. On the right we have defined the _augmented drag matrix_
with respect to the kinematic relation $\bm{\psi}$ as

$M^{\psi}_{ij}\left(\gamma\right)=M_{ij}+\gamma\psi_{i}\psi_{j},~{}~{}\textrm{or}~{}~{}\textbf{M}^{\psi}\left(\gamma\right)=\textbf{M}+\gamma\bm{\psi}\bm{\psi}^{\top}.$
(4.9)

This augmentation was first considered in the context of multicomponent
transport by Helfand [27]; much later it was developed independently by
Giovangigli [35, 36, 37].

Significantly, the augmented transport matrix
$\textbf{M}^{\psi}\left(\gamma\right)$ is positive definite [30], and is
consequently invertible. Upon inversion one finds that the excess velocity of
$i$ is

$\vec{v}_{i}^{\hskip
1.0pt\psi}=\sum_{j=1}^{n}\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}_{ij}\vec{d}_{j},$
(4.10)

the desired flux-explicit form of the transport equations. In its present
state this formulation is not satisfactory, however. First, the ‘penalty term’
involving $\gamma$ is artificial. Second, and of deeper physical importance,
the kinematic relation from equation (4.6) is not obviously enforced for every
collection of driving forces.

An unambiguous formulation of inverted transport laws is found by defining

$\textbf{L}^{\psi}=\lim_{\gamma\rightarrow\infty}\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}.$
(4.11)

This limit exists for any valid choice of M and $\bm{\psi}$. Indeed,

$\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}=\textbf{L}^{\psi}+\mathscr{O}(\gamma^{-1}),$
(4.12)

and the unique matrix $\textbf{L}^{\psi}$, which we call the _Onsager
diffusivity matrix_ relative to the $\psi$-average velocity, is symmetric
whenever M is. This argument is substantiated in A, which also provides two
algebraic formulas useful to compute $\textbf{L}^{\psi}$ from M without
resorting to a limit process.

Since equation (4.10) holds for any nonzero $\gamma$, one can derive the flux-
explicit form of the transport laws by calculating the limit of equation
(4.10) as $\gamma\to\infty$, giving

$\vec{v}^{\hskip 1.0pt\psi}_{i}=\sum_{j=1}^{n}L^{\psi}_{ij}\vec{d}_{j}$ (4.13)

as a consequence of the asymptotic behavior in equation (4.12). The matrix
$\textbf{L}^{\psi}$ so formed also satisfies

$\sum_{j=1}^{n}L^{\psi}_{ij}\psi_{j}=0,$ (4.14)

that is to say, the kinematic relation $\bm{\psi}$ is a null eigenvector of
$\textbf{L}^{\psi}$. Analogous to the spectral structure of M enforcing the
Gibbs–Duhem equation, this spectral structure of $\textbf{L}^{\psi}$ enforces
kinematic constraint (4.6).

The procedure outlined above encompasses all methods currently used to invert
isothermal Onsager–Stefan–Maxwell equations. Newman’s inverse, relative to
$\vec{v}_{k}$, is produced by choosing a kinematic relation such that
$\psi_{i}=\delta_{ik}$, where $\delta_{ij}$ is the Kronecker delta. Given the
diffusion driving forces in equation (2.8), which make 1 the null eigenvector
of M, Helfand’s inverse is formed by taking $\bm{\psi}=\textbf{1}$, an
unorthodox choice where the arithmetic mean of the species velocities drives
convection. Curtiss and Bird’s inverse derives from the kinematic relation for
the mass-average velocity, $\psi_{i}=\rho_{i}$. Crucially, the need for
kinematic relations here indicates that one cannot consider a generic matrix L
in a flux-explicit thermodiffusion model. Instead there is a class of equally
valid Onsager diffusivity matrices, from which one is determined by specifying
a convective velocity.

Sometimes it may be of interest to calculate M given an $\textbf{L}^{\psi}$,
as discussed by Monroe and Newman in a specific case [25]. This can be
accomplished generally by augmenting the singular matrix $\textbf{L}^{\psi}$
to make

$\textbf{L}^{\psi}\left(\gamma\right)=\textbf{L}^{\psi}+\gamma\textbf{1}\textbf{1}^{\top},$
(4.15)

again with $\gamma>0$. Then, similar to the process given above, one can pass
into a limit to get

$\textbf{M}=\lim_{\gamma\rightarrow\infty}\left[\textbf{L}^{\psi}\left(\gamma\right)\right]^{-1}.$
(4.16)

The M form may be preferable when reporting experimental measurements, since
it is independent of the convective velocity used for data processing.
Although there exist many legitimate $\textbf{L}^{\psi}$ matrices, there is
just one valid choice of M for a proper set of diffusion driving forces.

###  4.3 Inversion of the non-isothermal drag matrix

Problems posed by the convection invariance of M carry over into the non-
isothermal drag matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$:
the null eigenvector
$[\begin{array}[]{cc}0&\textbf{1}^{\top}\end{array}]^{\top}$ induced by
convection invariance implies that
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
cannot be directly inverted. It is, however, straightforward to adapt the
pseudo-inversion procedure developed earlier to the non-isothermal case.

Thermal properties in the generalized Stefan–Maxwell equations merit some
additional discussion here. Observe that only differences between Soret
diffusivities appear in equation (2.9). Since the physical content of that
equation is invariant to a variable change where the same scalar is added to
every Soret diffusivity, it follows that the collection
$\left\\{\mathscr{D}_{i}^{T}\right\\}_{i=1}^{n}$ must be linearly dependent.
To quantify Soret diffusivities in isolation, one must adopt an additional
constraint to close this degree of freedom. Newman circumvented the ambiguity
by defining a new property,

$\mathscr{A}_{ij}=\frac{D_{i}^{T}}{\rho_{i}}-\frac{D_{j}^{T}}{\rho_{j}}=\mathscr{D}_{i}^{T}-\mathscr{D}_{j}^{T},$
(4.17)

which we call the _Newman–Soret diffusivity_. As well as having units
$\mathrm{m}^{2}\,\mathrm{s}^{-1}$, the fact that $\mathscr{A}_{ii}=0$ ensures
that the nontrivial $\mathscr{A}_{ij}$ are linearly independent for a given
$j$. Newman further asserts that these diffusivities are more approximately
constant with respect to composition in binary systems [11]. Newman’s choice
of symbol emphasizes antisymmetry in the indices, that is, that
$\mathscr{A}_{ji}=-\mathscr{A}_{ij}$.

It is standard to enforce a kinematic relation to close the degree of freedom
when specifying the coefficients of thermal diffusion. Here we impose the
general constraint

$\sum_{j=1}^{n}\psi_{j}\mathscr{D}_{j}^{T}=0,$ (4.18)

which is agnostic to the choice of convective velocity. Given any complete set
of $\mathscr{A}_{ij}$ with $j$ fixed, the identities
$\mathscr{A}_{kl}=\mathscr{A}_{kj}-\mathscr{A}_{lj}$ and
$\mathscr{A}_{lk}=-\mathscr{A}_{kl}$ fill out the whole matrix of Newman–Soret
diffusivities. Thus, any particular Soret diffusivity can be computed with the
formula

$\mathscr{D}_{i}^{T}=\frac{1}{\psi_{\textrm{T}}}\sum_{k=1}^{n}\psi_{k}\mathscr{A}_{ik},$
(4.19)

so long as one row or column of the Newman–Soret matrix is known. Given the
kinematic relation $\bm{\psi}$, equations (4.17) and (4.19) establish a
bijective map between the Soret diffusivities and the Newman–Soret matrix.

In their implementation of Chapman–Enskog theory, Hirschfelder et al. state
[15] that

$\sum_{i=1}^{n}D_{i}^{T}=0.$ (4.20)

This implies that the properties $D_{i}^{T}/\rho_{i}$ satisfy a kinematic
relation $\psi_{i}=\rho_{i}$; the Soret diffusivities stand relative to the
mass-average velocity. Brenner argued that some constitutive laws are more
appropriately formulated with respect the volume-average velocity, however
[38]. Rather than advocating one view or the other, we adopt Newman–Soret
diffusivities here, which emphasizes that thermal diffusion can be fully
parameterized without choosing a reference velocity for convection.

One can write a new equivalent of equation (3.3) in terms of the augmented
drag matrix $\textbf{M}^{\psi}\left(\gamma\right)$ from equation (4.9), as

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0i}^{\psi}\left(\gamma\right)=-\frac{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}{\rho\hat{C}_{p}T}\sum_{j=1}^{n}M^{\psi}_{ij}\left(\gamma\right)\mathscr{D}_{j}^{T}.$
(4.21)

These coefficients can be assembled into an augmented column matrix
$\widetilde{\textbf{m}}_{0}^{\psi}\left(\gamma\right)=[\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{0i}^{\psi}\left(\gamma\right)]^{i=1}_{n}$.
The constraint stated in equation (4.18) is natural because it implies that

$\widetilde{\textbf{m}}_{0}^{\psi}\left(\gamma\right)=\widetilde{\textbf{m}}_{0}$
(4.22)

for any $\gamma>0$. Thus augmentation of the non-isothermal drag matrix has no
impact on the thermal entries of
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$.

The non-isothermal transport matrix augmented by kinematic relation
$\bm{\psi}$ is defined as

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}^{\psi}\left(\gamma\right)=\begin{bmatrix}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}&\widetilde{\textbf{m}}_{0}^{\top}\\\
\widetilde{\textbf{m}}_{0}&\textbf{M}^{\psi}\left(\gamma\right)+\frac{\widetilde{\textbf{m}}_{0}\widetilde{\textbf{m}}_{0}^{\top}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\end{bmatrix},$
(4.23)

in which only the lower-right block depends on the penalty factor $\gamma$.
The inverse of this augmented matrix is partially calculable through Schur’s
complement. In block form,

$\left[\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}^{\psi}\left(\gamma\right)\right]^{-1}=\begin{bmatrix}\frac{1}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}+\frac{\widetilde{\textbf{m}}_{0}^{\top}\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}\widetilde{\textbf{m}}_{0}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}^{2}}&-\frac{\widetilde{\textbf{m}}_{0}^{\top}\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\\\
-\frac{\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}\widetilde{\textbf{m}}_{0}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}&\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}\end{bmatrix}.$
(4.24)

Interestingly, the only matrix inversion needed to carry out this computation
is that of the augmented isothermal drag matrix.

It is straightforward to rephrase equation (4.21) in a form that computes
Soret diffusivities from the Onsager drag coefficients. The construction

$\left[\mathscr{D}_{i}^{T}\right]^{i=1}_{n}=-\frac{\rho\hat{C}_{p}T}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}\widetilde{\textbf{m}}_{0}$
(4.25)

retrieves a column of Soret diffusivities constrained by equation (4.18). With
equation (4.25), one can define for $i\in\left\\{1,...,n\right\\}$ that

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}_{0i}^{\psi}=\frac{\mathscr{D}_{i}^{T}}{\rho\hat{C}_{p}T},$
(4.26)

and introduce a column
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}^{\hskip
1.0pt\psi}_{0}=[\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}_{0i}^{\psi}]^{i=1}_{n}$.
(The superscript $\psi$ on these quantities emphasizes kinematic relationship
(4.18), i.e. that
$\bm{\psi}^{\top}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}^{\hskip
1.0pt\psi}_{0}=0$.) Then

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}^{\hskip
1.0pt\psi}_{0}=-\frac{\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}\widetilde{\textbf{m}}_{0}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}$
(4.27)

for any $\gamma>0$, and the inverse augmented drag matrix transforms to

$\left[\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}^{\psi}\left(\gamma\right)\right]^{-1}=\begin{bmatrix}\frac{1}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\left(1-\widetilde{\textbf{m}}_{0}^{\top}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}^{\hskip
1.0pt\psi}_{0}\right)&\left(\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}^{\hskip
1.0pt\psi}_{0}\right)^{\top}\\\
\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}^{\hskip
1.0pt\psi}_{0}&\left[\textbf{M}^{\psi}\left(\gamma\right)\right]^{-1}\end{bmatrix}.$
(4.28)

Finally, the unique inverse of the non-isothermal drag matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
with respect to convective velocity $\vec{v}^{\hskip 1.0pt\psi}$ is found by
taking the limit of equation (4.28) as $\gamma\to\infty$. This produces the
_Onsager thermodiffusivity matrix_

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{L}}}^{\psi}=\begin{bmatrix}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}^{-1}\left(1-\widetilde{\textbf{m}}_{0}^{\top}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}^{\hskip
1.0pt\psi}_{0}\right)&\left(\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}^{\hskip
1.0pt\psi}_{0}\right)^{\top}\\\
\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}^{\hskip
1.0pt\psi}_{0}&\textbf{L}^{\psi}\end{bmatrix},$ (4.29)

a symmetric matrix with null eigenvector
$[\begin{array}[]{cc}0&\bm{\psi}^{\top}\end{array}]^{\top}$.

###  4.4 Onsager–Fick–Fourier laws

Equation (4.29) gives rise to a velocity-explicit form of the non-isothermal
flux laws similar to isothermal equation (4.13). A thermodynamically
consistent form of the generalized Fick’s law is given by

$\vec{v}_{i}^{\hskip
1.0pt\psi}=\sum_{j=0}^{n}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}^{\psi}_{ij}\vec{d}_{j}.$
(4.30)

Here it should be understood that $\vec{v}_{0}^{\hskip
1.0pt\psi}=\vec{v}_{0}$, because the thermal velocity $\vec{v}_{0}$ is not
included in the kinematic constraint from equation (4.6). Unlike the intuitive
extension of Fick’s law involving transport matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{D}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{D}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{D}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{D}}}$,
the derivation of
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{L}}}^{\hskip
1.0pt\psi}$ from the Onsager–Stefan–Maxwell thermodiffusion equations ensures
that the coefficient matrix in equation (4.30) has consistent units, as well
as being symmetric positive semidefinite and satisfying a kinematic
constraint.

Equation (4.30) substantiates a thermodynamically consistent set of
generalized Fourier–Fick thermodiffusion laws as follows. After inserting
equations (2.3), (2.4), and (4.26) into equation (4.30), and defining the
excess flux of species $i$ relative to the $\psi$-average velocity as

$\vec{J}_{i}^{\hskip 1.0pt\psi}=c_{i}\vec{v}_{i}^{\hskip 1.0pt\psi},$ (4.31)

one obtains flux-explicit constitutive equations

$\displaystyle\vec{q}\hskip 1.5pt^{\prime}$
$\displaystyle=-\rho^{2}\hat{C}_{p}^{2}T\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}_{00}^{\psi}\vec{\nabla}T+\sum_{j=1}^{n}\mathscr{D}_{j}^{T}\vec{d}_{j},$
(4.32) $\displaystyle\vec{J}_{i}^{\hskip 1.0pt\psi}$
$\displaystyle=-c_{i}\mathscr{D}_{i}^{T}\vec{\nabla}\ln
T+\sum_{j=1}^{n}c_{i}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}_{ij}^{\psi}\vec{d}_{j},$
(4.33)

for $i\in\left\\{1,...,n\right\\}$. Equation (4.32) is analogous to Fourier’s
law, with extra terms describing the Dufour effect; equations (4.33)
generalize Fick’s law, with an additional term for the Soret effect. These
equations satisfy Onsager’s principles of irreversible thermodynamics, in that
they involve all the proper fluxes and driving forces that appear in the
dissipation functional, as well as including a positive-semidefinite transport
matrix that satisfies a symmetric reciprocal relation and has a single null
eigenvector to enforce convection invariance. To emphasize their thermodynamic
rigour, we henceforth call equations (4.32) and (4.33) the
_Onsager–Fick–Fourier equations_.

###  4.5 Perspectives on thermal conductivity

To elucidate simple heat conduction, it is worth returning to the dissipation
functional (2.2). Insertion of definitions (2.3) and (2.4) produces the
compact expression

$T\dot{s}=\sum_{i=0}^{n}\vec{d}_{i}\cdot\vec{v}_{i}.$ (4.34)

Putting the Onsager–Stefan–Maxwell thermodiffusion laws (2.5)-(2.6) into
equation (4.34) produces equation (3.8), which shows that the dissipation
associated with temperature gradients when the species velocities vanish is
proportional to
$1/\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}$.
Past analyses have asserted on this basis that
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}$
should be inversely proportional to the thermal conductivity [20, 39]; in that
case force-explicit equation (2.5) reduces to Fourier’s law of heat conduction
in the absence of relative species velocities.

Alternatively, inserting Onsager–Fick–Fourier equations (4.30) into equation
(4.34) and applying Gibbs–Duhem relation (2.14) leads to

$T\dot{s}=\sum_{i,j=0}^{n}\vec{d}_{i}\cdot\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}^{\psi}_{ij}\vec{d}_{j},$
(4.35)

showing that the dissipation associated with temperature gradients in the
absence of diffusion driving forces is proportional to
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}_{00}$.
One could alternatively argue that thermal conductivity should be the
coefficient of $-\vec{\nabla}T$ in equation (4.32), so that it reduces to
Fourier’s law in the absence of diffusion driving forces.

Newman pointed out that these two perspectives on thermal conductivity are
mutually exclusive [11]. Indeed, after inserting the definitions of
$\widetilde{\textbf{m}}_{0}$ and
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{l}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{l}}}_{0}^{\hskip
1.0pt\psi}$ into equation (4.32) and applying an identity similar to the one
presented in equation (3.9), one finds that

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}_{00}^{\hskip
1.0pt\psi}=\frac{1}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}+\frac{R}{2\rho^{2}\hat{C}_{p}^{2}T}\sum_{i,j=1}^{n}\frac{c_{i}c_{j}\mathscr{A}_{ij}^{2}}{c_{\textrm{T}}\mathscr{D}_{ij}},$
(4.36)

so the thermal conductivity cannot be proportional to both
$1/\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}$
and
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}_{00}$
simultaneously. Note that we have included Newman–Soret diffusivities to
emphasize that this relationship is independent of the convective velocity,
addressing another concern raised by Newman [11].

The appropriate definition of thermal conductivity is decided by adopting a
standard approach to its experimental measurement. Thermal conductivity in
fluids is typically measured by inducing a temperature gradient across a one-
dimensional cell with closed ends, enforcing no-molar-flux boundary
conditions. When such an apparatus reaches a steady state, the differential
material balances demand that species velocities vanish throughout the cell,
as well as at its ends. For this experiment, equation (2.5) is appropriate,
and it reduces to Fourier’s law if

$k=\frac{\rho^{2}\hat{C}_{p}^{2}T}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}$
(4.37)

defines the thermal conductivity $k$. The dissipation functional from equation
(3.8) requires that $k\geq 0$.

One consequence of this definition of $k$ is that the generalized Fourier’s
law contains temperature-gradient-driven heat flux in addition to
$-k\vec{\nabla}T$. The extra contribution, which Bird, Stewart, and Lightfoot
call the _Dufour flux_ [16], can be quantified by inserting equations (4.36)
and (4.37) into equation (4.32). One practical problem the experimental
definition of $k$ raises is that thermodynamic consistency demands _a priori_
knowledge of the Soret diffusivities when applying the Onsager–Fick–Fourier
laws.

Our choice of thermal conductivity also has implications for thermodynamic
stability. Since the second law of thermodynamics requires that $k\geq 0$, the
condition that dissipation must be positive in equation (4.36) without
diffusion driving forces places a lower bound on
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}_{00}^{\hskip
1.0pt\psi}$, namely

$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{L}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{L}}_{00}^{\hskip
1.0pt\psi}\geq\frac{R}{2\rho^{2}\hat{C}_{p}^{2}T}\sum_{i,j=1}^{n}\frac{c_{i}c_{j}\mathscr{A}_{ij}^{2}}{c_{\textrm{T}}\mathscr{D}_{ij}}.$
(4.38)

Here inclusion of the Newman–Soret diffusivity emphasizes this condition’s
convection invariance.

##  5 Properties and balances for viscous fluids

Thermodiffusion is simulated by solving a system of differential balance
equations that governs the transient spatial distributions of temperature and
composition. As well as necessitating a parameterization of the transport
coefficients that make up
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{M}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{M}}}$
or
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{L}}}^{\psi}$,
closure of this system requires various equilibrium material properties.
Generally these parameters are not constants, but rather state functions
dependent on local temperature, pressure, and composition. A key advantage of
the irreversible thermodynamics methodology is that every equilibrium state
function ultimately derives from the system’s free energy.

###  5.1 Diffusion driving forces

Since molar Gibbs energy depends on temperature, pressure, and composition,
all of the material properties depend locally on these basis variables, as
well. Thus, elementary thermodynamic dependences imply that each diffusion
driving force in a viscous fluid depends on the natural basis variables $T$,
$p$, and $\left\\{c_{i}\right\\}_{i=1}^{n}$ as

$\vec{d}_{i}=-RTc_{\textrm{T}}\sum_{j=1}^{n}\chi_{ij}\vec{\nabla}\left(\frac{c_{j}}{c_{\textrm{T}}}\right)+c_{i}\left(\frac{\overline{m}_{i}}{\rho}-\overline{V}_{i}\right)\vec{\nabla}p,$
(5.1)

where $\overline{V}_{i}$ is the partial molar volume of species $i$ and
$\chi_{ij}$ quantifies how the change in chemical activity of species $i$
varies with the mole fraction of species $j$ (i.e., $c_{j}/c_{\textrm{T}}$) at
fixed $T$ and $p$. Notably, the temperature-gradient term in the original
definition of $\vec{d}_{i}$, equation (2.8), cancels out when chemical-
potential gradients are expanded over the natural basis. Terms involving
composition gradients in equation (5.1) express the forces that drive mass
diffusion; the remaining term drives pressure diffusion.

Entries within the $n\times n$ activity-gradient matrix $\bm{\chi}$ must
adhere to many constraints. When mixing is thermodynamically ideal,
$\bm{\chi}=\textbf{I}$, where I is the identity matrix. Therefore, one
conventionally lets $\bm{\chi}=\textbf{I}+\bm{\Gamma}$, such that the $n\times
n$ matrix $\bm{\Gamma}$ represents deviations from ideal mixing behaviour
under isothermal, isobaric conditions; Gibbs–Duhem relationships mandate that
$\textbf{1}^{\top}\bm{\Gamma}=\textbf{o}^{\top}$. A constraint that
$\bm{\Gamma}\textbf{1}=\textbf{o}$ ensures that nonphysical composition
changes lie in the nullspace of $\bm{\Gamma}$. Finally, although $\bm{\Gamma}$
is not generally symmetric, Maxwell relations still reduce its number of
independent entries to $n\left(n-1\right)/2$, one for each distinct pair of
species. B leverages results of a recent paper by Goyal and Monroe [40] to
produce a formula for $\bm{\Gamma}$ in terms of experimentally measurable
properties, and C shows how Maxwell relations constrain the entries of
$\bm{\chi}$.

###  5.2 Balance equations

Every multicomponent mass-transport process is governed by a set of species
material balances,222Generally, the material balances can be augmented by
generation terms, for example, to account for the formation or consumption of
species by homogeneous chemical reactions. Such phenomena, which produce
commensurate generation terms in the thermal energy balance as well, are
neglected here for simplicity.

$\frac{Dc_{i}}{Dt}=-\vec{\nabla}\cdot\left(c_{i}\vec{v}_{i}^{\hskip
1.0pt\rho}\right)-c_{i}\vec{\nabla}\cdot\vec{v}\hskip 1.0pt^{\rho}$ (5.2)

for $i\in\left\\{1,...,n\right\\}$. Here the operator

$\frac{D}{Dt}=\frac{\partial}{\partial t}+\vec{v}\hskip
1.0pt^{\rho}\cdot\vec{\nabla}$ (5.3)

represents the substantial derivative with respect to the mass-average
velocity.

Thermodynamic considerations reveal that the material balances imply a
connection between changes of temperature and pressure. A state equation for
the volume derives directly from the Euler equation for the Gibbs energy, by
differentiation with respect to pressure. The molar volume,
$1/c_{\textrm{T}}$, is an intensive equilibrium property and consequently
depends generally on $T$, $p$, and at most $n-1$ species fractions [19]. Thus
the local concentrations satisfy

$\sum_{i=1}^{n}\overline{V}_{i}c_{i}=1,$ (5.4)

in which the partial molar volumes depend on the same basis variables as
$c_{\textrm{T}}$. Since $\phi_{i}=\overline{V}_{i}c_{i}$ represents the volume
fraction occupied by species $i$ within a phase, this can be seen as a rather
trivial statement, that species volume fractions sum to unity. But it also
establishes the Gibbs phase rule, a more substantive physical conclusion.

Equation (5.4) leads to the Gibbs phase rule as follows. A substantial
derivative of equation (5.4) combines with the Gibbs–Duhem relation derived
from the extensivity of volume [19] to give

$\frac{Dp}{Dt}=K\alpha_{V}\frac{DT}{Dt}+K\sum_{i=1}^{n}\overline{V}_{i}\frac{Dc_{i}}{Dt},$
(5.5)

where $K$ is the bulk modulus (inverse isothermal compressibility) and
$\alpha_{V}$ is the volumetric coefficient of thermal expansion. Therefore, in
an $n$-ary single phase, one of the intensive natural variables $T$, $p$,
$\left\\{c_{i}\right\\}_{i=1}^{n}$ necessarily depends on the others.

To close the governing system’s degree of freedom associated with temperature,
a balance of thermal energy is added alongside the material balances. In any
nonreactive single-phase viscous fluid, this heat balance is [19]

$\displaystyle\rho\hat{C}_{p}\frac{DT}{Dt}=-\vec{\nabla}\cdot\vec{q}\hskip
1.5pt^{\prime}-\sum_{i=1}^{n}c_{i}\vec{v}_{i}^{\hskip
1.0pt\rho}\cdot\vec{\nabla}\overline{H}_{i}+T\alpha_{V}\frac{Dp}{Dt}$
$\displaystyle-\
\left(\vec{\vec{\tau}}+\sum_{i=1}^{n}c_{i}\overline{m}_{i}\vec{v}_{i}^{\hskip
1.0pt\rho}\vec{v}_{i}^{\hskip 1.0pt\rho}\right):\vec{\nabla}\vec{v}\hskip
1.0pt^{\rho},$ (5.6)

in which $\overline{H}_{i}$ is the partial molar enthalpy of species $i$ and
$\vec{\vec{\tau}}$ is deformation stress, the latter following Bird, Stewart,
and Lightfoot’s sign convention [16]. The terms on the right of equation (5.2)
respectively quantify local heat accumulation from: net influx of irreversible
heat; excess enthalpy convection due to diffusion; externally imposed
compressive power; and viscous friction.

###  5.3 Species enthalpies

The heat balance apparently introduces a new set of properties — the partial
molar enthalpies — but much of the information they contain is embedded in
quantities that have already been defined. In light of the standard Legendre
transformation $\mu_{i}=\overline{H}_{i}-T\overline{S}_{i}$, the definition of
the diffusion driving force and equation (5.1) give

$\vec{\nabla}\overline{H}_{i}=\overline{V}_{i}\vec{\nabla}p+\frac{RTc_{\textrm{T}}}{c_{i}}\sum_{j=1}^{n}\chi_{ij}\vec{\nabla}y_{j}+T\vec{\nabla}\overline{S}_{i}.$
(5.7)

Also, $\overline{S}_{i}$ connects to the chemical potential through a Maxwell
relation, so the last term depends on existing parameters. As discussed in B,
one ultimately finds

$\displaystyle\vec{\nabla}\overline{H}_{i}=\overline{C}_{p,i}\vec{\nabla}T+\left[1-\left(\frac{\partial\ln\overline{V}_{i}}{\partial\ln
T}\right)_{p}\right]\overline{V}_{i}\vec{\nabla}p-$
$\displaystyle\frac{RTc_{\textrm{T}}}{c_{i}}\sum_{j=1}^{n}\left(\frac{\partial\chi_{ij}}{\partial\ln
T}\right)_{p}\vec{\nabla}\left(\frac{c_{j}}{c_{\textrm{T}}}\right),$ (5.8)

in which $\overline{C}_{p,i}$ represents the partial molar constant-pressure
heat capacity of species $i$.

Substituting material balances (5.2) into equation (5.5) shows that the time
derivative of pressure in equation (5.2) depends only on spatial derivatives
and the time derivative of temperature. Insertion of equation (5.8) then
yields an explicit equation for the time evolution of temperature as a
function of the instantaneous spatial distributions of temperature, pressure,
composition, the set of excess species velocities $\left\\{\vec{v}_{i}^{\hskip
1.0pt\rho}\right\\}_{i=1}^{n}$, and the mass-average velocity $\vec{v}\hskip
1.0pt^{\rho}$. Since the Euler equation for entropy implies that

$\rho\hat{C}_{p}=\sum_{i=1}^{n}c_{i}\overline{C}_{p,i},$ (5.9)

this heat balance is fully parameterized by specifying the temperature,
pressure, and composition dependences of $\alpha_{V}$, $K$,
$\left\\{\overline{C}_{p,i}\right\\}_{i=1}^{n}$,
$\left\\{\overline{V}_{i}\right\\}_{i=1}^{n}$, the molar masses
$\left\\{\overline{m}_{i}\right\\}_{i=1}^{n}$, and the $n(n-1)/2$ independent
state functions that underpin $\bm{\chi}$. Bearing in mind that mass density
satisfies

$\rho=\sum_{i=1}^{n}\overline{m}_{i}c_{i}$ (5.10)

by definition, the parameters that have already been stated suffice to
quantify the diffusion driving forces, as well.

###  5.4 Dissipation

The only quantity that remains free in heat balance (5.2) is the deformation
stress. In general, the dissipation function for a nonreactive viscous fluid
takes the form [19]

$T\dot{s}=\sum_{i=0}^{n}\vec{d}_{i}\cdot\vec{v}_{i}-\left(\vec{\vec{\tau}}+\sum_{i=1}^{n}c_{i}\overline{m}_{i}\vec{v}_{i}^{\hskip
1.0pt\rho}\vec{v}_{i}^{\hskip 1.0pt\rho}\right):\vec{\nabla}\vec{v}\hskip
1.0pt^{\rho},$ (5.11)

which augments equation (4.34) to account for viscosity. Viscous dissipation
was neglected up to now to keep the discussion simpler: as a consequence of
Curie’s principle, which states that transport processes are decoupled if
their associated fluxes differ in tensor rank by one [41], the addition of
viscous dissipation to the functional does not affect the derivation of
Onsager–Stefan–Maxwell laws. Closing a model that includes $\vec{\vec{\tau}}$
requires a distinct constitutive law for stress, such as Newton’s law of
viscosity.

##  6 Thermodiffusion in monatomic-gas mixtures

Numerous measurements have been made to quantify thermodiffusion in binary
noble-gas mixtures. To put this body of experimental data in context, we first
lay out equilibrium properties and the thermal balance equation for
multicomponent gas mixtures. We subsequently discuss the physical assumptions
that underpin experimental data processing, and tabulate the measured
thermodiffusion properties for various noble-gas pairs.

###  6.1 Equilibrium properties

Here we restrict the development to ideal mixtures of monatomic gases, for
which a substantial amount of experimental data is available. The equilibrium
material properties take simpler forms in this case, greatly simplifying both
the diffusion driving forces and the energy balance.

As well as satisfying $\bm{\chi}=\textbf{I}$, an ideal gas is defined in part
by the condition that $\overline{V}_{i}=RT/p$ for all $i$. Thus equation (5.4)
yields

$c_{\textrm{T}}=\frac{p}{RT},$ (6.1)

a relationship which can be used to demonstrate that $K=p$ and
$\alpha_{V}=1/T$, both independent of composition. Also, this form of the
partial molar volumes reduces equation (5.8) to

$\vec{\nabla}\overline{H}_{i}=\overline{C}_{p,i}\vec{\nabla}T.$ (6.2)

Generally the ideal-gas state is also defined by a condition that each
species’ partial molar heat capacity $\overline{C}_{p,i}$ equals that of pure
$i$, and further, that each pure-species molar heat capacity depends on
temperature alone.

In an ideal gas, equation (5.1) for the force driving diffusion of species $i$
reduces to

$\vec{d}_{i}=-p\nabla
y_{i}+y_{i}\left(\frac{\overline{m}_{i}c_{\textrm{T}}}{\rho}-1\right)\vec{\nabla}p,$
(6.3)

where we have eliminated $c_{\textrm{T}}$ with state equation (6.1) and
defined the mole fraction $y_{i}=c_{i}/c_{\text{T}}$ to simplify notation.
Thus the driving force for pressure diffusion is proportional to the
difference between the molar mass of species $i$ and the number-average molar
mass of the mixture,
$\rho/c_{\textrm{T}}=\sum_{i=1}^{n}\overline{m}_{i}y_{i}$.

The equilibrium thermal properties simplify even more if all the species
comprising an ideal gas are assumed to be monatomic. In a monatomic ideal gas,
$\overline{C}_{p,i}=\frac{5}{2}R$ for each $i$. Through equation (5.9), the
volumetric heat capacity of a monatomic ideal-gas mixture is therefore
$\rho\hat{C}_{p}=\frac{5}{2}c_{\textrm{T}}R$, independent of composition.

###  6.2 Practical constraints

As it presently stands, the governing system is one vector equation short of
closure. Letting $d$ be the number of spatial dimensions, transport is
governed by pressure continuity equation (5.5), heat balance (5.2), $n$
material balances (5.2), $d(n+1)$ Onsager–Fourier–Fick flux-law components
(4.30), and $d^{2}$ viscous-fluid constitutive-law components. These govern
the transient spatial distributions of $T$, $p$,
$\left\\{c_{i}\right\\}_{i=1}^{n}$, $\vec{q}\hskip 1.5pt^{\prime}$,
$\\{\vec{v}_{i}^{\hskip 1.0pt\rho}\\}_{i=1}^{n}$, $\vec{v}\hskip
1.0pt^{\rho}$, and $\vec{\vec{\tau}}$. Thus the model comprises
$\left[2+n+d\left(1+n+d\right)\right]$ scalar equations in
$\left[2+n+d\left(2+n+d\right)\right]$ scalar unknowns, leaving $d$ degrees of
freedom.

In practice, experimentalists studying gases have simplified the model with an
isobaric assumption, taking $p$ to be spatially uniform. For isobaric gas
thermodiffusion, the condition that $\vec{\nabla}p=\vec{0}$ simplifies
diffusion driving forces (6.3) to

$\vec{d}_{i}=-p\vec{\nabla}y_{i}.$ (6.4)

The assumed spatial uniformity still allows a time dependence of pressure,
however, so the substantial derivatives in equations (5.5) and (5.2) reduce to
$dp/dt$, rather than vanishing. Thus, although the isobaric assumption
simplifies the governing equations, it does not fully resolve the closure
problem, because it leaves the balance of equations and unknowns unchanged.

Whether a system is isobaric or not, the last degrees of freedom can be closed
naturally by including a detailed local momentum balance, such as Cauchy’s
equation, to specify how the convective velocity evolves in response to the
pressure and deformation-stress distributions. This has not been done in
practice, however, in part because most thermodiffusion experiments probe
steady states, use closed containers, and assume one-dimensional geometries.
Thus no-molar-flux conditions at the ends of the container enforce a condition
that all species velocities vanish everywhere.

In the sequel, we will neglect viscous effects for consistency with the
thermodiffusion literature, and will assume an isobaric condition with
$\vec{v}\hskip 1.0pt^{\rho}=\vec{0}$ in place of a formal momentum balance.
This simplification is consistent with literature assumptions for one-
dimensional systems. For simulations of closed containers with more spatial
dimensions at steady state, it amounts to an additional stipulation that the
mass-average velocity is irrotational.

###  6.3 Thermodiffusion property measurements

Experiments to parameterize thermodiffusion typically employ binary gas
mixtures. The measurement is performed in an apparatus designed to be one-
dimensional spatially, which is subjected to a temperature difference across
its ends. Designating one end of the apparatus by $a$ and the other by $b$,
experimentalists record the steady-state value of the parameter

$\alpha_{12}=\frac{\ln\left(y_{2}^{b}/y_{1}^{b}\right)-\ln\left(y_{2}^{a}/y_{1}^{a}\right)}{\ln
T^{b}/T^{a}},$ (6.5)

known as the thermal diffusion factor [42]. Conventionally, index $1$ is
assigned to the heavier of the two species.

An expression relating $\alpha_{12}$ to transport properties derives naturally
from the Onsager–Stefan–Maxwell laws. In a steady state, the condition
$\vec{v}^{\rho}=0$ simplifies material balances (5.2) dramatically, to

$0=-\vec{\nabla}\cdot\left(c_{i}\vec{v}_{i}\right)$ (6.6)

for each species $i$. In a one-dimensional apparatus with closed ends, these
balances imply that $\vec{v}_{i}=\vec{0}$ for every species $i$. Thus the
force-explicit flux law (2.5) can be used to eliminate $\vec{v}_{0}$ from
equation (2.6), showing that

$\vec{d}_{i}=\frac{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{i0}}{\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{00}}\vec{d}_{0}$
(6.7)

for each species. Substitution of driving forces (2.4) and (6.4) into relation
(6.7), followed by simplification with equation of state (6.1) and
rearrangement based on the fact that $y_{2}=1-y_{1}$, give

$\vec{\nabla}\ln\left(\frac{y_{2}}{y_{1}}\right)=\frac{\mathscr{A}_{12}}{\mathscr{D}_{12}}\vec{\nabla}\ln
T$ (6.8)

for an isobaric binary mixture of monatomic gases in a closed, one-dimensional
apparatus subjected to a steady temperature polarization. Integration from end
$a$ to end $b$ yields

$\alpha_{12}=\frac{\mathscr{A}_{12}}{\mathscr{D}_{12}},$ (6.9)

assuming that the ratio $\mathscr{A}_{12}/\mathscr{D}_{12}$ is constant with
respect to temperature and composition. Thus the thermal diffusion factor
relates directly to the Newman–Soret coefficient and the Stefan–Maxwell
diffusivity.

Table 1 presents the binary diffusion coefficients and thermal diffusion
factors for various pairs of noble gases.

Gas pair |  $p\mathscr{D}_{12}$ ($\mathrm{a}\mathrm{t}\mathrm{m}\,\mathrm{c}\mathrm{m}^{2}\mathrm{s}^{-1}$) |  $\alpha_{12}$ (unitless)
---|---|---
Ne - He |  1.1079 |  0.3432
Ar - He |  0.7560 |  0.3984
Kr - He |  0.6553 |  0.4279
Xe - He |  0.5611 |  0.4250
Ar - Ne |  0.3250 |  0.1741
Kr - Ne |  0.2647 |  0.2710
Xe - Ne |  0.2231 |  0.2951
Kr - Ar |  0.1398 |  0.0705
Kr - Ar |  0.1144 |  0.0803
Xe - Kr |  0.0745 |  0.0262
Table 1: Thermodiffusion transport coefficients measured for noble-gas pairs
at $T=300~{}$\mathrm{K}$$, gathered from systems which comprise uniform
equimolar mixtures when relaxed to an isothermal equilibrium [43, 44].

In the first-order approximation of the Chapman–Enskog expansion, the terms
$p\mathscr{D}_{12}$ and $\alpha_{12}$ are independent of pressure [15]. Hence
the reported values are nominally valid for a range of pressures. Binary
diffusion coefficients were computed with the least-squares equations compiled
in the comprehensive review by Mason and Marrero [43]. The thermal separation
factors were taken from temperature-dependent least-squares correlations
reported by Taylor [44].

Many additional measurements of thermal gas-mixture separation exist in the
literature. The reported data cover hydrogen/noble-gas mixtures [45], mixtures
of hydrogen isotopes [46], halogen/noble-gas mixtures [47, 48], and gas
mixtures that include methane [49].

Complete specification of the thermodiffusion transport matrix requires a
value for the thermal conductivity. Thermal conductivities of pure noble gases
are reported in Table 2, taken from the survey of Kestin et al. [50].

Gas |  $k$ ($\mathrm{W}\,\mathrm{c}\mathrm{m}^{-1}\,\mathrm{K}^{-1}$)
---|---
He |  15.66
Ne |  4.98
Ar |  1.78
Kr |  0.95
Xe |  0.55
Table 2: Thermal conductivities $k$ for pure noble gases at
$300~{}$\mathrm{K}$$ [50].

Even simulating binary systems leads to complications because the thermal
conductivity is a strong function of composition; for example, the thermal
conductivities of xenon and helium differ by well over an order of magnitude.
Thus it is necessary to incorporate a mixing rule for thermal conductivities,
a topic for which there is a significant body of literature. Options include
the Wassijewa model [51] based on the kinetic theory of gases, the empirical
Kennard mixing rule [52], and the Mason–Saxena mixing rule [53]. Recently, an
accurate and easily implemented mixing rule based on equivalent-circuit
analysis was proposed by Udoetok [54].

##  7 Numerical computation

Practical modelling of thermodiffusion has been impeded by the difficulty of
numerical implementation. On one hand, deriving Onsager–Fick–Fourier laws for
multicomponent mixtures is a complex procedure, and, as we have illustrated,
the parameterization of the Onsager thermodiffusivity matrix
$\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{\textbf{L}}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{\textbf{L}}}^{\psi}$
depends strongly on the choice of reference velocity. On the other hand, in
the generalized Stefan–Maxwell equations (2.9), the mixture of fluxes and
driving forces apparently yields no useful mathematical structure to
facilitate numerical algorithms.

Recently the authors proposed an effective finite element discretization for
the isothermal Onsager–Stefan–Maxwell equations [30], formulated with an
arbitrary number of species. Simulations were designed to solve the force-
explicit transport equations for gases directly, coupling them to mass-
continuity equations for each species. Although that work used the mass-
average velocity as convective reference, it is straightforward to make
alternative choices. For further analysis of the finite-element solver, we
direct the reader to reference [30].

###  7.1 Extension of binary data

The Onsager–Stefan–Maxwell equations allow a natural extrapolation of
isothermal property data from binary systems to multicomponent systems.
Reasoning physically that species/species drag is dominated by pairwise
interactions to a first approximation, one can use diffusion coefficients
$\mathscr{D}_{ij}$ measured from binary systems directly in simulations with
more than two species. A similar process cannot be used to infer the set of
Newman–Soret diffusivities $\mathscr{A}_{ij}$ from binary measurements,
however, as a consequence of the structure required by equation (4.17).
Producing a consistent form of the multicomponent Newman–Soret matrix requires
that the binary data be regularized to ensure kinematic consistency. D
presents a regularization process underpinned by the method of least squares.

###  7.2 Steady thermodiffusion in a ternary noble-gas mixture

The Onsager–Stefan–Maxwell formalism casts thermodiffusion in such a way that
heat can be handled as another species within the transport constitutive laws,
allowing ready extension of the force-explicit solver developed previously
[30]. In the extended solver, flux laws (2.5) and (2.6) are coupled not only
to $n$ material balances in the form of equation (5.2), but also to a thermal
energy balance like equation (5.2).

Here we consider steady, three-dimensional thermodiffusion in an isobaric
ternary mixture of noble gases. In this situation, flux laws (2.5)-(2.6) for
species $i\in\left\\{0,1,2,3\right\\}$ are written in compact form as

$\vec{d}_{i}=\sum_{j=0}^{3}\mathchoice{\accentset{\displaystyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\textstyle\text{\smash{\raisebox{-6.45831pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptstyle\text{\smash{\raisebox{-4.52083pt}{$\widetildesym$}}}}{M}}{\accentset{\scriptscriptstyle\text{\smash{\raisebox{-3.22916pt}{$\widetildesym$}}}}{M}}_{ij}\vec{v}_{j},$
(7.1)

where it is understood that
$\vec{d}_{0}=-\frac{5}{2}c_{\textrm{T}}R\vec{\nabla}T$ and
$\vec{v}_{0}=\vec{q}\hskip 1.5pt^{\prime}/(\tfrac{5}{2}c_{\textrm{T}}RT)$. To
close the model we adopt the condition $\vec{v}^{\rho}=0$, which simplifies
the material balances to the form of equation (6.6) for each
$i\in\left\\{1,2,3\right\\}$. For a monatomic-gas mixture, the steady-state
heat balance simplifies to

$0=-\frac{5}{2}\vec{\nabla}\cdot(c_{\textrm{T}}RT\vec{v}_{0}),$ (7.2)

and equation (6.3) describes the diffusion driving forces.

###  7.3 Example

We illustrate the numerical method by simulating the separation of an
equimolar ternary gas mixture of helium, argon and krypton in a three
dimensional separation chamber. The conceptual schematic for this numerical
experiment is shown in Figure 1.

Figure 1: Schematic diagram of the separation chamber. The figure shows a
cross section of the axially symmetric device.

Values of binary Stefan–Maxwell coefficients $\mathscr{D}_{ij}$ were taken
from Table 1. Newman–Soret coefficients $\mathscr{A}_{ij}$ computed using
binary data for $\alpha_{12}$ are tabulated in Table 3, and the regularized
Newman–Soret matrix used for simulations is presented in Table 4. Thermal
conductivities for each pure component were taken from Table 2; to evaluate
the thermal conductivity of the mixture, the mixing rule proposed by Udoetok
[54] was employed. This mixing rule was chosen due to its good balance between
accuracy and simplicity, the latter of which is key because conductivity must
be computed dynamically from the local composition and temperature at each
point in the domain. The exact molar masses used for the simulation were
$4.00\>$\mathrm{g}\,\mathrm{m}\mathrm{o}\mathrm{l}^{-1}$$,
$39.95\>$\mathrm{g}\,\mathrm{m}\mathrm{o}\mathrm{l}^{-1}$$, and
$83.80\>$\mathrm{g}\,\mathrm{m}\mathrm{o}\mathrm{l}^{-1}$$ for helium, argon,
and krypton, respectively.

The simulation was performed with Firedrake software [55], using the MUMPS
direct linear solver [56, 57] via PETSc [58, 59]. The mesh of the geometry was
constructed using the Gmsh software [60]. Each linear system had 3,895,568
degrees of freedom and convergence was achieved in 11 non-linear iterations.
The code can be found at https://bitbucket.org/vanbrunt/consolidated-
thermodiffusion-repository along with the mesh used for the separation
chamber.

$\mathscr{A}_{ij}$ from binary measurements
---
|  He |  Ar |  Kr
He |  $0$ |  $-0.3012$ |  $-0.2804$
Ar |  $0.3012$ |  $0$ |  $-0.0099$
Kr |  $0.2804$ |  $0.0099$ |  $0$
Table 3: Table of Newman–Soret diffusivities taken directly from binary
measurements. Inferred $\mathscr{A}_{ij}$ entries
---
|  He |  Ar |  Kr
He |  $0$ |  $-0.2910$ |  $-0.2906$
Ar |  $0.2910$ |  $0$ |  $0.004$
Kr |  $0.2906$ |  $-0.004$ |  $0$
Table 4: Table of regularized Newman–Soret diffusivities for a He–Ar–Kr
mixture, determined from the data in Table 3 with the procedure detailed in D.

For the experiment, the cooling bath was maintained at 300 $\mathrm{K}$, a
condition reflected in the simulation by imposing Dirichlet boundary
conditions specifying the wall temperatures of the left chamber. To drive the
temperature gradient, an evenly-distributed fixed heat flux was fed normally
into the rightmost face of the device. The remaining boundary conditions were
set as homogeneous Neumann conditions to reflect the presence of insulation
and gas-impermeable walls. Last, it was assumed that the total molar amounts
of each gas are constants of the closed system. The constant total amount of
material in the separation chamber, $n_{\textrm{T}}$, necessarily satisfies

$\int_{V}c_{\textrm{T}}dV=n_{\textrm{T}},$ (7.3)

where $V$ indicates the chamber volume. Through equation of state (6.1), this
implies that

$p=\frac{Rn_{\textrm{T}}}{\int_{V}1/TdV},$ (7.4)

which was used to determine a self-consistent steady-state pressure from the
temperature distribution. For simulations, $n_{\textrm{T}}$ was chosen to be
consistent with an average total concentration given by equation (6.1) with
$p=1\,$\mathrm{a}\mathrm{t}\mathrm{m}$$ and $T=300\,$\mathrm{K}$$.

A heat input of $135\,$\mathrm{W}$$ induces a temperature difference of
approximately $100\,$\mathrm{K}$$ across the separation chamber, as can be
seen from Figure 2 alongside arrows that indicate the direction and magnitude
of local irreversible heat flux.

Figure 2: Temperature profile within the separation chamber. Arrows indicate
the irreversible-heat-flux ($\mathrm{W}\,\mathrm{c}\mathrm{m}^{-2}$) vector
field.

The separation of helium effected by the temperature gradient is plotted in
Figure 3. As expected, helium, the lightest component, tends to build up
towards the hotter end [61]; a $7.5\%$ separation of helium is induced by the
heat flux.

Figure 3: Mole fraction profile of helium.

The steady-state polarization of krypton induced by the heat flux is shown in
Figure 4. Krypton tends to move towards the cooler end, but the separation is
only about $4.5\%$.

Figure 4: Mole fraction profile of krypton.

This numerical example demonstrates the potential of the
Onsager–Stefan–Maxwell framework for modelling thermodiffusion in multispecies
phases with three-dimensional geometries. From this perspective, one can
investigate a rich variety of phenomena for a number of practical
applications. To the best of the authors’ knowledge, no other multidimensional
simulations of thermodiffusion have been performed using the linear
irreversible thermodynamics approach.

##  8 Conclusion

A unified framework for modelling fluid thermodiffusion was presented. The new
force-explicit Onsager–Stefan–Maxwell formulation was shown to be compatible
with the generalized Stefan–Maxwell equations implied by kinetic theory, as
well as the commonly employed flux-explicit Onsager–Fick–Fourier constitutive
laws. A general method for inverting isothermal transport laws to obtain
velocity-explicit forms from force-explicit forms, and vice versa, was
proposed, and shown to encompass the disparate prior approaches to flux-law
inversion proposed by Newman, Helfand, and Curtiss and Bird. The general
inversion process was further extended to cover non-isothermal situations,
proving equivalence of the various constitutive formulations.

The Onsager–Stefan–Maxwell thermodiffusion equations involve a non-isothermal
transport matrix that proves to have a simple spectral structure. It was shown
how to infer entries of multicomponent transport matrices given binary
measurements, allowing for straightforward practical implementation of the
theory. Further mathematical relationships necessary to implement a
thermodiffusion model, such as balance equations and expressions of the
driving forces in terms of natural thermodynamic basis variables, were
provided for thermodynamically non-ideal viscous fluids and ideal gases.
Finally, an effective and novel numerical method to simulate thermodiffusion
was deployed, in which heat was treated as a pseudospecies. Multidimensional
simulations of steady thermodiffusion in a ternary noble-gas mixture
demonstrated the effectiveness of the proposed approach.

Although the numerical methodology was shown to be sound, the absence of
convection (i.e., a zero mass-average velocity) is generally not a valid
physical assumption. Addressing this issue requires consideration of a local
momentum balance alongside the balances of material and heat. Future work is
needed to confront the full dynamics of convective thermodiffusion from the
irreversible-thermodynamics perspective.

## Acknowledgements

This work was supported by the Engineering and Physical Sciences Research
Council Centre for Doctoral Training in Partial Differential Equations:
Analysis and Applications (EP/L015811/1), the Engineering and Physical
Sciences Research Council (grants EP/R029423/1 and EP/V001493/1); the
Clarendon fund scholarship; and the Faraday Institution Multiscale Modelling
project (subaward FIRG003 under grant EP/P003532/1).

##  Appendix A Inversion of transport matrices

Both the procedure that maps M to $\textbf{L}^{\psi}$ for a given kinematic
relation $\bm{\psi}$, and the inverse process, by which $\textbf{L}^{\psi}$
maps to M given the Gibbs–Duhem relation 1, require some further
justification. Here we present formulas to construct $\textbf{L}^{\psi}$ from
M and M from $\textbf{L}^{\psi}$, as well as demonstrating that these two
mappings are inverses of each other.

Given a symmetric matrix R with a single null eigenvector r, and an additional
column matrix p such that $\textbf{p}^{\top}\textbf{r}\neq 0$, the augmented
matrix $\textbf{R}^{p}\left(\gamma\right)$, defined as

$\textbf{R}^{p}\left(\gamma\right)=\textbf{R}+\gamma\textbf{p}\textbf{p}^{\top},$
(A.1)

has rank one higher than R whenever $\gamma\neq 0$. This construction is
amenable to one of the generalized Sherman–Morrison formulas derived by
Baksalary et al. [62]. We employ the notation $\textbf{A}^{\sharp}$ to
represent the Moore–Penrose pseudoinverse of A. Further note that since there
is just one null eigenvector of R, the matrix Baksalary and colleagues call
$\textbf{Q}_{\textbf{R}^{\top}}$, which represents the orthogonal projector
onto the orthogonal complement of the column space of R, can be formed
explicitly:

$\textbf{Q}_{\textbf{R}^{\top}}=\frac{\textbf{r}\textbf{r}^{\top}}{\textbf{r}^{\top}\textbf{r}}.$
(A.2)

Equation (2.5) from reference [62] then yields

$\left[\textbf{R}^{p}\left(\gamma\right)\right]^{\sharp}=\left(\textbf{I}-\frac{\textbf{r}\textbf{p}^{\top}}{\textbf{r}^{\top}\textbf{p}}\right)\textbf{R}^{\sharp}\left(\textbf{I}-\frac{\textbf{p}\textbf{r}^{\top}}{\textbf{r}^{\top}\textbf{p}}\right)+\frac{1}{\gamma}\frac{\textbf{r}\textbf{r}^{\top}}{(\textbf{r}^{\top}\textbf{p})^{2}}$
(A.3)

after some algebraic rearrangement.

Next define the matrix P such that

$\textbf{P}=\left(\textbf{I}-\frac{\textbf{r}\textbf{p}^{\top}}{\textbf{r}^{\top}\textbf{p}}\right)\textbf{R}^{\sharp}\left(\textbf{I}-\frac{\textbf{p}\textbf{r}^{\top}}{\textbf{r}^{\top}\textbf{p}}\right).$
(A.4)

Direct calculation shows that P affords p as a null eigenvector, explaining
its notation. The matrix P so formed is unique because $\textbf{R}^{\sharp}$
is, and the fact that symmetry of R implies symmetry of $\textbf{R}^{\sharp}$
means $\textbf{P}=\textbf{P}^{\top}$ as well.

Since R has a single null eigenvector and $\textbf{R}^{p}\left(\gamma\right)$
has rank one higher than R, it follows that
$\textbf{R}^{p}\left(\gamma\right)$ is invertible. Therefore the pseudoinverse
$\left[\textbf{R}^{p}\left(\gamma\right)\right]^{\sharp}$ simplifies to the
standard matrix inverse; equation (A.3) rearranges to give

$\textbf{P}=\left[\textbf{R}^{p}\left(\gamma\right)\right]^{-1}-\frac{1}{\gamma}\frac{\textbf{r}\textbf{r}^{\top}}{(\textbf{r}^{\top}\textbf{p})^{2}}.$
(A.5)

This formula for P is independent of the choice of nonzero $\gamma$ and
circumvents the computation of $\textbf{R}^{\sharp}$. A trivial rearrangement
of this result,

$\left[\textbf{R}^{p}\left(\gamma\right)\right]^{-1}=\textbf{P}+\frac{1}{\gamma}\frac{\textbf{r}\textbf{r}^{\top}}{(\textbf{r}^{\top}\textbf{p})^{2}},$
(A.6)

substantiates the asymptotic behaviour asserted in equation (4.12). Direct
calculation then shows that

$\textbf{P}=\lim_{\gamma\to\infty}\left[\textbf{R}^{p}\left(\gamma\right)\right]^{-1},$
(A.7)

justifying the claim that the limit in equation (4.11) exists.

Finally we consider uniqueness of the mapping between P and R. Suppose we have
a new matrix, $\hat{\textbf{P}}$, whose sole null eigenvector is also the
previously defined column matrix p, and let r continue to be the sole null
eigenvector of R. For any $\phi\neq 0$ the augmented matrix

$\hat{\textbf{P}}^{r}\left(\phi\right)=\hat{\textbf{P}}+\phi\textbf{r}\textbf{r}^{\top}$
(A.8)

is invertible with rank one higher than $\hat{\textbf{P}}$. As before, the
formula of Baksalary et al. [62] constructs a unique matrix

$\hat{\textbf{R}}=\left[\hat{\textbf{P}}^{r}\left(\phi\right)\right]^{-1}-\frac{1}{\phi}\frac{\textbf{p}\textbf{p}^{\top}}{\left(\textbf{r}^{\top}\textbf{p}\right)^{2}}.$
(A.9)

Letting $\gamma=\phi^{-1}\left(\textbf{r}^{\top}\textbf{p}\right)^{-2}$ in
equation (A.1) and using that result to eliminate the second term on the right
produces an equality which proves that $\hat{\textbf{R}}=\textbf{R}$ if and
only if

$\left[\hat{\textbf{P}}^{r}\left(\phi\right)\right]^{-1}=\textbf{R}^{p}\left(\frac{1}{\phi\left(\textbf{r}^{\top}\textbf{p}\right)^{2}}\right).$
(A.10)

Since the matrices here are nonsingular, one can invert both sides, then
insert equations (A.6) and (A.8) to show that $\hat{\textbf{P}}=\textbf{P}$.
Thus our flux-law inversion establishes a bijective mapping between the space
of matrices $\textbf{P}=\textbf{P}^{\top}$ whose nullspaces are spanned by p
and the space of matrices $\textbf{R}=\textbf{R}^{\top}$ whose nullspaces are
spanned by r.

The transport-law inversions developed in section 4.2 can be implemented
directly using the general equations provided here. Specifically, to get flux-
explicit isothermal transport laws from the Onsager–Stefan–Maxwell
formulation, take $\textbf{r}=\textbf{1}$, $\textbf{p}=\bm{\psi}$, and
$\textbf{R}=\textbf{M}$, then use equation (A.4), (A.5), or (A.7) to produce
$\textbf{P}=\textbf{L}^{\psi}$. Alternatively, to retrieve force-explicit
isothermal transport laws from the isothermal Onsager–Fick–Fourier
formulation, let $\textbf{r}=\bm{\psi}$, $\textbf{p}=\textbf{1}$, and
$\textbf{R}=\textbf{L}^{\psi}$, in which case any of the same equations will
produce $\textbf{P}=\textbf{M}$.

##  Appendix B Non-ideal Darken factors

This appendix focuses on connecting measurable thermodynamic factors to the
properties $\chi_{ij}$ in equation (5.1). Goyal and Monroe analyzed nonideal
thermodynamic factors, emphasizing the practical constraint that only $n-1$
composition descriptors can be varied independently, and that species activity
coefficients are constrained by Gibbs–Duhem relations [40]. They developed
several results in terms of Darken factors [63] over the basis of mole
fractions $\left\\{y_{i}\right\\}_{i=1}^{n}$, where
$y_{i}=c_{i}/c_{\textrm{T}}$. This choice of basis makes clear the intrinsic
composition constraint, because

$\sum_{i=1}^{n}y_{i}=1$ (B.1)

everywhere, independent of $T$ and $p$. The Darken factor $\overline{Q}_{ij}$
is defined through the partial derivative

$\overline{Q}_{ij}=y_{i}\left(\frac{\partial\mu_{i}}{\partial
y_{j}}\right)_{T,p,y_{k\neq j,n}}.$ (B.2)

If $y_{j}$ varies, then mole-fraction sum (B.1) requires at least one other
mole fraction to vary. By requiring that the particular mole fraction $y_{n}$
be left free when taking any other composition derivative, Darken factors
defined in this way both ensure compatibility with the mole-fraction sum and
clarify what conditions to enforce during characterization experiments. Note
that $\overline{Q}_{ij}=RT\delta_{ij}$ if mixing is ideal.

This discussion requires matrices with various sizes and shapes, so it will
help to notate matrix dimensions. Here, a bold capital letter with a single
subscript index will imply a square matrix, with the index indicating the
number of rows; a bold lower-case letter with a single subscript index will
imply a column matrix, again with the index denoting total rows; row matrices
will always be annotated as transposes of columns; and rectangular matrices
with $m>1$ rows and $n>1$ columns will be annotated with capital letters and
subscript indices $m\times n$.

Our task here is to relate the nonideal activity-gradient matrix
$\bm{\Gamma}_{n}$ introduced in section 5 to the smaller Darken matrix
$\overline{\textbf{Q}}_{n-1}$. As we shall see shortly, this can be done
explicitly through matrix operations that exploit an invertible, symmetric
matrix $\textbf{G}_{n}$ with constant entries

$\left(\textbf{G}_{n}\right)_{ij}=\delta_{ij}-\delta_{in}-\delta_{nj}.$ (B.3)

One can write the inverse of $\textbf{G}_{n}$ as

$\left(\textbf{G}_{n}^{-1}\right)_{ij}=\delta_{ij}-\delta_{in}\delta_{nj}-\frac{1}{n}.$
(B.4)

The truncated matrix $\textbf{G}_{n\times(n-1)}$ formed by striking the
$n^{\textrm{th}}$ column of $\textbf{G}_{n}$ represents the linear
transformation process Goyal and Monroe call Gibbs–Duhem extension [40].

In terms of the Darken matrix $\overline{\textbf{Q}}_{n-1}$, the first $n-1$
diffusion driving forces are given by the matrix equation

$\left[\vec{d}_{i}\right]^{i=1}_{n-1}=-c_{\textrm{T}}\overline{\textbf{Q}}_{n-1}\left[\vec{\nabla}y_{i}\right]^{i=1}_{n-1}+\left[\frac{\rho_{i}}{\rho}-\phi_{i}\right]^{i=1}_{n-1}\vec{\nabla}p,$
(B.5)

an expression that only involves the first $n-1$ mole-fraction gradients.
Alternatively, all $n$ driving forces are written through equation (5.1) in a
form that includes all $n$ mole-fraction gradients. Multiplication by
$\textbf{G}_{n\times(n-1)}$ can extend the column of driving forces in
equation (B.5) as