# A multiscale thermodynamic generalization of Maxwell-Stefan diffusion
equations and of the dusty gas model

Petr Vágner<EMAIL_ADDRESS>Michal Pavelka  Jürgen Fuhrmann
Václav Klika  Weierstrass Institute, Mohrenstr. 39, 10117 Berlin, Germany
Mathematical Institute, Faculty of Mathematics and Physics, Charles University
in Prague, Sokolovská 83, 186 75 Prague, Czech Republic  Department of
Mathematics, FNSPE, Czech Technical University in Prague, Trojanova 13, Prague
2, 120 00, Czech Republic

###### Abstract

Despite the fact that the theory of mixtures has been part of non-equilibrium
thermodynamics and engineering for a long time, it is far from complete. While
it is well formulated and tested in the case of mechanical equilibrium (where
only diffusion-like processes take place), the question how to properly
describe homogeneous mixtures that flow with multiple independent velocities
that still possess some inertia (before mechanical equilibrium is reached) is
still open. Moreover, the mixtures can have several temperatures before they
relax to a common value. In this paper, we derive a theory of mixtures from
Hamiltonian mechanics in interaction with electromagnetic fields. The
resulting evolution equations are then reduced to the case with only one
momentum (classical irreversible thermodynamics), providing a generalization
of the Maxwell-Stefan diffusion equations. In a next step, we reduce that
description to the mechanical equilibrium (no momentum) and derive a non-
isothermal variant of the dusty gas model. These reduced equations are solved
numerically, and we illustrate the results on efficiency analysis, showing
where in a concentration cell efficiency is lost. Finally, the theory of
mixtures identifies the temperature difference between constituents as a
possible new source of the Soret coefficient. For the sake of clarity, we
restrict the presentation to the case of binary mixtures; the generalization
is straightforward.

###### keywords:

##  1 Introduction

Transport models are a hallmark application of non-equilibrium thermodynamics
(NET), backed up by a vast experimental evidence for the applicability of
linear closures between thermodynamic fluxes and forces [29]. But should one
consider a single or multiple constituents for a plausible description of the
transport? And how many independent velocities should be present in the
models? And what if temperatures of the constituents are not the same [8]?

Let us pause on the famous Duncan-Toor experiment where a rather unexpected
behaviour of diffusion was observed if the mixture is seen as a single
continuum. Various “regimes” of diffusion are reported which are not found in
the classical Fickian diffusion, where the diffusive flux occurs solely in the
direction of the negative chemical potential gradient, see [28, 24]. However,
if one just considers the same problem as a mixture, linear nonequilibrium
thermodynamics offers a resolution of these new phenomena by very standard
means, namely by coupling among the well-known driving forces of transport.
This model is known as Maxwell-Stefan model of diffusion [22]. Hence, there
are potentially key phenomena in transport that can be captured only on the
detailed mixture level. On one hand, Maxwell-Stefan models, as well as other
transport models, can be viewed as a special case of a more general NET
framework, the mixture theory 111We refer the interested reader to reviews [2,
4, 48, 23] and to classical texts [6, 9, 10, 11, 14, 17, 41, 44].. On the
other hand, there is a recognised distinguished limit of the Maxwell-Stefan
model, the dusty gas model of porous media transport.

We can order the various transport theories as follows:

  * 1.

no coupling of thermodynamic forces (no friction, i.e. transfer of momentum,
among the constituents), isothermal case – Fick’s law of transport

  * 2.

coupling among thermodynamic driving forces for friction, isothermal case –
Maxwell-Stefan model

  * 3.

coupling among thermodynamic driving forces, nonisothermal case – this
article.

Moreover, in this paper we go beyond by adding further state variables to the
classical densities and temperature, namely

  * 1.

electromagnetic fields,

  * 2.

one or two velocities 222we restrict ourselves to two component mixtures
without loss of generality

  * 3.

and possibly also two temperatures.

In order to derive such complex theories, however, we start with the most
detailed level, which can be derived from the Liouville equation [38], and
then we reduce it by letting the additional state variables relax.

In particular, we derive a nonisothermal generalisation of the Maxwell-Stefan
model that has the overall velocity as an independent state variable, so that
the evolution is not purely dissipative. Because we start from the Liouville
equation, we obtain also highly non-trivial terms in the evolution equation
for momentum density, which can not be found in previous works [5, 40, 42]. In
other words, by starting from the Liouville equation, we reveal mechanical
terms that are not visible when one starts from the standard balance equations
only. Moreover, the framework leads to a nonisothermal dusty gas model, which
we include as an application of the presented theory, as well as to a new
origin of the Soret coefficient [20].

The paper is organised as follows. Section 2 contains Hamiltonian and
thermodynamic evolution of binary homogeneous mixtures, derived from the
microscopic dynamics (Liouville equation), and Hamiltonian vacuum
electrodynamics. We rely on the General Equation for Non-Equilibrium
Reversible-Irreversible Coupling (GENERIC) [16, 33, 34, 38], which combines
Hamiltonian mechanics with gradient dynamics. Section 3 couples the motion of
mixtures with electromagnetic fields in a framework of Galilean invariant
electrodynamics of mixtures. Section 4 then contains a reduction of the
description to a model with a single velocity, which is thus in the realm of
Classical Irreversible Thermodynamics [10]. Section 5 further reduces the
model to the mechanical equilibrium, where no inertia is present anymore.
Section 6 contains a particular application of the reduced model – a
generalization of the dusty gas model. Section 7 shows an alternative and more
precise route towards the reduced models, where a detailed knowledge of the
entropy (either in terms of partial entropies or in terms of partial
temperatures) is not necessary. Finally, Section 8 shows the implication for
the modeling of the Soret coefficient.

##  2 GENERIC evolution of mixtures

Our goal in this section is to formulate models for a mixture at various
levels, with various levels of details, but in a systematic and
thermodynamically consistent way. To this end, we shall formulate our approach
within GENERIC framework [16, 33, 34, 38], which combines Hamiltonian
mechanics (generated by a Poisson bracket and energy) and gradient dynamics
(generated by a dissipation potential and entropy). Apart from its convenience
for multiscale problems, GENERIC clearly identifies and separates reversible
and irreversible processes in the evolution of state variables. We shall start
by the former.

###  2.1 Hamiltonian mechanics of mixtures

The reversible Hamiltonian evolution of the mixture is generated by a Poisson
bracket and energy (or free energy due to the assumption of isothermality), as
shown for example in [35]. The Poisson bracket is derived from the Liouville
equation [38] and reads

$\\{A,B\\}=\sum_{\alpha}\\{A,B\\}^{(FM)_{\alpha}}\leavevmode\nobreak\ ,$ (2.1)

where $A$ and $B$ are two arbitrary sufficiently regular functionals of the
state variables: partial density $\rho_{\alpha}$, partial momentum density
${u_{\alpha\,i}}$, and partial entropy density ${s_{\alpha}}$. The fluid
mechanics Poisson bracket for species $\alpha$, that is the bracket
corresponding to the fluid mechanics level for each individual species, is
given by

$\displaystyle\\{A,B\\}^{(FM)_{\alpha}}$ $\displaystyle=$
$\displaystyle\int\,\mathrm{d}\mathbf{r}\rho_{\alpha}\left(\partial_{i}A_{\rho_{\alpha}}B_{u_{\alpha\,i}}-\partial_{i}B_{\rho_{\alpha}}A_{u_{\alpha\,i}}\right)$
(2.2)
$\displaystyle+\int\,\mathrm{d}\mathbf{r}{u_{\alpha\,i}}\left(\partial_{j}A_{{u_{\alpha\,i}}}B_{{u_{\alpha\,j}}}-\partial_{i}B_{{u_{\alpha\,i}}}A_{{u_{\alpha\,j}}}\right)$
$\displaystyle+\int\,\mathrm{d}\mathbf{r}{s_{\alpha}}\left(\partial_{j}A_{s_{\alpha}}B_{u_{\alpha\,j}}-\partial_{i}B_{s_{\alpha}}A_{u_{\alpha\,j}}\right)\leavevmode\nobreak\
.$

Note that for instance expression $A_{\rho_{\alpha}}$ stands for functional
derivative of functional $A$ with respect to the field
${\rho_{\alpha}}(t,\mathbf{r})$, $\frac{\delta
A}{\delta{\rho_{\alpha}}(t,\mathbf{r})}$, see for instance [34, 38], and
$\partial_{i}$ stands for spatial derivative $\frac{\partial}{\partial
r^{i}}$. Einstein summation of spatial indexes is employed in the whole
manuscript.

Given total energy of the system $E$, reversible evolution of an arbitrary
functional $A$ can be then expressed in two ways, see [38],

$\displaystyle\left(\frac{\mathrm{d}A}{\mathrm{d}t}\right)_{rev}$
$\displaystyle=$ $\displaystyle\\{A,E\\}$ (2.3a) $\displaystyle=$
$\displaystyle\int\,\mathrm{d}\mathbf{r}\sum_{\alpha}A_{\rho_{\alpha}}\frac{\partial{\rho_{\alpha}}}{\partial
t}+\sum_{\alpha}A_{\mathbf{u}_{\alpha}}\cdot\frac{\partial{\mathbf{u}_{\alpha}}}{\partial
t}+\sum_{\alpha}A_{s_{\alpha}}\cdot\frac{\partial{s_{\alpha}}}{\partial t}.$
(2.3b)

By rewriting expression (2.3a) in the form of Eq. (2.3b), we can read the
evolution equations of ${\rho_{\alpha}}$, ${\mathbf{u}_{\alpha}}$ and
${s_{\alpha}}$,

$\displaystyle\left(\frac{\partial{\rho_{\alpha}}}{\partial t}\right)_{rev}$
$\displaystyle=$
$\displaystyle-\partial_{i}({\rho_{\alpha}}{u^{\dagger\,i}_{\alpha}})\leavevmode\nobreak\
,$ (2.4a) $\displaystyle\left(\frac{\partial{u_{\alpha\,i}}}{\partial
t}\right)_{rev}$ $\displaystyle=$
$\displaystyle-{\rho_{\alpha}}\partial_{i}{\rho^{\dagger}_{\alpha}}-{u_{\alpha\,j}}\partial_{i}{u^{\dagger\,j}_{\alpha}}-\partial_{j}\left({u_{\alpha\,i}}{u^{\dagger\,j}_{\alpha}}\right)\leavevmode\nobreak\
,$ (2.4b) $\displaystyle\left(\frac{\partial{s_{\alpha}}}{\partial
t}\right)_{rev}$ $\displaystyle=$
$\displaystyle-\partial_{i}({s_{\alpha}}{u^{\dagger\,i}_{\alpha}})$ (2.4c)
where the conjugate densities, ${\rho^{\dagger}_{\alpha}}$, and conjugate
momenta, ${\mathbf{u}^{\dagger}_{\alpha}}$, and conjugate entropies,
${s^{\dagger}_{\alpha}}$, are identified with the respective derivatives of
energy,
${\rho^{\dagger}_{\alpha}}=E_{\rho_{\alpha}},{\mathbf{u}^{\dagger}_{\alpha}}=E_{\mathbf{u}_{\alpha}}\mbox{
and }{s^{\dagger}_{\alpha}}=E_{s_{\alpha}}\leavevmode\nobreak\ .$ (2.4e)

###  2.2 Hamiltonian evolution of electromagnetic field in vacuum

The reversible evolution of the displacement field $\mathbf{D}$ and magnetic
field $\mathbf{B}$ in vacuum is generated by the Poisson bracket

$\displaystyle\\{A,B\\}^{(\textrm{EM})}=\int
A_{D^{i}}\varepsilon^{ijk}\partial_{j}B_{B^{k}}-B_{D^{i}}\varepsilon^{ijk}\partial_{j}A_{B^{k}}\leavevmode\nobreak\
,$ (2.5)

where we denoted with $\varepsilon^{ijk}$ the Levi-Civita symbol. For an
arbitrary energy E, the evolution equation of $\mathbf{D}$ and $\mathbf{B}$
read

$\displaystyle\partial_{t}D^{i}$
$\displaystyle=\varepsilon^{ijk}\partial_{j}E_{B_{k}}\leavevmode\nobreak\ ,$
(2.6a) $\displaystyle\partial_{t}B^{i}$
$\displaystyle=-\varepsilon^{ijk}\partial_{j}E_{D_{i}}\leavevmode\nobreak\ ,$
(2.6b) where $E_{\mathbf{D}}=\mathbf{E}$ is the electric field and
$E_{\mathbf{B}}=\mathbf{H}$ is the magnetic intensity, see [46], [38].

###  2.3 Irreversible evolution

Within the used GENERIC framework, irreversible evolution of the mixture is
generated by a dissipation potential [36]. Namely

$\displaystyle\left(\frac{\partial{\rho_{\alpha}}}{\partial t}\right)_{irr}$
$\displaystyle=$
$\displaystyle\frac{\delta\Xi}{\delta{\rho^{*}_{\alpha}}}\leavevmode\nobreak\
,$ (2.7a) $\displaystyle\left(\frac{\partial{\mathbf{u}_{\alpha}}}{\partial
t}\right)_{irr}$ $\displaystyle=$
$\displaystyle\frac{\delta\Xi}{\delta{\mathbf{u}^{*}_{\alpha}}}\leavevmode\nobreak\
,$ (2.7b) $\displaystyle\left(\frac{\partial{s_{\alpha}}}{\partial
t}\right)_{irr}$ $\displaystyle=$
$\displaystyle\frac{\delta\Xi}{\delta{s^{*}_{\alpha}}}\leavevmode\nobreak\ ,$
(2.7c) where the conjugates with respect to entropy appear
${\rho^{*}_{\alpha}}=S_{\rho_{\alpha}},{\mathbf{u}^{*}_{\alpha}}=S_{\mathbf{u}_{\alpha}}\mbox{
and }{s^{*}_{\alpha}}=S_{s_{\alpha}}.$ (2.7d) We assume the simplest possible
form of dissipation potential [38, Ch 6] being algebraic and quadratic in the
conjugate variables
$\Xi=\frac{1}{2}\int\,\mathrm{d}\mathbf{r}\frac{1}{\frac{1}{n}\sum_{\gamma}s_{\gamma}^{\dagger}}\frac{1}{2}\sum_{\stackrel{{\scriptstyle\alpha,\beta=1}}{{\beta\neq\alpha}}}^{n}\zeta_{\alpha\beta}({\mathbf{u}_{\alpha}}^{\dagger}-{\mathbf{u}_{\beta}}^{\dagger})^{2}+\frac{1}{\frac{1}{n}\sum_{\gamma}s_{\gamma}^{\dagger}}\frac{1}{2}\sum_{\stackrel{{\scriptstyle\alpha,\beta=1}}{{\beta\neq\alpha}}}^{n}K_{\alpha\beta}(s_{\alpha}^{\dagger}-s_{\beta}^{\dagger})^{2}$
(2.7e) where the prefactor
$\frac{1}{n}\sum_{\alpha=1}^{n}s_{\alpha}^{\dagger}$ corresponds to the
average temperature of the mixture ($s^{\dagger}$, the conjugate of the total
entropy $s=\sum_{\alpha}s_{a}$), as will be shown below.

Using Eq. (2.7e), the irreversible evolution becomes explicitly

$\displaystyle\left(\frac{\partial{\rho_{\alpha}}}{\partial t}\right)_{irr}$
$\displaystyle=0$ (2.8a)
$\displaystyle\left(\frac{\partial{\mathbf{u}_{\alpha}}}{\partial
t}\right)_{irr}$
$\displaystyle=\sum_{\beta}\zeta_{\alpha\beta}({\mathbf{u}^{\dagger}_{\beta}}-{\mathbf{u}^{\dagger}_{\alpha}})$
(2.8b) $\displaystyle\left(\frac{\partial{s_{\alpha}}}{\partial
t}\right)_{irr}$
$\displaystyle=\frac{1}{\sum_{\gamma=1}^{n}s_{\gamma}^{\dagger}}\left[\sum_{\beta,\gamma}\zeta_{\beta\gamma}({\mathbf{u}^{\dagger}_{\beta}}-\mathbf{u}_{\gamma}^{\dagger})^{2}+\sum_{\beta\neq\alpha}K_{\alpha\beta}(s_{\beta}^{\dagger}-s_{\alpha}^{\dagger})(s_{\beta}^{\dagger}+\sum_{\gamma\neq\alpha}s_{\gamma}^{\dagger})\right.$
$\displaystyle\phantom{sssssssssssss}\left.+\frac{1}{2}\sum_{\stackrel{{\scriptstyle\beta,\gamma=1}}{{\beta,\gamma\neq\alpha}}}^{n}K_{\beta\gamma}(s_{\beta}^{\dagger}-s_{\gamma}^{\dagger})^{2}\right]$
(2.8c)

The zero-th component was separated as it will be considered as the dust in
Sec. 6.

Coefficients $\zeta_{\alpha\beta}$ are to be interpreted as mutual friction
coefficients as in [3], and they are assumed to be symmetric, i.e.
$\zeta_{\alpha\beta}=\zeta_{\beta\alpha}$. This symmetry can be regarded as
the validity of Onsager reciprocal relations [31, 32, 37] which has been
recognised and accepted in particular in transport problems [25, 26].
Similarly, the heat transfer coefficients $K_{\alpha\beta}$ are also assumed
to be symmetric, i.e. $K_{\alpha\beta}=K_{\beta\alpha}$. In the following
section, we add coupling between mechanics of mixtures and electrodynamics.

##  3 Binary electromagnetohydrodynamics with friction

Since charged matter within electromagnetic fields is subject to the Lorentz
force, we have to combine the building blocks introduced in the previous
section to formulate the evolution equations of the binary
electromagnetohydrodynamic mixture with friction. First, we couple the fluid
mechanics Poisson bracket (2.1) with the electromagnetic Poisson bracket (2.5)
using the semi-direct product technique, see [46]. Second, we endow this
coupled Poisson bracket with the irreversible terms generated by the
dissipation potential (2.7e). We present the resulting evolution equations
below.

###  3.1 Evolution equations

Evolution equations describing a binary charged mixture read

$\displaystyle\partial_{t}\rho_{a}$
$\displaystyle=-\partial_{i}\left(\rho_{a}v_{a}^{i}\right)\leavevmode\nobreak\
,$ (3.1a) $\displaystyle\partial_{t}u_{a,i}$
$\displaystyle=-\rho_{a}\partial_{i}\frac{\partial\varepsilon}{\partial\rho_{a}}-s_{a}\partial_{i}\frac{\partial\varepsilon}{\partial
s_{a}}-\partial_{j}\left(u_{a,i}v_{a}^{j}\right)$
$\displaystyle+\rho_{a}\frac{q_{a}}{m_{a}}\left(E_{i}+\epsilon_{ijk}v_{a}^{j}B^{k}\right)+\zeta_{ab}(v_{b}^{i}-v_{a}^{i})\leavevmode\nobreak\
,$ (3.1b) $\displaystyle\partial_{t}s_{a}$
$\displaystyle=-\partial_{i}\left(s_{a}v_{a}^{i}\right)+\frac{1}{T_{a}+T_{b}}\left(\zeta_{ab}(\mathbf{u}_{b}^{\dagger}-\mathbf{u}_{a}^{\dagger})^{2}+2K_{ab}T_{b}(T_{b}-T_{a})\right)\leavevmode\nobreak\
,$ (3.1c) $\displaystyle\partial_{t}\rho_{b}$
$\displaystyle=-\partial_{i}\left(\rho_{b}v_{b}^{i}\right)\leavevmode\nobreak\
,$ (3.1d) $\displaystyle\partial_{t}u_{b,i}$
$\displaystyle=-\rho_{b}\partial_{i}\frac{\partial\varepsilon}{\partial\rho_{b}}-s_{b}\partial_{i}\frac{\partial\varepsilon}{\partial
s_{b}}-\partial_{j}\left(u_{b,i}v_{b}^{j}\right)\leavevmode\nobreak\ ,$
$\displaystyle+\rho_{b}\frac{q_{b}}{m_{b}}\left(E_{i}+\epsilon_{ijk}v_{b}^{j}B^{k}\right)+\zeta_{ab}(v_{a}^{i}-v_{b}^{i})\leavevmode\nobreak\
,$ (3.1e) $\displaystyle\partial_{t}s_{b}$
$\displaystyle=-\partial_{i}\left(s_{b}v_{b}^{i}\right)+\frac{1}{T_{a}+T_{b}}\left(\zeta_{ab}(\mathbf{u}_{b}^{\dagger}-\mathbf{u}_{a}^{\dagger})^{2}+2K_{ab}T_{a}(T_{a}-T_{b})\right)\leavevmode\nobreak\
,$ (3.1f) $\displaystyle\partial_{t}D^{i}$
$\displaystyle=\epsilon^{ijk}\partial_{j}H_{k}-\left(\frac{q_{a}\rho_{a}}{m_{a}}v_{a}^{i}+\frac{q_{b}\rho_{b}}{m_{b}}v_{b}^{i}\right)\leavevmode\nobreak\
,$ (3.1g) $\displaystyle\partial_{t}B^{i}$
$\displaystyle=-\epsilon^{ijk}\partial_{j}E_{k}\leavevmode\nobreak\ ,$ (3.1h)

where we denoted electric $\mathbf{E}$ and magnetic intensities $\mathbf{H}$,
and electric $\mathbf{D}$ and magnetic field $\mathbf{B}$, internal energy
$\varepsilon$ and charge flux $q_{\alpha}$ [38] for more details and relation
to the kinetic theory. By integration of the equations for $\mathbf{D}$ and
$\mathbf{B}$ with respect to time, the static Maxwell equations are obtained,

$\displaystyle\mathrm{div}\mathbf{D}=\frac{q_{a}\rho_{a}}{m_{a}}+\frac{q_{b}\rho_{b}}{m_{b}}\qquad\mbox{and}\qquad\mathrm{div}\mathbf{B}=0.$
(3.2)

Note that the total energy is conserved and the total entropy density,

$\displaystyle s=s_{a}+s_{b},$ (3.3)

grows. Volumetric total energy density can be prescribed, for example, as

$\displaystyle e$
$\displaystyle=\frac{1}{2}\rho_{a}\mathbf{v}_{a}^{2}+\frac{1}{2}\rho_{b}\mathbf{v}_{b}^{2}+\varepsilon(\rho_{a},\rho_{b},s_{a},s_{b})+\frac{1}{2\epsilon}\mathbf{D}^{2}+\frac{1}{2\mu}\mathbf{B}^{2}+\frac{\rho_{a}\mathbf{v}_{a}+\rho_{b}\mathbf{v}_{b}}{\rho_{a}+\rho_{b}}\cdot(\mathbf{D}\times\mathbf{B})\leavevmode\nobreak\
,$ (3.4)

so that it consists of the total kinetic energy density, internal energy and
we have $\mathbf{E}=\frac{\partial e}{\partial\mathbf{D}}$,
$\mathbf{H}=\frac{\partial e}{\partial\mathbf{B}}$, see [46]. In addition, we
adopt the standard notation for
$\mu_{a,b}=\frac{\partial\varepsilon}{\partial\rho_{a,b}}$,
$T_{a,b}=\frac{\partial\varepsilon}{\partial s_{a,b}}$, and
$\mathbf{m}=\frac{\partial e}{\partial\mathbf{v}_{a}}+\frac{\partial
e}{\partial\mathbf{v}_{b}}$ for the chemical potentials, partial temperatures,
and mass momentum, respectively.

The mixture is equipped with two entropies $s_{a}$ and $s_{b}$. That means
that the mixture has also two temperatures $T_{a}$ and $T_{b}$. This situation
is usual in low-temperature plasma, where electrons have higher temperature
than ions because they absorb energy of rapidly varying electromagnetic field
more efficiently [8].

###  3.2 Two-temperature local equilibrium thermodynamics

When summing the two equations for $\mathbf{u}_{a}$ and $\mathbf{u}_{b}$, one
can read the total pressure (noting that pressure follows from the definition
of the state variables and energy and hence is not subjected to any
constitutive choices),

$\displaystyle-p=-\epsilon+\rho_{a}\mu_{a}+\rho_{b}\mu_{b}+s_{a}T_{a}+s_{b}T_{b},$
(3.5)

and the Gibbs-Duhem relation

$\displaystyle-\mathrm{d}p=\rho_{a}\mathrm{d}\mu_{a}+\rho_{b}\mathrm{d}\mu_{b}+s_{a}\mathrm{d}T_{a}+s_{b}\mathrm{d}T_{b}$
(3.6)

or

$\displaystyle-\nabla
p=\rho_{a}\nabla\mu_{a}+\rho_{b}\nabla\mu_{b}+s_{a}\nabla T_{a}+s_{b}\nabla
T_{b}.$ (3.7)

Differential of the local equilibrium energy (standing for the internal energy
per a small volume) can then be written as

$\displaystyle\mathrm{d}E=T_{a}\mathrm{d}S_{a}+T_{b}\mathrm{d}S_{b}-p\mathrm{d}V+\mu_{a}\mathrm{d}M_{a}+\mu_{b}\mathrm{d}M_{b},$
(3.8)

the differential of the Gibbs free energy $G=E-T_{a}S_{a}-T_{b}S_{b}+pV$
becomes

$\displaystyle\mathrm{d}G=-S_{a}\mathrm{d}T_{a}-S_{b}\mathrm{d}T_{b}+V\mathrm{d}p+\mu_{a}\mathrm{d}M_{a}+\mu_{b}\mathrm{d}M_{b}.$
(3.9)

Therefore, we obtain the Maxwell relations

$\displaystyle s_{a}=\rho_{a}\left(\frac{\partial S_{a}}{\partial
M_{a}}\right)_{T_{a},T_{b},p,M_{b}}=-\rho_{a}\left(\frac{\partial\mu_{a}}{\partial
T_{a}}\right)_{T_{b},p,M_{a},M_{b}}.$ (3.10)

This definition of $s_{a}$ is indeed compatible with the foregoing definition,
i.e. the local volumetric entropy density. The Euler theorem for 1-homogeneous
functions implies for the total entropy, which is extensive, that

$\displaystyle S=S_{a}+S_{b}\qquad\mbox{for}\qquad
S_{a}=s_{a}\frac{M_{a}}{\rho_{a}},S_{b}=s_{b}\frac{M_{b}}{\rho_{b}}.$ (3.11)

Therefore, $s_{a}$ defined in (3.10) represent volumetric densities of
entropies of the species (when assuming no entropy attributed to the
microscopic interfaces between constituents).

Moreover, we can define the partial enthalpy by

$\displaystyle
h_{\alpha}=\mu_{\alpha}+T\frac{s_{\alpha}}{\rho_{\alpha}}=\mu_{\alpha}+T\left(\frac{\partial
S_{\alpha}}{\partial M_{\alpha}}\right)_{T,p,M_{\beta\neq\alpha}},$ (3.12)

which plays a role in the energy conservation below. In the next Section, we
reduce the rather complicated models from the current Section to less detailed
descriptions.

##  4 Binary electrohydrodynamics within Classical Irreversible Thermodynamics

Let us, in two steps, reduce model (3.1) to a less detailed model on the level
of the Classical Irreversible Thermodynamics (CIT) [10], that is to a
description with only one velocity. In the first step, the temperatures of the
constituents $T_{a}$ and $T_{b}$ will be equalized to a single temperature
$T_{a}=T=T_{b}$ using the MaxEnt principle. In the second step, the partial
momenta ${\mathbf{u}_{\alpha}}$ and ${\mathbf{u}_{\beta}}$ will be
approximated to the total momentum
$\mathbf{u}={\mathbf{u}_{\alpha}}+{\mathbf{u}_{\beta}}$ and the diffusion
velocity $\mathbf{w}^{\dagger}$ using the Dynamic MaxEnt. However, before
proceeding to the reductions, we let for simplicity the magnetic field relax
to $\mathbf{B}=0$.

###  4.1 MaxEnt reduction to one-temperature continuum

In the two-temperature treatment we were equipped with two entropies, $s_{a}$
and $s_{b}$. After the two temperatures have locally relaxed to their common
value $T$, we are equipped with only one entropy $s=s_{a}+s_{b}$. What is the
least-biased estimate of the partial entropies $s_{a}$ and $s_{b}$ in that
case? We shall use the MaxEnt principle [15, 38, 27].

Since we are working with entropies as state variables, it is more comfortable
to minimize energy rather than maximize entropy.333Both approaches are
equivalent as far as we have entropy or energy as a state variable [7]. They
are not equivalent for instance in kinetic theory, where neither energy nor
entropy are among the state variables and where only the distribution function
plays the role of state variable [38]. Consider a detailed energy
${{}^{\uparrow}}E(s_{a},s_{b},\rho_{a},\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b})$.
The reduction to the total entropy $s=s_{a}+s_{b}$ is carried out by the
reducing Legendre transformation, see [15],

$\displaystyle\frac{\partial}{\partial s_{a}}\left(-{{}^{\uparrow}}E+\langle
T,s_{a}+s_{b}\rangle\right)=0=\frac{\partial}{\partial
s_{b}}\left(-{{}^{\uparrow}}E+\langle T,s_{a}+s_{b}\rangle\right),$ (4.1)

which is equivalent to the equations

$\displaystyle\frac{\partial{{}^{\uparrow}}E}{\partial
s_{a}}=T=\frac{\partial{{}^{\uparrow}}E}{\partial s_{b}},$ (4.2)

solution of which are the dependencies $\tilde{s}_{a}(T,\rho_{a},\rho_{b})$
and $\tilde{s}_{b}(T,\rho_{a},\rho_{b})$, which we will further refer to as
_reduced partial entropy densities_. Note that these functions can be
explicitly obtained only after an explicit formula for ${{}^{\uparrow}}E$ has
been provided.

We can proceed in the Legendre transformation by defining the negative
Helmholtz free energy

$\displaystyle{{}^{\downarrow}}E^{*}(T,\rho_{a},\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b})=-{{}^{\uparrow}}E(\tilde{s}_{a}(T,\rho_{a},\rho_{b}),\tilde{s}_{b}(T,\rho_{a},\rho_{b}),\rho_{a},\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b})+T\cdot(\tilde{s}_{a}(T,\rho_{a},\rho_{b})+\tilde{s}_{a}(T,\rho_{a},\rho_{b})).$
(4.3)

Finally, the Legendre transformation (this time not reducing, but invertible)

$\displaystyle\frac{\partial}{\partial
T}\left(-{{}^{\downarrow}}E^{*}+Ts\right)=0,$ (4.4)

which leads to a dependency
$T(s,\rho_{a},\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b})$ and the final formula
for the reduced energy

$\displaystyle{{}^{\downarrow}}E(s,\rho_{a},\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b})=-{{}^{\downarrow}}E^{*}(T(s,\rho_{a},\rho_{b}),\rho_{a},\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b})+T(s,\rho_{a},\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b})s.$
(4.5)

Differential of the free energy reads

$\displaystyle
d{{}^{\downarrow}}E^{*}=-sdT+\mu_{a}d\rho_{a}+\mu_{b}d\rho_{b}+v_{a}^{i}du_{a}^{i}+v_{b}^{i}du_{b}^{i},$
(4.6)

which leads to the Maxwell relations

$\displaystyle\left(\frac{\partial
s}{\partial\rho_{a,b}}\right)_{T,\mathbf{u}_{a},\mathbf{u}_{b}}=-\left(\frac{\partial\mu_{a,b}}{\partial
T}\right)_{\rho_{a},\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b}}.$ (4.7)

Recalling that $s=s_{a}+s_{b}$, we can write, for instance, that

$\displaystyle\left(\frac{\partial
s_{a}}{\partial\rho_{a}}\right)_{T}+\left(\frac{\partial
s_{b}}{\partial\rho_{a}}\right)_{T}=-\left(\frac{\partial\mu_{a}}{\partial
T}\right)_{\rho_{a},\rho_{b}},$ (4.8)

for fixed partial momenta $\mathbf{u}_{a},\leavevmode\nobreak\
\mathbf{u}_{b}$. If the second term on the left hand side vanishes, as in the
case of non-interacting ideal gases [38], we obtain

$\displaystyle\left(\frac{\partial
s_{a}}{\partial\rho_{a}}\right)_{T}=-\left(\frac{\partial\mu_{a}}{\partial
T}\right)_{\rho_{a},\rho_{b}},$ (4.9)

which will become useful later on. Note however, that this formula ceases to
hold in case of interacting gases and one has to rely on (4.8).

###  4.2 DynMaxEnt reduction of the relative momentum $\mathbf{w}$

In order to reduce the dynamical equations (3.1), we first rewrite the two
partial momenta into the total momentum $\mathbf{u}$ and relative momentum
$\mathbf{w}$, as follows,

$\displaystyle\mathbf{u}$
$\displaystyle=\mathbf{u}_{a}+\mathbf{u}_{b}\leavevmode\nobreak\ ,$ (4.10a)
$\displaystyle\mathbf{w}$
$\displaystyle=\frac{1}{2}(\mathbf{u}_{a}-\mathbf{u}_{b})\leavevmode\nobreak\
,$ (4.10b) and conversely, $\displaystyle\mathbf{u}_{a}$
$\displaystyle=\frac{\mathbf{u}}{2}+\mathbf{w}\leavevmode\nobreak\ ,$ (4.10c)
$\displaystyle\mathbf{u}_{b}$ $\displaystyle=\frac{\mathbf{u}}{2}-\mathbf{w}.$
(4.10d)

The corresponding conjugate variables with respect to energy

$\displaystyle\mathbf{u}^{\dagger}$
$\displaystyle=\frac{1}{2}(\mathbf{u}^{\dagger}_{a}+\mathbf{u}^{\dagger}_{b})\leavevmode\nobreak\
,$ (4.11a) $\displaystyle\mathbf{w}^{\dagger}$
$\displaystyle=\mathbf{u}^{\dagger}_{a}-\mathbf{u}^{\dagger}_{b}\leavevmode\nobreak\
,$ (4.11b) and conversely, $\displaystyle\mathbf{u}^{\dagger}_{a}$
$\displaystyle=\mathbf{u}^{\dagger}+\frac{1}{2}\mathbf{w}^{\dagger}\leavevmode\nobreak\
,$ (4.11c) $\displaystyle\mathbf{u}^{\dagger}_{b}$
$\displaystyle=\mathbf{u}^{\dagger}-\frac{1}{2}\mathbf{w}^{\dagger}\leavevmode\nobreak\
.$ (4.11d) Note that the conjugate to the total momentum can be interpreted as
the barycentric velocity, $\mathbf{u}^{\dagger}=\mathbf{v}$.

In order to identify the reduced version of dynamics, we have to specify the
state variables on the lower (less-detailed) level together with its link to
the upper (more-detailed) level state variables. Consequently, entropy and
MaxEnt principle [27] provides the least biased estimate of the reduced
dynamics.

Our aim is to identify the dynamics of the system with a single macroscopic
velocity/momentum and temperature. Hence, the lower state variables are
$(\rho_{a},\rho_{b},\mathbf{u},s)$. The MaxEnt estimate of the relative
momentum $\mathbf{w}$ (while keeping the total momentum unchanged) reads

$\displaystyle\mathbf{u}_{a}=\frac{\rho_{a}}{\rho_{a}+\rho_{b}}\mathbf{u}\leavevmode\nobreak\
,\qquad$
$\displaystyle\mathbf{u}_{b}=\frac{\rho_{b}}{\rho_{a}+\rho_{b}}\mathbf{u}\leavevmode\nobreak\
,$ (4.12a) $\displaystyle\mathbf{u}=\mathbf{u}\leavevmode\nobreak\ ,\qquad$
$\displaystyle\mathbf{w}=\frac{\rho_{a}-\rho_{b}}{2(\rho_{a}+\rho_{b})}\mathbf{u}\leavevmode\nobreak\
,$ (4.12b)

see e.g. [38], which we shall substitute into the equations. The conjugate
variable $\mathbf{w}^{\dagger}$ is to be determined as the solution to the
stationary evolution equation for $\mathbf{w}$ [27]. Finally after the
reduction, Eqs. (3.1) become

$\displaystyle\partial_{t}\rho_{a}$
$\displaystyle=-\partial_{i}\left(\rho_{a}\left(v^{i}+\frac{w^{\dagger
i}}{2}\right)\right)\leavevmode\nobreak\ ,$ (4.13a)
$\displaystyle\partial_{t}\rho_{b}$
$\displaystyle=-\partial_{i}\left(\rho_{b}\left(v_{i}-\frac{w^{\dagger
i}}{2}\right)\right)\leavevmode\nobreak\ ,$ (4.13b)
$\displaystyle\partial_{t}u_{i}$
$\displaystyle=-\rho_{a}\partial_{i}\mu_{a}-\tilde{s}_{a}\partial_{i}T-\rho_{b}\partial_{i}\mu_{b}-\tilde{s}_{b}\partial_{i}T$
(4.13c)
$\displaystyle-\partial_{j}\left(\frac{\rho_{a}}{\rho_{a}+\rho_{b}}u_{i}\left(v^{j}+\frac{w^{\dagger
j}}{2}\right)+\frac{\rho_{b}}{\rho_{a}+\rho_{b}}u_{i}\left(v^{j}-\frac{w^{\dagger
j}}{2}\right)\right)+\left(\rho_{a}\frac{q_{a}}{m_{a}}+\rho_{b}\frac{q_{b}}{m_{b}}\right)E_{i}$
$\displaystyle=-\rho_{a}\partial_{i}\mu_{a}-\tilde{s}_{a}\partial_{i}T-\rho_{b}\partial_{i}\mu_{b}-\tilde{s}_{b}\partial_{i}T-\partial_{j}\left(u_{i}v^{j}+\frac{\rho_{a}-\rho_{b}}{2(\rho_{a}+\rho_{b})}u_{i}w^{\dagger
j}\right)+\left(\rho_{a}\frac{q_{a}}{m_{a}}+\rho_{a}\frac{q_{b}}{m_{b}}\right)E_{i}\leavevmode\nobreak\
,$ $\displaystyle\partial_{t}s$
$\displaystyle=-\partial_{i}\left(\tilde{s}_{a}\left(v^{i}+w^{\dagger
i}/2\right)+\tilde{s}_{b}\left(v^{i}-w^{\dagger
i}/2\right)\right)+\frac{1}{2T}\zeta_{ab}(\mathbf{v}_{a}-\mathbf{v}_{b})^{2}\leavevmode\nobreak\
,$ (4.13d) $\displaystyle\zeta_{ab}w^{\dagger i}$
$\displaystyle=\frac{1}{2}\left(-\rho_{a}\partial_{i}\mu_{a}-\tilde{s}_{a}\partial_{i}T+\frac{\rho_{a}}{m_{a}}q_{a}E_{i}+\rho_{b}\partial_{i}\mu_{b}+\tilde{s}_{b}\partial_{i}T-\frac{\rho_{b}}{m_{b}}q_{b}E_{i}\right)+\frac{1}{4}\partial_{j}\left(u_{i}w^{\dagger
j}\right)$ (4.13e) $\displaystyle\partial_{t}D^{i}$
$\displaystyle=-\left(\frac{q_{a}\rho_{a}}{m_{a}}v_{a}^{i}+\frac{q_{b}\rho_{b}}{m_{b}}v_{b}^{i}\right),$
(4.13f) where the total entropy $s=\tilde{s}_{a}+\tilde{s}_{b}$ and the
reduced partial entropies $\tilde{s}_{a,b}(T)$ are defined as functions of
temperature $T$; see Eqs. (4.2). We refer the reader to the section 5.2 where
we work out an example of the reduced entropies $\tilde{s}_{a,b}$ .

Neither the term in the equation for $\mathbf{u}$ that depends on
$\mathbf{w}^{\dagger}$, nor the terms in the equation for
$\mathbf{w}^{\dagger}$ that depend on $\mathbf{u}$ could be seen without
starting from the Liouville equation, which shows the advantages of our
approach to the standard approach based on balance equations [5].

####  4.2.1 Approximate diffusion velocity $\widehat{\mathbf{w}}^{\dagger}$

The velocity ${\mathbf{w}}^{\dagger}$ appears with alternating signs in the
continuity equations (4.13a) and (4.13b) besides the barycentric velocity
$\mathbf{v}=\mathbf{u}^{\dagger}$. It is thus tempting to interpret
$\mathbf{w}^{\dagger}$ as a diffusion velocity. How can one see that? To this
end, let us neglect the divergence term $\partial_{j}(u_{i}w^{\dagger j})$ in
equation (4.13e). We use the remaining terms to define the approximate
velocity,

$\displaystyle\widehat{w}^{\dagger
i}:=\frac{1}{2\zeta_{ab}}\left(-\rho_{a}\partial_{i}\mu_{a}-\tilde{s}_{a}\partial_{i}T+\frac{\rho_{a}}{m_{a}}q_{a}E_{i}+\rho_{b}\partial_{i}\mu_{b}+\tilde{s}_{b}\partial_{i}T-\frac{\rho_{b}}{m_{b}}q_{b}E_{i}\right)\leavevmode\nobreak\
.$ (4.14)

Clearly, the velocity $\widehat{\mathbf{w}}^{\dagger}$ is proportional to the
difference of the thermodynamic forces444gradients of the chemical potential
and temperature, and electric field acting upon the mixture species. Hence,
together with its role in the continuity equations, (4.13a) and (4.13b),
$\widehat{\mathbf{w}}^{\dagger}$ is the diffusion velocity. This results is
consistent with the canonical CIT formulation [10].

This observation that $\mathbf{w}^{\dagger}$ corresponds to diffusion velocity
is in line with the definition of the direct variable $\mathbf{w}$ itself. In
particular, $\mathbf{u}_{\alpha}^{\dagger}$ is the partial velocity as can be
seen from the standard choice of energy above, (3.4), and the conjugate
variables are related via
$\mathbf{u}_{\alpha}^{\dagger}=\mathbf{u}^{\dagger}\pm\frac{1}{2}\mathbf{w}^{\dagger}$
and hence, as $\mathbf{u}^{\dagger}$ is the barycentric velocity, we expect
the $\mathbf{w}^{\dagger}$ variable to represent the diffusion velocity.

However, in contrast to the usual formulations of CIT, such as [10], the
diffusive velocity appears in the evolution equation of the total momentum
(4.13c). What is its role therein? Let us integrate the momentum equation
(4.13c) over a volume $V$. We get

$\displaystyle\partial_{t}\int_{V}u_{i}$ $\displaystyle=\int_{\partial
V}-\left(-\varepsilon+\rho_{a}\mu_{a}+\rho_{b}\mu_{b}+(\widetilde{s}_{a}+\widetilde{s}_{b})T\right)\nu_{i}+\int_{V}\left(\rho_{a}\frac{q_{a}}{m_{a}}+\rho_{a}\frac{q_{b}}{m_{b}}\right)E_{i}$
$\displaystyle\phantom{--}-\int_{\partial
V}u_{i}\left(v^{j}+\frac{\rho_{a}-\rho_{b}}{2(\rho_{a}+\rho_{b})}\widehat{w}^{\dagger
j}\right)\nu_{j}\leavevmode\nobreak\ ,$ (4.15)

where $\nu_{i}$ denote components of the outer normal. The first integral on
the right hand side represents the material pressure; the second expresses the
force acting upon $V$ due to the electric field $\mathbf{E}$; and the last
integral is the transport of the total momentum $\mathbf{u}$ across the
boundary $\partial V$. In other words, the total momentum $\mathbf{u}$ is on
top of the (barycentric)-velocity field $\mathbf{v}$ also advected by the
density-weighted diffusion velocity
$\frac{\rho_{a}-\rho_{b}}{2(\rho_{a}+\rho_{b})}\widehat{\mathbf{w}}^{\dagger}$.
In the case that $\rho_{a}\gg\rho_{b}$, for example a dilute aqueous solution
of a salt, the $\widehat{\mathbf{w}}^{\dagger}$ could become in certain flow
regimes comparable to $\mathbf{v}$.

##  5 Reduction to mechanical equilibrium

In this Section, we further reduce the description to the level of mechanical
equilibrium, where no velocity plays the role of state variable and thus no
inertial effects are present. We start again with system (3.1) and use the
MaxEnt single-temperature reduction from Section 4.1. The partial momenta are
put equal to their MaxEnt values, $\mathbf{u}_{a}=0=\mathbf{u}_{b}$; and as it
follows from the Dynamic MaxEnt, the partial velocities $\mathbf{v}_{a}$ and
$\mathbf{v}_{b}$ are governed by the reduced momentum balances. As in the
previous section, we assume that magnetic field is weak enough, putting
$\mathbf{B}=0=\mathbf{H}$. The above-described reduction of the model (3.1)
reads

$\displaystyle\partial_{t}\rho_{a}$
$\displaystyle=-\partial_{i}\left(\rho_{a}v_{a}^{i}\right)\leavevmode\nobreak\
,$ (5.1a) $\displaystyle\partial_{t}\rho_{b}$
$\displaystyle=-\partial_{i}\left(\rho_{b}v_{b}^{i}\right)\leavevmode\nobreak\
,$ (5.1b) $\displaystyle 0$
$\displaystyle=-\rho_{a}\partial_{i}\mu_{a}-\tilde{s}_{a}\partial_{i}T+\rho_{a}\frac{q_{a}}{m_{a}}E_{i}+\zeta_{ab}(v_{b}^{i}-v_{a}^{i})$
(5.1c) $\displaystyle 0$
$\displaystyle=-\rho_{b}\partial_{i}\mu_{b}-\tilde{s}_{b}\partial_{i}T+\rho_{b}\frac{q_{b}}{m_{b}}E_{i}+\zeta_{ab}(v_{a}^{i}-v_{b}^{i})$
(5.1d) $\displaystyle\partial_{t}D^{i}$
$\displaystyle=-\left(\frac{q_{a}\rho_{a}}{m_{a}}v_{a}^{i}+\frac{q_{b}\rho_{b}}{m_{b}}v_{b}^{i}\right)$
(5.1e) $\displaystyle\partial_{t}s$
$\displaystyle=-\partial_{i}\left(\tilde{s}_{a}v_{a}^{i}+\tilde{s}_{b}v_{b}^{i}\right)+\frac{1}{2T}\zeta_{ab}(\mathbf{v}_{a}-\mathbf{v}_{b})^{2},$
(5.1f) where $\mathbf{D}=\epsilon\mathbf{E}$ and $\mathbf{E}=-\nabla\varphi$.
Again, $\tilde{s}_{a,b}(T)$ are given by Eqs. (4.2) and
$s=\tilde{s}_{a}+\tilde{s}_{b}$.

Note that the two flux-force relations (following from the linear momentum
balance) do yield a single equation for the unknown flux
$\rho_{\alpha}\mathbf{v}_{\alpha}$ while the second relation follows from the
condition for mixture mechanical equilibrium. Hence
$\rho_{\alpha}\mathbf{v}_{\alpha}=-\rho_{\beta}\mathbf{v}_{\beta}$ and the sum
of (5.1)c,d yields Gibbs-Duhem relation. In addition, Galilean invariance has
been sacrificed as momentum is no longer among the state variables. The
equations are thus valid only in the chosen lab frame, where the total
momentum vanishes, that is,
$\rho_{\alpha}\mathbf{v}_{\alpha}+\rho_{\beta}\mathbf{v}_{\beta}=0$.

In the case of non-interacting gases, we use formula (4.9) to rewrite the
gradient of chemical potential of species $a$ at constant temperature:

$\displaystyle\nabla(\mu_{a})_{T}$
$\displaystyle=\left(\nabla\mu_{a}+\frac{s_{a}}{\rho_{a}}\nabla
T_{a}\right)_{T_{a}=T_{b}=T}=\nabla\mu_{a}+\left(\frac{\partial S}{\partial
M_{a}}\right)_{T,p,M_{b}}\nabla T$
$\displaystyle=\nabla\mu_{a}-\left(\frac{\partial\mu_{a}}{\partial
T}\right)_{p,M_{a},M_{b}}\nabla T\leavevmode\nobreak\ ,$ (5.2)

to further simplify the Eqs. (5.1c) and (5.1d),.

Equations (5.1c) and (5.1d) can be interpreted as the constitutive relations
for the unknown velocities $\mathbf{v}_{a}$ and $\mathbf{v}_{b}$. Sum of the
two equations expresses the momentum balance, using the Gibbs-Duhem relation
(3.7),

$0=-\nabla
p+\left(\rho_{a}\frac{q_{a}}{m_{a}}+\rho_{b}\frac{q_{b}}{m_{b}}\right)\mathbf{E}.$
(5.3)

The approximation of momentum balance in the mechanical equilibrium is thus
implied. The total energy density is of the following form
$e(\rho_{\alpha},\mathbf{D},s)=\int\mathbf{D}^{2}/(2\epsilon)+\epsilon(T,\rho_{a},\rho_{b})$.
Entropy is clearly produced in Eqs. (5.1) and energy is conserved. Indeed, one
obtains by chain rule for the total energy density
$e(\rho_{\alpha},\mathbf{D},s)$ that

$\partial_{t}e=-\nabla\cdot\mathbf{J}_{e},\qquad\mbox{for}\qquad\mathbf{J}_{e}=\sum_{\alpha=a,b}\rho_{\alpha}h_{\alpha}\mathbf{v}_{\alpha}.$
(5.4)

The equations are thus compatible with the standard conservation laws of mass
and energy and with the balance of entropy.

The evolution equation for $\mathbf{D}$ can be integrated in time to obtain
the Poisson equation (3.2), or

$-\mathrm{div}\left(\epsilon\nabla\varphi\right)=\frac{q_{a}\rho_{a}}{m_{a}}+\frac{q_{b}\rho_{b}}{m_{b}},$
(5.5)

which can be solved instead of Eq. (5.1e). The advantage of having
$\mathbf{D}$ among the state variables is the explicit conservation of energy.
On the other hand, the Poisson equation is more suitable for actual
calculation because boundary conditions are typically given by setting the
electrostatic potential.

The friction coefficient $\boldsymbol{\zeta}_{ab}$ should disappear when
$\rho_{a}\rightarrow 0$ or $\rho_{b}\rightarrow 0$, since then one of the
components is missing. Therefore, it is often expressed as

$\displaystyle\boldsymbol{\zeta}_{ab}=\frac{k_{B}T}{\mathfrak{D}_{ab}}\frac{\frac{\rho_{a}}{m_{a}}\frac{\rho_{b}}{m_{b}}}{\frac{\rho_{a}}{m_{a}}+\frac{\rho_{b}}{m_{b}}},$
(5.6)

where $\mathfrak{D}_{ab}$ is the Maxwell-Stefan interdiffusion coefficient as
defined in [28]. The advantage is that the interdiffusion coefficients are
typically more or less independent of composition [28]. Equations (5.1c) and
(5.1d) can be also rewritten as

$\displaystyle\rho_{a}\frac{\nabla({\tilde{\mu}}_{\alpha})_{T}}{k_{B}T}=\frac{1}{\frac{\rho_{a}}{m_{a}}+\frac{\rho_{b}}{m_{b}}}\sum_{\beta\neq\alpha}\frac{\rho_{\alpha}}{m_{\alpha}}\frac{\rho_{\beta}}{m_{\beta}}\frac{1}{\mathfrak{D}_{\alpha\beta}}\left(\mathbf{v}_{\beta}-\mathbf{v}_{\alpha}\right)$
(5.7)

for $\alpha,\beta=a,b$. These equations are also valid for mixtures with more
than two components. Notice the presence of the electrochemical potentials

$\displaystyle{\tilde{\mu}}_{\alpha}=\mu_{\alpha}+\frac{z_{a}e_{0}}{m_{a}}\varphi,$
(5.8)

where $z_{a}=q_{a}/e_{0}$ is the charge number and $e_{0}$ the elementary
charge. Equations (5.7) are referred to as the Maxwell-Stefan relations.
Equations (5.1a), (5.1b), (5.7), (5.5) are referred to as the generalized
Poisson-Nernst-Planck equations, and they serve as the basic equations
governing electrochemical processes. Note that, in addition, the total
momentum vanishes, $0=\sum\rho_{\alpha}\mathbf{u}_{\alpha}$, in mechanical
equilibrium, which accompanies the above velocity relations.

###  5.1 Heat conduction

We have seen that energy can be transported by a motion of matter and by the
electromagnetic field. Another means of energy transport is the heat
conduction. One can again formulate a detailed evolution equation for heat
including the inertia of the heat flux [27], but we shall restrict ourselves
to the simplest possible description of heat flux, the Fourier law. The heat
flux, which represents a part of energy flux, is proportional to the negative
of the gradient of temperature so that heat flows to the lower temperatures,
in accordance with the second law of thermodynamics,

$\mathbf{J}_{q}=-\lambda\nabla T,$ (5.9)

where $\lambda$ is the coefficient of thermal conductivity. The balance of
energy then becomes

$\partial_{t}e=-\nabla\cdot(J_{e}+J_{q}).$ (5.10)

The corresponding balance of entropy is

$\displaystyle\partial_{t}s=-\nabla\cdot\left(\sum_{\alpha}s_{\alpha}\mathbf{v}_{\alpha}+\frac{\mathbf{J}_{q}}{T}\right)+\mathbf{J}_{q}\cdot\nabla
T^{-1}+\frac{1}{2T}\zeta_{ab}(\mathbf{v}_{a}-\mathbf{v}_{b})^{2},$ (5.11)

so that Eq. (5.10) is satisfied. Note that the part of entropy balance that is
not under the divergence is called entropy production,

$\sigma_{s}=\frac{\lambda}{T^{2}}\nabla T\cdot\nabla
T+\frac{1}{4T}\zeta_{ab}(\mathbf{v}_{a}-\mathbf{v}_{b})^{2},$ (5.12)

a clearly non-negative quantity.

The overall heat flux into the system is denoted by $Q$,

$Q=-\int_{\partial\Omega}\mathbf{J}_{q}\cdot\mathrm{d}S.$ (5.13)

###  5.2 Example: free energy model of a regular binary solution

Let us consider a free energy model of a regular solution: a binary mixture of
ideal gasses with two temperatures and a heat-of-solution term. The free
energy, expressed in the number densities
$n_{\alpha}=\frac{\rho_{\alpha}}{m_{\alpha}}$, reads

$\displaystyle f(n_{a},n_{b},T_{a},T_{b})$
$\displaystyle=\sum_{\alpha\in\\{a,b\\}}\underbrace{c_{\alpha}^{\textrm{V}}n_{\alpha}T_{\alpha}\left(1-\ln\left(c_{\alpha}^{\textrm{V}}n_{\alpha}^{1-\gamma_{\alpha}}T_{\alpha}\Phi_{\alpha}\right)\right)}_{\textrm{ideal
gas}}+\underbrace{(T_{a}\kappa_{a}+T_{b}\kappa_{b})\frac{n_{a}n_{b}}{n_{a}+n_{b}}}_{\textrm{generalized
regular solution}}\leavevmode\nobreak\ ,$ (5.14)

where $c_{\alpha}^{\textrm{V}}$ is the isochoric heat capacity;
$\gamma_{\alpha}$ is the adiabatic factor; $\Phi_{\alpha}$ is a constant that
follows from statistical physics underpinning; and finally $\kappa_{a}$ and
$\kappa_{b}$ are generalized heats of solution.

To illustrate the above theory, we identify the generalised Maxwell-Stefan
governing equations. To this end, we identify the $\tilde{s}_{a}$ as we assume
a detailed knowledge on the macroscale (5.14). Although the partial entropies
$s_{\alpha}$ have not been identified, the role of individual temperatures
$T_{\alpha}$ is equivalent as $s_{\alpha}=-\frac{\partial f}{\partial
T_{\alpha}}$.

The condition (4.2) defining the unknown functions $\tilde{s}_{\alpha}$ can be
simplified in the current case when the total free energy is a known function
of the partial temperatures as
$\tilde{s}_{\alpha}=s_{\alpha}|_{T_{\alpha}=T}$. Hence, we have

$\tilde{s}_{\alpha}=c_{\alpha}^{\textrm{V}}n_{\alpha}\ln\left(c_{\alpha}^{\textrm{V}}Tn_{\alpha}^{1-\gamma_{\alpha}}\Phi_{\alpha}\right)-\frac{n_{a}n_{b}}{n_{a}+n_{b}}\kappa_{\alpha}.$
(5.15)

Further, the chemical potential (per kilogram, as defined above) reads

$\displaystyle\mu_{\alpha}=$ $\displaystyle\frac{\partial
f}{\partial\rho_{\alpha}}\bigg{|}_{T_{a}=T_{b}=T}=\frac{\partial
f\big{|}_{T_{a}=T_{b}=T}}{\partial\rho_{\alpha}}=$
$\displaystyle=\frac{T}{m_{\alpha}}\left[c_{\alpha}^{\textrm{V}}\gamma_{\alpha}+\frac{\kappa_{a}+\kappa_{b}}{(n_{a}+n_{b})^{2}}n_{\beta}^{2}-c_{\alpha}^{\textrm{V}}\ln\left(c_{\alpha}^{\textrm{V}}Tn_{\alpha}^{1-\gamma_{\alpha}}\Phi_{\alpha}\right)\right]$

where $\alpha,\beta\in\\{a,b\\}$ and $\alpha\neq\beta$.

The governing equations read

$\displaystyle\partial_{t}\rho_{a}$
$\displaystyle=-\partial_{i}\left(\rho_{a}v_{a}^{i}\right)\leavevmode\nobreak\
,$ $\displaystyle\partial_{t}\rho_{b}$
$\displaystyle=\partial_{i}\left(\rho_{a}v_{a}^{i}\right)\leavevmode\nobreak\
,$ $\displaystyle v_{a}^{i}$
$\displaystyle=\frac{1}{\zeta_{ab}}\frac{\rho_{b}}{\rho}\left(-\rho_{a}\partial_{i}\mu_{a}-\tilde{s}_{a}\partial_{i}T\right)$
$\displaystyle=$
$\displaystyle\frac{1}{\zeta_{ab}}\frac{\rho_{b}}{\rho}\bigg{(}\frac{T}{(n_{a}+n_{b})^{3}}\left[c_{a}^{\textrm{V}}(1-\gamma_{a})n_{a}^{3}+3c_{a}^{\textrm{V}}(1-\gamma_{a})n_{a}^{2}n_{b}+(3c_{a}^{\textrm{V}}(1-\gamma_{a})+2(\kappa_{a}+\kappa_{b}))n_{a}n_{b}^{2}+c_{a}^{\textrm{V}}(1-\gamma_{a})n_{b}^{3}\right]\partial_{i}n_{a}$
$\displaystyle+\frac{2(\kappa_{a}+\kappa_{b})Tn_{a}^{2}n_{b}}{(n_{a}+n_{b})^{3}}\partial_{i}n_{b}$
$\displaystyle+\frac{n_{a}}{(n_{a}+n_{b})^{2}}\left[c_{a}^{\textrm{V}}(1-\gamma_{a})n_{a}^{2}+(2c_{a}^{\textrm{V}}(1-\gamma_{a})+\kappa_{a})n_{a}n_{b}+(c_{a}^{\textrm{V}}(1-\gamma_{a})-\kappa_{b})n_{b}^{2}\right]\partial_{i}T\bigg{)}$
$\displaystyle\partial_{t}s$
$\displaystyle=-\partial_{i}\left[\left(\tilde{s}_{a}-\frac{m_{a}}{m_{b}}\tilde{s}_{b}\right)v_{a}^{i}\right]+\frac{1}{4T}\zeta_{ab}\left(1-\frac{m_{a}}{m_{b}}\right)^{2}\mathbf{v}_{a}^{2}.$

To this end, we calculate the thermodynamic force of the component $\alpha$,
which contributes to the momentum equations (4.13c) and (4.13e), it reads

$\displaystyle-
n_{\alpha}\partial_{i}\mu_{\alpha}-\widetilde{s}_{\alpha}\partial_{i}T=$
$\displaystyle-(\gamma_{\alpha}-1)c_{\alpha}^{\textrm{V}}T\left(\partial_{i}n_{\alpha}+n_{\alpha}\partial_{i}\ln
T\right)-T(\kappa_{a}+\kappa_{b})\frac{2n_{a}n_{b}}{(n_{a}+n_{b})^{3}}\left(n_{\alpha}\partial_{i}n_{\beta}-n_{\beta}\partial_{i}n_{\alpha}\right)$
$\displaystyle+(\kappa_{a}+\kappa_{b})\frac{n_{a}n_{b}^{2}}{(n_{a}+n_{b})^{2}}\partial_{i}T-\kappa_{a}\frac{n_{a}n_{b}}{(n_{a}+n_{b})}\partial_{i}T\leavevmode\nobreak\
,$ (5.16)

where $\alpha,\beta\in\\{a,b\\}$ and $\alpha\neq\beta$. The coefficients
before the temperature gradient $\partial_{i}T$ in this case read

$\displaystyle-
n_{\alpha}(\gamma_{\alpha}-1)c_{\alpha}^{\textrm{V}}+\underbrace{\frac{n_{a}n_{b}}{(n_{a}+n_{b})^{2}}\left(n_{b}\kappa_{b}-n_{a}\kappa_{a}\right)}_{\textrm{non-
isothermal interaction}}$ (5.17)

In the following Section, we show another application, namely the dusty gas
model.

##  6 Dusty gas limit of the Maxwell-Stefan equations

Dusty gas is a widely applied model of porous medium, obtained as the limit of
Maxwell-Stefan transport equations [30], [28], [21]. To derive the dusty gas
model, let us consider the system (5.1) with an additional species, denoted
with subscript $D$. This additional species, representing the dust, is assumed
to be static and uniformly-distributed. In other words, we assume that
$n_{D}=const$ and $\mathbf{v}_{D}=0$. Note that the last assumption is the
main difference from the mechanical equilibrium case developed in Section 5.
We adjust the diffusion coefficients, similarly as in (5.6) as in [30], we
define
$\frac{\varepsilon}{\tau}nD_{\alpha\beta}=\mathfrak{D}_{\alpha\beta}(n+n_{D})$,
$\frac{\varepsilon}{\tau}D_{\alpha D}=\frac{n+n_{D}}{n_{D}}\mathfrak{D}_{aD}$,
where $\varepsilon$ is the dust’s volume fraction and $\tau$ stands for the
tortuosity. The force-flux relation yields

$\displaystyle
n_{\alpha}\partial_{i}\bar{\mu}_{\alpha}+s_{\alpha}\partial_{i}T+z_{\alpha}e_{0}n_{\alpha}\partial_{i}\varphi-{k_{B}T}\left(\sum_{\beta\neq
D}\frac{1}{\frac{\varepsilon}{\tau}nD_{\alpha\beta}}\left(n_{\alpha}J^{i}_{\beta}-n_{\beta}J^{i}_{\alpha}\right)-\frac{1}{\frac{\varepsilon}{\tau}D_{\alpha
D}}J^{i}_{\alpha}\right)$ $\displaystyle=0\leavevmode\nobreak\ ,\quad\mbox{for
}\alpha\neq D\leavevmode\nobreak\ ,$ (6.1)

where $n=\sum_{\alpha\neq D}n_{\alpha}$ and
$J_{\alpha}^{i}=n_{\alpha}v_{\alpha}^{i}$. The dusty gas limit of the Maxwell-
Stefan system reads

$\displaystyle\partial_{t}n_{\alpha}+\partial_{i}J_{\alpha}^{i}$
$\displaystyle=0\leavevmode\nobreak\ ,\quad\mbox{for }\alpha\neq
D\leavevmode\nobreak\ ,$ (6.2a)
$\displaystyle\partial_{i}\left(\varepsilon_{0}\varepsilon_{r}\partial_{i}\varphi\right)-e_{0}\sum_{\alpha}z_{\alpha}n_{\alpha}$
$\displaystyle=0\leavevmode\nobreak\ ,$ (6.2b)
$\displaystyle\partial_{t}\left(\epsilon(n_{\alpha},T)+\frac{\varepsilon_{0}\varepsilon_{r}}{2}\partial_{i}\varphi\partial_{i}\varphi\right)+\partial_{i}\left(-\lambda\partial_{i}T+\sum_{\alpha\neq
D}\left(\mu_{\alpha}+T\frac{\tilde{s}_{\alpha}}{n_{a}}\right)J_{\alpha}^{i}\right)$
$\displaystyle=0\leavevmode\nobreak\ .$ (6.2c)

###  6.1 Ideal gas and dust – Fick’s law

In particular for a mixture of neutral ideal gas and dust at constant
temperature, the dusty gas limit of the Maxwell-Stefan equations (6.1) can be
rewritten as

$\displaystyle-\frac{\varepsilon}{\tau}D_{aD}\nabla n_{a}=\mathbf{J}_{a},$
(6.3)

which is the Fick’s law, telling that flux of component $a$ is proportional to
negative of the gradient of concentration of the component.

###  6.2 Nonisothermal binary mixture of neutral ideal gasses in a porous
medium

Let us illustrate the system (6.2) on an example of a mixture of two neutral
ideal gasses and dust: $a$, $b$ and $D$. The free energy of the mixture reads

$\displaystyle\epsilon^{*}(n_{a},n_{b},T)$
$\displaystyle=c_{a}^{V}n_{a}T\left(1-\ln\left(c_{a}^{V}n_{a}^{1-\gamma_{a}}T\Phi_{a}\right)\right)+c_{b}^{V}n_{b}T\left(1-\ln\left(c_{b}^{V}n_{b}^{1-\gamma_{b}}T\Phi_{b}\right)\right)+c_{D}^{V}n_{D}T\leavevmode\nobreak\
,$ (6.4)

where $c_{\alpha}^{V}$ is the isochoric heat capacity, $\gamma_{\alpha}$ is
the adiabatic factor, and $\Phi_{\alpha}=e^{\gamma_{\alpha}}\frac{4\pi
m_{\alpha}}{3\hbar^{2}}$ is a constant that follows from statistical physics
underpinning, cf. Sackur-Tetrode entropy. The reduced momentum equations read

$\displaystyle\frac{c_{a}^{V}}{k_{B}}(\gamma_{a}-1)\left(\partial_{i}n_{a}+n_{a}\partial_{i}\ln
T\right)$
$\displaystyle=\frac{n_{a}}{n_{a}+n_{b}}\frac{1}{\frac{\varepsilon}{\tau}D_{ab}}J^{i}_{b}-\left(\frac{n_{b}}{n_{a}+n_{b}}\frac{1}{\frac{\varepsilon}{\tau}D_{ab}}+\frac{1}{\frac{\varepsilon}{\tau}D_{aD}}\right)J^{i}_{a}\leavevmode\nobreak\
,$ (6.5a)
$\displaystyle\frac{c_{b}^{V}}{k_{B}}(\gamma_{b}-1)\left(\partial_{i}n_{b}+n_{b}\partial_{i}\ln
T\right)$
$\displaystyle=\frac{n_{b}}{n_{a}+n_{b}}\frac{1}{\frac{\varepsilon}{\tau}D_{ab}}J^{i}_{a}-\left(\frac{1}{\frac{\varepsilon}{\tau}D_{bD}}+\frac{n_{a}}{n_{a}+n_{b}}\frac{1}{\frac{\varepsilon}{\tau}D_{ab}}\right)J^{i}_{b}\leavevmode\nobreak\
.$ (6.5b)

A stationary numerical solution of the system (6.2) with Dirichlet boundary
condition is illustrated in 1. The model formulation and the finite volume
discretization are summarized in the B.

####  6.2.1 Analysis of irreversibility of the nonisothermal binary mixture of
neutral ideal gasses in a porous medium

The particular form of the entropy production density follows from the general
formula (6.6d) and it reads

$\displaystyle\sigma_{\text{Fourier}}=$ $\displaystyle\
\lambda\left(\frac{\partial_{i}T}{T}\right)^{2}$ (6.6a)
$\displaystyle\sigma_{\text{MS}}=$ $\displaystyle\
\frac{n_{a}n_{b}}{n_{a}+n_{b}}\frac{k_{B}}{\frac{\varepsilon}{\tau}D_{ab}}\left(\frac{J_{a}^{i}}{n_{a}}-\frac{J_{b}^{i}}{n_{b}}\right)^{2}$
(6.6b) $\displaystyle\sigma_{\text{dust}}=$ $\displaystyle\
{n_{a}}\frac{k_{B}}{\frac{\varepsilon}{\tau}D_{aD}}\left(\frac{J_{a}^{i}}{n_{a}}\right)^{2}+{n_{b}}\frac{k_{B}}{\frac{\varepsilon}{\tau}D_{bD}}\left(\frac{J_{b}^{i}}{n_{b}}\right)^{2}$
(6.6c) $\displaystyle\sigma=$ $\displaystyle\
\sigma_{\text{Fourier}}+\sigma_{\text{MS}}+\sigma_{\text{dust}}$ (6.6d)

We can now plot a quantity assessing the locations of lost power. Once we know
what is a useful power in a given problem, typically determined by what can be
controlled in a given system, we can define the map of losses
$\textrm{MOL}(x)$–the density of lost power555energy dissipation rate. Let us
consider that we can control the mass fluxes across the boundary but not
temperature and temperature gradients. Then, the map of losses can be
retrieved from the stationary energy balance (6.2c) as follows

$\displaystyle 0=$
$\displaystyle\int_{\Omega}\partial_{i}\left(-\lambda\partial_{i}T+T\frac{\tilde{s}_{a}}{n_{a}}J^{i}_{a}+T\frac{\tilde{s}_{b}}{n_{b}}J^{i}_{b}\right)+\partial_{i}\left(\mu_{a}J^{i}_{a}+\mu_{b}J^{i}_{b}\right)$
$\displaystyle=$
$\displaystyle\int_{\Omega}\partial_{i}\left(T\left[-\lambda\frac{\partial_{i}T}{T}+\frac{\tilde{s}_{a}}{n_{a}}J^{i}_{a}+\frac{\tilde{s}_{b}}{n_{b}}J^{i}_{b}\right]\right)+\partial_{i}\left(\mu_{a}J^{i}_{a}+\mu_{b}J^{i}_{b}\right)$
$\displaystyle=$
$\displaystyle\underbrace{\int_{\partial\Omega}\left(\mu_{a}J^{i}_{a}+\mu_{b}J^{i}_{b}\right)\cdot\nu_{i}}_{\dot{\Delta}G}+\int_{\Omega}\underbrace{\left[\partial_{i}T\left(-\lambda\frac{\partial_{i}T}{T}+\frac{\tilde{s}_{a}}{n_{a}}J^{i}_{a}+\frac{\tilde{s}_{b}}{n_{b}}J^{i}_{b}\right)+T\sigma\right]}_{=:MOL(x)}\leavevmode\nobreak\
.$ (6.7)

The spatial profile of $\textrm{MOL}(x)$ is shown in Fig. 1, the entropy
production density overestimates the power losses which supports the findings
in [39] and [47].

Figure 1: Map of losses (MOL), see (6.7) shows the density of energy losses
due to the binary diffusion and heat conduction in a porous medium (B.1).

##  7 Full reduction to CIT

The reductions described in the preceding two sections are only approximate,
as the reduction is not complete, the reduced equations are not closed with
respect to the lower state variables (e.g. $s_{\alpha}$ still persist). One
requires to provide some information from the detailed level, detailed entropy
in our case, to be able to calculate $\tilde{s}_{a},\leavevmode\nobreak\
\tilde{s}_{b}$.

Here, we take a different route where the reduction is complete and the
reduced equations are closed with respect to the less detailed state
variables. The price we pay is the complexity of the obtained equations. To
this end, we consider a reduction in direct and conjugate state variables. As
we shall see, the reduction in conjugate variables is more natural and general
because it does not require the particular knowledge of the detailed entropy.
In this sense, the approach presented above and the projection to CIT in
direct variables is similar as they both require the specification of the
detailed entropy to proceed.

###  7.1 Projection to CIT in direct variables

Starting point is the binary mixture description without charge [38], i.e.
Eqns (3.1) where $\mathbf{E}=\mathbf{B}=0$. With an explicit choice of entropy
(entropy of ideal gases) we may use the dynamic MaxEnt [38, 15] to identify
the evolution equations on the CIT level in the mechanical equilibrium, i.e.
state variables $\rho_{a},\leavevmode\nobreak\ \rho_{b},\leavevmode\nobreak\
\mathbf{u}=\mathbf{u}_{a}+\mathbf{u}_{b},\leavevmode\nobreak\ e=e_{a}+e_{b}$
where momentum $\mathbf{u}$ is constant.

The dynamics on the microscale corresponding to binary fluid mechanics
$\rho_{a},\leavevmode\nobreak\ \rho_{b},\leavevmode\nobreak\
\mathbf{u}_{a},\leavevmode\nobreak\ \mathbf{u}_{b},\leavevmode\nobreak\
s_{a},\leavevmode\nobreak\ s_{b}$ is given by

$\frac{\partial\rho_{\alpha}}{\partial
t}=-\mathrm{div}(\rho_{\alpha}\mathbf{u}_{\alpha}^{\dagger}),$ (7.1a)
$\frac{\partial u_{\alpha,i}}{\partial
t}=-\rho_{\alpha}\partial_{i}\rho_{\alpha}^{\dagger}-u_{\alpha,j}\partial_{i}u_{\alpha,j}^{\dagger}-s_{\alpha}\partial_{i}s_{\alpha}^{\dagger}-\partial_{j}(u_{\alpha,i}u_{\alpha,j}^{\dagger})-(\delta_{a\alpha}-\delta_{b\alpha})\zeta_{ab}(u_{a,i}^{\dagger}-u_{b,i}^{\dagger})$
(7.1b) $\frac{\partial s_{\alpha}}{\partial
t}=-\mathrm{div}(s_{\alpha}\mathbf{u}_{\alpha}^{\dagger})+\frac{1}{2}\frac{s_{\alpha}^{\dagger}}{(s_{a}^{\dagger}+s_{b}^{\dagger})^{2}}\sum_{\beta,\gamma}\zeta_{\gamma\beta}(\mathbf{u}_{\gamma}^{\dagger}-\mathbf{u}_{\beta}^{\dagger})^{2}-\frac{2}{s_{a}^{\dagger}+s_{b}^{\dagger}}K_{ab}(s_{a}^{\dagger}-s_{b}^{\dagger})(\delta_{b\alpha}s_{a}^{\dagger}-\delta_{a\alpha}s_{b}^{\dagger}),$
(7.1c)

with a typical energy

$E=\int\mathrm{d}\mathbf{r}\left(\frac{\mathbf{u}_{a}^{2}}{2\rho_{a}}+\frac{\mathbf{u}_{b}^{2}}{2\rho_{b}}+\epsilon(\rho_{a},s_{a},\rho_{b},s_{b})\right).$
(7.2)

The effect of temperature equilibration and mutual friction is described via a
dissipation potential as above, eq. (2.7e) while, to complete the description,
we consider the entropy potential for a binary mixture of ideal gases

$S(\rho_{a},\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b},e_{a},e_{b})=S_{a}+S_{b}=\sum_{\alpha\in\\{a,b\\}}\int\mathrm{d}\mathbf{r}k_{B}\frac{\rho_{\alpha}}{m_{\alpha}}\left(\frac{5}{2}+\ln\left[\frac{m_{\alpha}}{\rho_{\alpha}}\left(\frac{4\pi
m_{\alpha}}{3h^{2}}\frac{e_{\alpha}-\frac{\mathbf{u}_{\alpha}^{2}}{2\rho_{\alpha}}}{\rho_{\alpha}/m_{\alpha}}\right)^{3/2}\right]\right).$
(7.3)

Finally, the classical irreversible thermodynamics (CIT) level is
characterised by the following choice of state variables
$\rho_{a},\leavevmode\nobreak\ \rho_{b},\leavevmode\nobreak\
\mathbf{u},\leavevmode\nobreak\ s=s_{a}+s_{b}$. It is again advantageous to
change the state variables in binary hydrodynamics prior to this projection to
CIT level so that one can easily track the effects of different temperatures
or dissipation related to the relative motion. To this end we again introduce

$\mathbf{u}=\mathbf{u}_{a}+\mathbf{u}_{b},\leavevmode\nobreak\
\mathbf{w}=1/2(\mathbf{u}_{a}-\mathbf{u}_{b})$

and define

$s=s_{a}+s_{b},\leavevmode\nobreak\ \sigma=1/2(s_{a}-s_{b}).$

Similarly as for $u^{\dagger},\leavevmode\nobreak\ w^{\dagger}$ discussed
above, one can show [38] that the conjugate variables are related via

$\mathbf{u}^{\dagger}=1/2(\mathbf{u}_{a}^{\dagger}+\mathbf{u}_{b}^{\dagger}),\leavevmode\nobreak\
\mathbf{w}^{\dagger}=\mathbf{u}_{a}^{\dagger}-\mathbf{u}_{b}^{\dagger},\leavevmode\nobreak\
s^{\dagger}=1/2(s_{a}^{\dagger}+s_{b}^{\dagger}),\leavevmode\nobreak\
\sigma^{\dagger}=s_{a}^{\dagger}-s_{b}^{\dagger}.$

With the above explicit choice of entropy, the MaxEnt (MaxEnt values of state
variables are denoted with ME in the upper index) gives

$\frac{\mathbf{u}_{\alpha}^{ME}}{\rho_{\alpha}}=\frac{\mathbf{u}}{\rho_{T}},\quad\rho_{T}=\rho_{a}+\rho_{b},$

hence the most probable way of distributing the total momentum among partial
momenta is via mass fractions. Similarly MaxEnt values of partial energies
satisfy

$\frac{e_{\alpha}^{ME}-(\mathbf{u}_{\alpha}^{ME})^{2}/(2\rho_{\alpha})}{\text{\sout{$\rho$}}_{\alpha}}=\frac{e-\mathbf{u}^{2}/(2\rho_{T})}{\text{\sout{$\rho$}}_{T}},$

where $\text{\sout{$\rho$}}_{\alpha}=\rho_{\alpha}/m_{\alpha}$.

The explicit knowledge of entropy can be used to identify $\mathbf{u}=0$ in
the mechanical equilibrium and similarly
$\mathbf{w}^{ME}=\frac{1}{2}(\mathbf{u}_{a}-\mathbf{u}_{b})=\frac{\mathbf{u}}{2}\frac{\rho_{a}-\rho_{b}}{\rho_{a}+\rho_{b}}=0$
while

$\sigma^{ME}(\rho_{a},\rho_{b},s)=\frac{1}{2}(s_{a}-s_{b})=-\frac{5}{2}k_{B}\frac{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}}{2}+2k_{B}\frac{\text{\sout{$\rho$}}_{a}\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}}\ln\left(\left(\frac{m_{a}}{m_{b}}\right)^{3/2}\frac{\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{a}}\right)+\frac{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}}s,$
(7.4)

where we introduced
$\text{\sout{$\rho$}}_{T}=\text{\sout{$\rho$}}_{a}+\text{\sout{$\rho$}}_{b}$.

Dynamic MaxEnt provides a relation for direct state variables of the relaxed
state variables from the higher, more detailed, level of description (in our
case $\mathbf{w},\sigma$). The values of their conjugate variables have to be
identified so that the dynamics of the relaxed evolution on the lower level
stays exactly on the MaxEnt values in the direct variables.

Hence $\mathbf{w}_{\alpha}^{\dagger}$ and $\sigma^{\dagger}$ have to be
identified from their evolution equations following from their Max Ent values.
In particular,

$\partial_{t}\mathbf{u}_{a}^{ME}=\boxed{0=-\rho_{a}\nabla\rho_{a}^{\dagger}-\zeta_{ab}(\mathbf{u}_{a}^{\dagger}-\mathbf{u}_{b}^{\dagger})}-\underbrace{(s/2+\sigma^{ME})}_{s_{a}}\nabla\underbrace{(s^{\dagger}+\sigma^{\dagger}/2)}_{s_{a}^{\dagger}},$
(7.5a)
$\partial_{t}\mathbf{u}_{b}^{ME}=\boxed{0=-\rho_{b}\nabla\rho_{b}^{\dagger}-\zeta_{ab}(\mathbf{u}_{b}^{\dagger}-\mathbf{u}_{a}^{\dagger})}-\underbrace{(s/2-\sigma^{ME})}_{s_{b}}\nabla\underbrace{(s^{\dagger}-\sigma^{\dagger}/2)}_{s_{b}^{\dagger}},$
(7.5b)

where the boxed terms correspond to the classical Maxwell-Stefan model. The
MaxEnt value of $\sigma$, $\sigma^{ME}$, can be used to get

$\displaystyle\partial_{t}\sigma^{ME}$
$\displaystyle=-\frac{1}{2}\mathrm{div}\left[(\mathbf{u}_{a}^{\dagger}+\mathbf{u}_{b}^{\dagger})\left(-\frac{5}{2}k_{B}\frac{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}}{2}+2k_{B}\frac{\text{\sout{$\rho$}}_{a}\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}}\ln\left[\left(\frac{m_{a}}{m_{b}}\right)^{3/2}\frac{\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{a}}\right]\right)\right]$
(7.6a)
$\displaystyle-\frac{1}{4}\mathrm{div}\left[\frac{s}{\text{\sout{$\rho$}}_{T}}\left(\mathbf{u}_{a}^{\dagger}(3\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b})-\mathbf{u}_{b}^{\dagger}(3\text{\sout{$\rho$}}_{b}-\text{\sout{$\rho$}}_{a})\right)\right]-K_{ab}\sigma^{\dagger}$
$\displaystyle=\partial_{t}\text{\sout{$\rho$}}_{a}\left[-\frac{5}{4}k_{B}+2k_{B}\frac{\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}}\ln\left[\left(\frac{m_{a}}{m_{b}}\right)^{3/2}\frac{\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{a}}\right]-2k_{B}\frac{\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}}+2\frac{\text{\sout{$\rho$}}_{a}}{\text{\sout{$\rho$}}_{T}}s\right]$
(7.6b)
$\displaystyle+\partial_{t}\text{\sout{$\rho$}}_{b}\left[\frac{5}{4}k_{B}+2k_{B}\frac{\text{\sout{$\rho$}}_{a}}{\text{\sout{$\rho$}}_{T}}\ln\left[\left(\frac{m_{a}}{m_{b}}\right)^{3/2}\frac{\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{a}}\right]+2k_{B}\frac{\text{\sout{$\rho$}}_{a}}{\text{\sout{$\rho$}}_{T}}-2\frac{\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}}s\right]+\partial_{t}s\frac{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}},$

where the first equality follows from the general evolution equations for
$\sigma$, i.e.
$\partial_{t}\sigma=\frac{1}{2}(\partial_{t}s_{a}-\partial_{t}s_{b})$, with
$s_{a}=\frac{s}{2}+\sigma^{ME}$, $s_{b}=\frac{s}{2}-\sigma^{ME}$, while the
second equality comes from the time derivative of the MaxEnt estimate of
$\sigma$, $\sigma^{ME}$ in Eq. (7.4).

These three relations, Eqs (7.5), (7.6), for the conjugate values of the
relaxed state variables, $\mathbf{u}_{a}^{\dagger},\leavevmode\nobreak\
\mathbf{u}_{b}^{\dagger},\leavevmode\nobreak\ \sigma^{\dagger}$, provide
closures for the evolution of the lower state variables

$\displaystyle\partial_{t}\rho_{\alpha}$
$\displaystyle=-\mathrm{div}(\rho_{\alpha}\mathbf{u}_{\alpha}^{\dagger}),$
(7.7) $\displaystyle\partial_{t}s$
$\displaystyle=-\mathrm{div}\left\\{(\mathbf{u}_{a}^{\dagger}-\mathbf{u}_{b}^{\dagger})\left(-\frac{5}{2}k_{B}\frac{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}}{2}+2k_{B}\frac{\text{\sout{$\rho$}}_{a}\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}}\ln\left[\left(\frac{m_{a}}{m_{b}}\right)^{3/2}\frac{\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{a}}\right]\right)+\frac{1}{2}\frac{s}{\text{\sout{$\rho$}}_{T}}\left(\mathbf{u}_{a}^{\dagger}(3\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b})+\mathbf{u}_{b}^{\dagger}(3\text{\sout{$\rho$}}_{b}-\text{\sout{$\rho$}}_{a})\right)\right\\}+$
$\displaystyle+\frac{1}{s^{\dagger}}\left[\frac{1}{4}\zeta_{ab}(\mathbf{u}_{a}^{\dagger}-\mathbf{u}_{b}^{\dagger})^{2}+K_{ab}\left(\sigma^{\dagger}\right)^{2}\right].$
(7.8)

The standard Maxwell-Stefan model is obtained (the difference in the conjugate
momenta provide a relation for the chemical potentials
$\rho_{\alpha}^{\dagger}$) when the entropy is allowed to equilibrate (when
the temperature, $s^{\dagger}$, is such that $\partial_{t}s=0$). In the
extended version with the temperature effect, one has to solve for
$\sigma^{\dagger}$ which decouples from the
$\mathbf{u}_{a}^{\dagger},\leavevmode\nobreak\ \mathbf{u}_{b}^{\dagger}$
dependency once the two relations (7.5) are added.

####  7.1.1 Fast dissipation

In the case of a fast dissipation (the time scale required for the entropy and
momentum dissipation is significantly shorter than the time scale of
reversible evolution), we may assume that the terms $K_{ab}$ and $\zeta_{ab}$
are significantly larger than the remaining ones (as in Chapman-Enskog
solution of Boltzmann equation). Then the leading order equation for
$\sigma^{\dagger}$ reads from (7.6) as

$0=K_{ab}\sigma^{\dagger}+\frac{1}{s^{\dagger}}K_{ab}\left(\sigma^{\dagger}\right)^{2}\frac{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}},$

with two leading order solutions

$\displaystyle\sigma^{\dagger}=0,$ (7.9)
$\displaystyle\sigma^{\dagger}=-s^{\dagger}\frac{\text{\sout{$\rho$}}_{T}}{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}}.$
(7.10)

The second solution, however, is not a physically admissible as one of the
corresponding partial temperatures would have to be negative

$s_{a}^{\dagger}=s^{\dagger}+\sigma^{\dagger}/2=-s^{\dagger}\frac{2\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}},\quad
s_{b}^{\dagger}=s^{\dagger}-\sigma^{\dagger}/2=s^{\dagger}\frac{2\text{\sout{$\rho$}}_{a}}{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}}.$

However, the first leading order solution is always plausible and natural as
it corresponds to the same partial temperatures in CIT mixture where one
measures just a single temperature without the possibility to distinguish
momenta of each constituent.

In this setting the extended version of Maxwell-Stefan model simplifies to

$\displaystyle\boxed{\partial_{t}\rho_{\alpha}=-\mathrm{div}(\rho_{\alpha}\mathbf{u}_{\alpha}^{\dagger})},$
(7.11)
$\displaystyle\boxed{0=-\rho_{a}\nabla\rho_{a}^{\dagger}-\zeta_{ab}(\mathbf{u}_{a}^{\dagger}-\mathbf{u}_{b}^{\dagger})}-(s/2+\sigma^{ME}(\rho_{a},\rho_{b},s))\nabla
s^{\dagger},$ (7.12)
$\displaystyle\boxed{0=-\rho_{b}\nabla\rho_{b}^{\dagger}-\zeta_{ab}(\mathbf{u}_{b}^{\dagger}-\mathbf{u}_{a}^{\dagger})}-(s/2-\sigma^{ME}(\rho_{a},\rho_{b},s))\nabla
s^{\dagger},$ (7.13)
$\displaystyle\partial_{t}s=-\mathrm{div}\left\\{(\mathbf{u}_{a}^{\dagger}-\mathbf{u}_{b}^{\dagger})\left(-\frac{5}{2}k_{B}\frac{\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b}}{2}+2k_{B}\frac{\text{\sout{$\rho$}}_{a}\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{T}}\ln\left[\left(\frac{m_{a}}{m_{b}}\right)^{3/2}\frac{\text{\sout{$\rho$}}_{b}}{\text{\sout{$\rho$}}_{a}}\right]\right)+\frac{1}{2}\frac{s}{\text{\sout{$\rho$}}_{T}}\left(\mathbf{u}_{a}^{\dagger}(3\text{\sout{$\rho$}}_{a}-\text{\sout{$\rho$}}_{b})+\mathbf{u}_{b}^{\dagger}(3\text{\sout{$\rho$}}_{b}-\text{\sout{$\rho$}}_{a})\right)\right\\}+$
(7.14)
$\displaystyle+\frac{1}{s^{\dagger}}\zeta_{ab}(\mathbf{u}_{a}^{\dagger}-\mathbf{u}_{b}^{\dagger})^{2}.$
(7.15)

Note again that, in addition, the total momentum vanishes,
$0=\sum\rho_{\alpha}\mathbf{u}_{\alpha}$, in mechanical equilibrium, which
closes the problem.

Finally note that $s^{\dagger}=T$ and using the assumption of the explicit
form of entropy one can rewrite both $s$ and $s^{\dagger}$ in terms
temperature ($e=\frac{3}{2}k_{B}\text{\sout{$\rho$}}_{T}T$ in MaxEnt):

$s=\frac{5}{2}\text{\sout{$\rho$}}_{T}k_{B}+k_{B}\left[\text{\sout{$\rho$}}_{a}\ln\left(\frac{m_{a}^{3/2}}{\text{\sout{$\rho$}}_{a}}\right)+\text{\sout{$\rho$}}_{b}\ln\left(\frac{m_{b}^{3/2}}{\text{\sout{$\rho$}}_{b}}\right)\right]+k_{B}\text{\sout{$\rho$}}_{T}\frac{3}{2}\ln\left(\frac{2\pi}{h^{2}}k_{B}T\right).$

It is a straightforward calculation to show that total energy is conserved via

$\partial_{t}e=s^{\dagger}\partial_{t}s+\rho_{a}^{\dagger}\partial_{t}\rho_{a}+\rho_{b}^{\dagger}\partial_{t}\rho_{b},$

as energy change occurs only via fluxes.

###  7.2 Projection to CIT in conjugate variables

An alternative method of extension is via a projection from the following
binary mixture state variables $\rho_{a},\leavevmode\nobreak\
\rho_{b},\leavevmode\nobreak\ \mathbf{u}_{a},\leavevmode\nobreak\
\mathbf{u}_{b},\leavevmode\nobreak\ s_{a}^{\dagger},\leavevmode\nobreak\
s_{b}^{\dagger}$ to the lower level $\rho_{a},\leavevmode\nobreak\
\rho_{b},\leavevmode\nobreak\ s$ where the projection is such that temperature
differences disappear, $\sigma^{\dagger}=0$, while reaching again mechanical
equilibrium.

This approach has the advantage of generality (it is not necessary to specify
entropy relation or look for the leading order solutions) however, as we shall
see, the resulting equations are somewhat more complex.

First, we need to rewrite the description of binary mixture to the new choice
of state variables. For this we use

$E^{*}=-E+s_{a}^{\dagger}s_{a}+s_{b}^{\dagger}s_{b},$

a free-energy. If we denote $\rho_{a}^{\ddagger}=\left(\frac{\partial
e^{*}}{\partial\rho_{a}}\right)_{\rho_{b},\mathbf{u}_{a},\mathbf{u}_{b},s_{a}^{\dagger},s_{b}^{\dagger}}$
then

$\displaystyle\rho_{\alpha}^{\ddagger}$
$\displaystyle=-\rho_{\alpha}^{\dagger},$
$\displaystyle\mathbf{u}_{\alpha}^{\ddagger}$
$\displaystyle=-\mathbf{u}_{\alpha}^{\dagger},$
$\displaystyle\left(s_{\alpha}^{\dagger}\right)^{\ddagger}$
$\displaystyle=s_{\alpha}.$

The evolution equations for conjugate variables are

$\displaystyle\partial_{t}s_{a}^{\dagger}$
$\displaystyle=\frac{1}{\det\mathrm{Hess}(e^{*})}\Bigg{[}\frac{\partial^{2}e^{*}}{\partial(s_{b}^{\dagger})^{2}}\underbrace{\left(\partial_{t}s_{a}-\frac{\partial^{2}e^{*}}{\partial
s_{a}^{\dagger}\rho_{a}}\partial_{t}\rho_{a}-\frac{\partial^{2}ej}{\partial
s_{a}^{\dagger}\rho_{b}}\partial_{t}\rho_{b}-\frac{\partial^{2}ej}{\partial
s_{a}^{\dagger}\mathbf{u}_{a}}\partial_{t}\mathbf{u}_{a}-\frac{\partial^{2}ej}{\partial
s_{a}^{\dagger}\mathbf{u}_{b}}\partial_{t}\mathbf{u}_{b}\right)}_{\boxed{1}}-$
$\displaystyle-\frac{\partial^{2}e^{*}}{\partial s_{a}^{\dagger}\partial
s_{b}^{\dagger}}\underbrace{\left(\partial_{t}s_{b}-\frac{\partial^{2}ej}{\partial
s_{b}^{\dagger}\rho_{a}}\partial_{t}\rho_{a}-\frac{\partial^{2}ej}{\partial
s_{b}^{\dagger}\rho_{b}}\partial_{t}\rho_{b}-\frac{\partial^{2}ej}{\partial
s_{b}^{\dagger}\mathbf{u}_{a}}\partial_{t}\mathbf{u}_{a}-\frac{\partial^{2}ej}{\partial
s_{b}^{\dagger}\mathbf{u}_{b}}\partial_{t}\mathbf{u}_{b}\right)}_{\boxed{2}}\Bigg{]},$
$\displaystyle\partial_{t}s_{b}^{\dagger}$
$\displaystyle=\frac{1}{\det\mathrm{Hess}(e^{*})}\left[\frac{\partial^{2}e^{*}}{\partial(s_{b}^{\dagger})^{2}}\boxed{2}-\frac{\partial^{2}e^{*}}{\partial(s_{a}^{\dagger})^{2}}\boxed{1}\right].$

Instead of using MaxEnt we prescribe the projection as equal partial
temperatures on the lower lever, i.e. $\sigma^{\dagger}=T_{a}-T_{b}=0$ while
$s^{\dagger}=\frac{1}{2}(s^{\dagger}_{a}+s^{\dagger}_{b})$.

With the assumption of mechanical equilibrium, the lower level is given by
$\sigma^{\dagger}=0,\mathbf{u}_{a}=\mathbf{u}_{b}=0$ (i.e. no dissipation due
to the mutual friction or temperature equilibration) while the conjugate
$(\sigma^{\dagger})^{\ddagger}=\sigma$ variable is given by

$\displaystyle
0=\partial_{t}\sigma^{\dagger}=\mathrm{div}\left(\left(\frac{s}{2}+\sigma\right)\mathbf{u}_{a}^{\ddagger}\right)-\mathrm{div}\left(\left(\frac{s}{2}-\sigma\right)\mathbf{u}_{b}^{\ddagger}\right)+\mathrm{div}(\rho_{a}\mathbf{u}_{a}^{\ddagger})\left(\frac{\partial^{2}e^{*}}{\partial
s_{b}^{\dagger}\partial\rho_{a}}-\frac{\partial^{2}e^{*}}{\partial
s_{a}^{\dagger}\partial\rho_{a}}\right)+\mathrm{div}(\rho_{b}\mathbf{u}_{b}^{\ddagger})\left(\frac{\partial^{2}e^{*}}{\partial
s_{b}^{\dagger}\partial\rho_{b}}-\frac{\partial^{2}e^{*}}{\partial
s_{a}^{\dagger}\partial\rho_{b}}\right).$ (7.16a)

This constitutive relation can then be used to close the extended Maxwell-
Stefan model

$\displaystyle\boxed{\partial_{t}\rho_{\alpha}=\mathrm{div}(\rho_{\alpha}\mathbf{u}_{\alpha}^{\ddagger})},$
(7.16b)
$\displaystyle\boxed{0=\rho_{a}\nabla\rho_{a}^{\ddagger}+\zeta_{ab}(\mathbf{u}_{a}^{\ddagger}-\mathbf{u}_{b}^{\ddagger})}-(s/2+\sigma)\nabla
s^{\dagger},$ (7.16c)
$\displaystyle\boxed{0=\rho_{b}\nabla\rho_{b}^{\ddagger}+\zeta_{ab}(\mathbf{u}_{b}^{\ddagger}-\mathbf{u}_{a}^{\ddagger})}-(s/2-\sigma)\nabla
s^{\dagger},$ (7.16d)
$\displaystyle\partial_{t}s^{\dagger}=\frac{1}{\det\mathrm{Hess}(e^{*})}\left(\frac{\partial^{2}e^{*}}{\partial(s_{b}^{\dagger})^{2}}-\frac{\partial^{2}e^{*}}{\partial(s_{a}^{\dagger})^{2}}\right)\left[\mathrm{div}\left(\left(\frac{s}{2}+\sigma\right)\mathbf{u}_{a}^{\ddagger}\right)-\mathrm{div}(\rho_{a}\mathbf{u}_{a}^{\ddagger})\frac{\partial^{2}e^{*}}{\partial
s_{a}^{\dagger}\partial\rho_{a}}-\mathrm{div}(\rho_{b}\mathbf{u}_{b}^{\ddagger})\frac{\partial^{2}e^{*}}{\partial
s_{a}^{\dagger}\partial\rho_{b}}\right],$ (7.16e)

where $s^{\dagger}=T$ is temperature and $s$ is related to it via the free
energy $e^{*}$ as $s=(s^{\dagger})^{\ddagger}$.

##  8 Soret coefficient as appearing from the multiscale analysis

An important outcome of our multiscale Maxwell-Stefan models is the ability to
assess the Soret coefficient based on a detailed description. The Soret
coefficient is a coefficient between the density or composition gradient and
the temperature gradient in a steady state characterised by a zero mass flux

$\nabla\rho_{\alpha}=-L_{\rho_{\alpha}T}\nabla T,$ (8.1)

where $L_{\rho_{\alpha}T}$ is the Soret coefficient describing the thermal
diffusion of component $\alpha$ (positive Soret coefficient have such
components which concentrate in colder regions).

Although the Soret coefficient has been studied for a long time in non-
equilibrium thermodynamics (being an off-diagonal entry of the Onsager matrix
relating thermodynamic fluxes and forces [10]), it is much less explored than
for instance the diffusion coefficient. A comprehensive review can be found in
[22]. The Soret coefficient can be approximated by means of kinetic theory
[18, 43], where it comes from the collision integral or activation energies,
or by using non-equilibrium statistical mechanics [19]. Finally, a
comprehensive comparison between the theoretical results and experiments can
be found in [20], where two origins of Soret coefficient are identified: the
kinetic origin (different properties of collisions of the constituents), and
the thermodynamic origin (different attraction/repulsion between the
constituents). In this work, we show another possible origin – different
behavior in the two-temperature state (before relaxation of temperatures of
the constituents to their common value).

In our study, the flux-force relations provide (noting that the flux vanish as
long as one of the components disappears)

$\rho_{a}v_{a}^{i}=-\frac{\rho_{a}\rho_{b}}{\zeta_{ab}\rho}\left(\rho_{a}\partial_{i}\mu_{a}+s_{a}\partial_{i}T\right)=0,$

where the expression for $s_{a}$ depends on the chosen reduction. In the first
approximate reduction where we assume detailed knowledge of entropy, $s_{a}$
is given by $\tilde{s}_{a}$ while in the full reductions they follow from
MaxEnt estimate either in the direct or conjugate state variables.

By splitting the chemical potential gradient into the temperature gradient
contribution and chemical potential potential of species $a$ at a constant
temperature, we have

$(\partial_{i}\mu_{a})_{T}=-\left[\frac{\partial\mu_{a}}{\partial
T}+\frac{s_{a}}{\rho_{a}}\right]\partial_{i}T.$

It is instructive to compare the various reductions. First, let us consider
the approximative reduction in the first part of this paper, i.e.
$s_{a}=\tilde{s}_{a}=s_{a}\big{|}_{T_{a}=T}$. Hence, Soret coefficient
$L_{\rho_{\alpha}T}$ satisfies

$L_{\rho_{a}T}=\frac{\frac{\partial\mu_{\beta}}{\partial
T}+\frac{\tilde{s}_{\beta}}{\rho_{\beta}}}{\frac{\partial\mu_{\beta}}{\partial\rho_{a}}}=\frac{1}{\frac{\partial\mu_{\beta}}{\partial\rho_{a}}}\left[\frac{\partial\mu_{\beta}}{\partial
T}-\frac{\partial\mu_{\beta}}{\partial
T_{\beta}}\bigg{|}_{T_{\beta}=T}+\frac{1}{\rho_{\beta}}(\tilde{s}_{\beta}-s_{\beta}\big{|}_{T_{\beta}=T})\right]=\frac{1}{\frac{\partial\mu_{\beta}}{\partial\rho_{a}}}\left[\frac{\partial\mu_{\beta}}{\partial
T}-\frac{\partial\mu_{\beta}}{\partial
T_{\beta}}\bigg{|}_{T_{\beta}=T}\right].$ (8.2)

Hence, a nontrivial Soret coefficient is appearing only as an effect of a
nonzero difference between partial temperatures. If we consider the particular
example described above in Section 5.2, we obtain

$L_{\rho_{a}T}=\frac{m_{b}}{n_{a}(n_{a}+n_{b})^{2}T}\left(n_{b}\frac{(n_{a}+n_{b})^{3}\kappa_{b}}{\kappa_{a}+\kappa_{b}}-\frac{2\kappa_{a}n_{a}^{3}}{\frac{(1-\gamma_{b})c_{b}^{\textrm{V}}}{n_{b}}-\frac{2n_{a}^{2}(\kappa_{a}+\kappa_{b})}{(n_{a}+n_{b})^{3}}}\right).$

On the other hand, if we use the results of the full reduction, we obtain
direct expressions for $s_{a}$ in terms of the less detailed state variables
without a clear link to the necessity of having distinct partial temperatures.
Hence, our study suggests that a temperature difference on a microscale
manifests on a macroscale as a Soret effect of thermodiffusion, which provides
an alternative origin of the Soret coefficient.

##  9 Discussion and Conclusion

Transport models are ubiquitous, attract a lot of attention, and are widely
applicable, see for example the recent monograph [48]. Here we focused on the
extension of the successful diffusion Maxwell-Stefan model to the
nonisothermal case to obtain a minimalistic and yet plausible description of
thermodiffusion. To this end, we consider momentum exchange among the mixture
components and keep track of the total entropy. To do so, we adapt a
thermodynamic multiscale framework and employ reduction techniques to reach
the desired level of detail.

In order to derive the hierarchy of theories of mixtures with two momenta and
two temperatures, one momentum and one temperature, and finally without any
momenta (mechanical equilibrium), we start from the Liouville equation and
reduce this complex behavior by Hamiltonian reductions. This way, we obtain a
hierarchy of models for multiscale binary homogeneous mixtures interacting
with electrodynamics. In particular we obtain new terms in the equations for
binary mixtures that can not be obtained by the standard approach based on
balance equations.

The first approach, Sections 2-6, is the most intuitive one and easiest to
use. It extends the classical Maxwell-Stefan to the non-isothermal case and,
in addition, includes the effect of electric field via electro-chemical
potential. Moreover it is straightforwardly extendable to a multicomponent
mixture. However, the drawback is that in this intuitive approach we
circumvented the problem of reducing the number of state variables entailing
the presence of partial entropies in the final set of equations. Hence the
problem is not closed and some form of further approximation have to take
place.

The second approach, Section 7, invoking MaxEnt to properly carry out the
reduction to the lower level is more precise and yields a closed system of
evolution equations. However, in order to explicitly carry out the reduction
we needed to employ a particular choice of the upper entropy. If the ideal gas
mixture is not appropriate, one has to repeat the reduction. Additionally, the
temperatures of the two constituents are subdued to the lower state variables
and are, in general, not equal. Note that when dissipation is dominant in the
lower level evolution, one can show that the difference in temperatures
vanishes.

Finally, the last approach used again MaxEnt dynamic reduction to yield the
lower level evolution being the extension of Maxwell-Stefan equations. This
time, the MaxEnt is carried out in the conjugate rather than direct state
variables as to control the temperature difference of the two constituents
being zero. Additionally, the yielded evolution equations are general – one
can chose free energy as the last step and thus specify the evolution
equations.

Note that in principle all these approaches are extendable to a mixtures with
a higher number of constituents but in the latter two approaches we do not
include such extensions explicitly here. In addition, our approach reminds, at
least in principle, the model of Lorentz diffusive gas [1], where the need to
include momentum exchange among constituents resulted in evolution equations
for new state variables in the spirit of Extended irreversible thermodynamics.
The obtained governing equations are also Maxwell-Stefan-like equations.

We apply the obtained Maxwell-Stefan theory to a dusty gas model, which we
numerically solve to illustrate the novelty of the nonisothermal extension of
the theory. One can see such effects directly from the equations themselves
as, for example from (7.16) one can expect nonfickian behaviour even in binary
mixtures due to the presence of temperature gradients. Another application of
our models is a non-isothermal generalization of the dusty gas model, its
generalized thermodynamic efficiency analysis, and finally the identification
of a new origin of the Soret coefficient.

In future, we would like test the obtained equation in fast processes where
the inertial effects and new terms might play a significant role.

## Acknowledgments

MP was supported by project No. UNCE/SCI/023 of Charles University Research
program. VK and MP were also supported by the Czech Science Foundation
(project no. 20-22092S). PV was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – The Berlin Mathematics Research Center MATH+
(EXC-2046/1, project ID: 390685689).

## References

  * [1] X Alvarez and D Jou.  Extended entropy and irreversible thermodynamics of a lorentz diffusive gas.  Physica A: Statistical Mechanics and its Applications, 377(1):79–83, 2007.
  * [2] RJ Atkin and RE Craine.  Continuum theories of mixtures: basic theory and historical development.  The Quarterly Journal of Mechanics and Applied Mathematics, 29(2):209–244, 1976.
  * [3] RJ Bearman.  On Molecular Basis of Some Current Theories of Diffusion.  Journal of Physical Chemistry, 65(11):1961–1968, 1961\.
  * [4] A Bedford and Do S Drumheller.  Theories of immiscible and structured mixtures.  International Journal of Engineering Science, 21(8):863–960, 1983\.
  * [5] D. Bothe and W. Dreyer.  Continuum thermodynamics of chemically reacting fluid mixtures.  ArXiv e-prints, 2014.
  * [6] Ray M Bowen.  Theory of mixtures.  In AC Eringen, editor, Continuum Physics, volume 3. Academic Press, New York, 1976.
  * [7] H.B. Callen.  Thermodynamics: an introduction to the physical theories of equilibrium thermostatics and irreversible thermodynamics.  Wiley, 1960.
  * [8] F.F. Chen.  Introduction to plasma physics.  Plenum Press, 1974.
  * [9] Reint de Boer.  Trends in continuum mechanics of porous media.  Springer, Dordrecht, The Netherlands, 2005.
  * [10] S. R. de Groot and P. Mazur.  Non-equilibrium Thermodynamics.  Dover Publications, New York, 1984.
  * [11] Donald A Drew and Stephen L Passman.  Theory of multicomponent fluids, volume 135.  Springer Science & Business Media, 2006.
  * [12] J. Fuhrmann and contributors.  VoronoiFVM.jl: Finite volume solver for coupled nonlinear partial differential equations.  https://github.com/j-fu/VoronoiFVM.jl, 2019-2021.
  * [13] J. W. Gibbs.  Collected Works.  Longmans; Green and Comp. New York, 1984.
  * [14] William G Gray and Cass T Miller.  Introduction to the thermodynamically constrained averaging theory for porous medium systems.  Springer, 2014.
  * [15] M. Grmela, V. Klika, and M. Pavelka.  Reductions and extensions in mesoscopic dynamics.  Phys. Rev. E, 92(032111), 2015.
  * [16] Miroslav Grmela and Hans Christian Öttinger.  Dynamics and thermodynamics of complex fluids. I. Development of a general formalism.  Phys. Rev. E, 56:6620–6632, Dec 1997.
  * [17] Morton E Gurtin, Eliot Fried, and Lallit Anand.  The mechanics and thermodynamics of continua.  Cambridge University Press, 2010.
  * [18] J.O. Hirschfelder, C.F. Curtiss, R.B. Bird, and University of Wisconsin. Theoretical Chemistry Laboratory.  Molecular theory of gases and liquids.  Structure of matter series. Wiley, 1954.
  * [19] L. J. T. M. Kempers.  A thermodynamic theory of the soret effect in a multicomponent liquid.  J. Chem. Phys., 90:6541, 1989.
  * [20] L. J. T. M. Kempers.  A comprehensive thermodynamic theory of the soret effect in a multicomponent gas, liquid, or solid.  Journal of chemical physics, 115(14), 2001.
  * [21] Piet JAM Kerkhof and Marcel AM Geboers.  Analysis and extension of the theory of multicomponent fluid diffusion.  Chemical Engineering Science, 60(12):3129–3167, 2005.
  * [22] S. Kjelstrup and D. Bedeaux.  Non-Equilibrium Thermodynamics of Heterogeneous Systems.  Series on Advances in Statistical Mechanics. World Scientific, 2008.
  * [23] V. Klika.  A guide through available mixture theories for applications.  Critical Reviews in Solid State and Material Science, 2013.  Accepted Manuscript.
  * [24] Václav Klika.  Modelling of Biomaterials as an Application of the Theory of Mixtures, pages 105–139.  Springer International Publishing, Cham, 2021.
  * [25] Václav Klika and Andrew L Krause.  Beyond onsager–casimir relations: shared dependence of phenomenological coefficients on state variables.  The journal of physical chemistry letters, 9(24):7021–7025, 2018\.
  * [26] Václav Klika, Michal Pavelka, and Jay B Benziger.  Functional constraints on phenomenological coefficients.  Physical Review E, 95(2):022125, 2017.
  * [27] Václav Klika, Michal Pavelka, Petr Vágner, and Miroslav Grmela.  Dynamic maximum entropy reduction.  Entropy, 21(715), 2019.
  * [28] R Krishna and JA Wesselingh.  Review article number 50 - The Maxwell-Stefan approach to mass transfer.  Chemical Engineering Science, 52(6):861–911, MAR 1997\.
  * [29] G. Lebon, D. Jou, and J.C. Vázquez.  Understanding Non-Equilibrium Thermodynamics: Foundations, Applications, Frontiers.  SpringerLink: Springer e-Books. Springer London, Limited, 2008.
  * [30] E.A. Mason and A.P. Malinauskas.  Gas transport in porous media: the dusty-gas model.  Number v. 17 in Chemical engineering monographs. Elsevier, 1983.
  * [31] Lars Onsager.  Reciprocal relations in irreversible processes. I.  Phys. Rev., 37:405–426, Feb 1931.
  * [32] Lars Onsager.  Reciprocal relations in irreversible processes. ii.  Phys. Rev., 38:2265–2279, Dec 1931.
  * [33] Hans Christian Öttinger and Miroslav Grmela.  Dynamics and thermodynamics of complex fluids. II. Illustrations of a general formalism.  Phys. Rev. E, 56:6633–6655, Dec 1997.
  * [34] H.C. Öttinger.  Beyond Equilibrium Thermodynamics.  Wiley, 2005.
  * [35] M. Pavelka, V. Klika, O. Esen, and M. Grmela.  A hierarchy of Poisson brackets in non-equilibrium thermodynamics.  Physica D: Nonlinear phenomena, 335:54–69, 2016.
  * [36] M. Pavelka, V. Klika, and M. Grmela.  Time reversal in nonequilibrium thermodynamics.  Phys. Rev. E, 90(062131), 2014.
  * [37] Michal Pavelka, Václav Klika, and Miroslav Grmela.  Time reversal in nonequilibrium thermodynamics.  Phys. Rev. E, 90:062131, Dec 2014.
  * [38] Michal Pavelka, Václav Klika, and Miroslav Grmela.  Multiscale Thermo-Dynamics.  de Gruyter (Berlin), 2018.
  * [39] Michal Pavelka, Václav Klika, Petr Vágner, and František Maršík.  Generalization of exergy analysis.  Applied Energy, 137(0):158 – 172, 2015.
  * [40] Michal Pavelka, František Maršík, and Václav Klika.  Consistent theory of mixtures on different levels of description.  International Journal of Engineering Science, 78(0):192 – 217, 2014\.
  * [41] Kumbakonam R Rajagopal and Luoyi Tao.  Mechanics of mixtures, volume 35.  World scientific, 1995.
  * [42] O. Souček, V. Průša, J. Málek, and K. R. Rajagopal.  On the natural structure of thermodynamic potentials and fluxes in the theory of chemically non-reacting binary mixtures.  Acta Mechanica, 2013.  Accepted manuscript.
  * [43] L. J. Tichacek, W. S. Kmak, and H. G. Drickamer.  Thermal diffusion in liquids: the effect of non-ideality and association.  J. Phys. Chem., 60:660, 1956.
  * [44] C Truesdell.  Rational Thermodynamics.  Springer-Verlag, New York, 1984.
  * [45] P. Vágner and J. Fuhrmann.  Finite-volume implementation of non isothermal dusty gas model.  https://github.com/fafroo/Nonisothermal-dusty-gas-model/, 2022.
  * [46] Petr Vágner, Michal Pavelka, and Oğul Esen.  Multiscale thermodynamics of charged mixtures.  Continuum Mechanics and Thermodynamics, Jul 2020.
  * [47] Petr Vágner, Michal Pavelka, and František Maršík.  Pitfalls of exergy analysis.  Journal of Non-Equilibrium Thermodynamics, 42(2):201–216, 2017\.
  * [48] David C Venerus and Hans Christian Öttinger.  A modern course in transport phenomena.  Cambridge University Press, 2018.

##  Appendix A MaxEnt is equivalent to MinEne

Let us now recall that maximization of concave entropy $S(E,x)$ with respect
to $x$ is equivalent to minimization of energy $E(S,x)$, as shown for instance
in [7, 13].

Assume a stationary MaxEnt point

$\left(\frac{\partial S}{\partial x}\right)_{E}\Big{|}_{\bar{x},E}=0.$ (A.1)

By the implicit function theorem one gets

$\left(\frac{\partial E}{\partial x}\right)_{S}=-\frac{\left(\frac{\partial
S}{\partial x}\right)_{E}}{\left(\frac{\partial S}{\partial E}\right)_{x}}=0$
(A.2)

and thus

$\left(\frac{\partial E}{\partial
x}\right)_{S}\Big{|}_{\bar{x},S(E,\bar{x})}=0.$ (A.3)

Energy thus also has a stationary point at $\bar{x}$.

Positivity of temperature means that entropy grows as energy grows (keeping
$x$ constant).

Imagine now a the surface $\sigma$ described by $S(x,E)=\mathrm{const}$ within
the three-dimensional space $(E,x,S)$. Plane $\sigma_{E}$ is a plane
satisfying $E=\mathrm{const}$. The point with $\bar{x}$ is the point with the
highest value $S$ from the curve $\gamma_{E}=\sigma\cap\sigma_{E}$. Plane
$\sigma_{S}$ is the plane characterized by $S=\mathrm{const}$ containing the
point $\bar{x}$. Because entropy grows with energy, the curve
$\gamma_{S}=\sigma_{S}\cap\sigma$ has its point with lowest energy located at
$\bar{x}$ while other points of that curve have higher energy. Curves
$\gamma_{E}$ and $\gamma_{S}$ have only one common point, namely $\bar{x}$.
Maximization of entropy along the curve $\gamma_{E}$ is equivalent to
miminization of energy along the curve $\gamma_{S}$. Therefore, MaxEnt is
equivalent with MinEne for convex entropies.

##  Appendix B Finite volume discretization of binary DGM

We summarize the binary, neutral, non-isothermal dusty gas model for ideal
gasses described by Equations (6.2), (6.4), and (6.5). It is given by

$\displaystyle\partial_{t}n_{\alpha}+\partial_{i}J_{\alpha}^{i}$
$\displaystyle=0\leavevmode\nobreak\ ,\quad\mbox{for
}\alpha\in\\{a,b\\}\leavevmode\nobreak\ ,$ (B.1a)
$\displaystyle\partial_{t}\left(T\left(c_{a}^{\textrm{V}}n_{a}+c_{b}^{\textrm{V}}n_{b}+c_{D}^{\textrm{V}}n_{D}\right)\right)+\partial_{i}\left(-\lambda\partial_{i}T+T\left(\gamma_{a}c_{a}^{\textrm{V}}J_{a}^{i}+\gamma_{b}c_{b}^{\textrm{V}}J_{b}^{i}\right)\right)$
$\displaystyle=0\leavevmode\nobreak\ .$ (B.1b) where the force-flux relation
are given as follows,
$\displaystyle\begin{bmatrix}\partial_{i}n_{a}+n_{a}\partial_{i}\ln T\\\
\partial_{i}n_{b}+n_{b}\partial_{i}\ln T\\\
\end{bmatrix}=\underbrace{\begin{bmatrix}-\left(\frac{1}{n_{a}+n_{b}}\frac{n_{b}}{\frac{\varepsilon}{\tau}D_{ab}}+\frac{1}{\frac{\varepsilon}{\tau}D_{aD}}\right)&\frac{n_{a}}{n_{a}+n_{b}}\frac{1}{\frac{\varepsilon}{\tau}D_{ab}}\\\
\\\
\frac{n_{b}}{n_{a}+n_{b}}\frac{1}{\frac{\varepsilon}{\tau}D_{ab}}&-\left(\frac{1}{\frac{\varepsilon}{\tau}D_{bD}}+\frac{1}{n_{a}+n_{b}}\frac{n_{a}}{\frac{\varepsilon}{\tau}D_{ab}}\right)\\\
\end{bmatrix}}_{\mathbf{M}}\begin{bmatrix}J^{i}_{a}\\\ J^{i}_{b}\\\
\end{bmatrix}\leavevmode\nobreak\ .$ (B.1c)

Let us assume a Voronoi tessellation $\\{\omega_{k}\\}_{k=1}^{N}$ of the
domain $\Omega$. For a control volume $\omega_{k}$, the discrete values of the
densities and temperature in the collocation point
$\mathbf{x_{k}}\in\omega_{k}$ are defined as

$\displaystyle n_{\alpha,k}$
$\displaystyle:=\frac{1}{|\omega_{k}|}\int_{\omega_{k}}n_{\alpha}\mathrm{d}\mathbf{x}\leavevmode\nobreak\
,\quad\alpha\in\\{a,b\\}\leavevmode\nobreak\ ,$ (B.2) $\displaystyle T_{k}$
$\displaystyle:=\frac{1}{|\omega_{k}|}\int_{\omega_{k}}T\mathrm{d}\mathbf{x}\leavevmode\nobreak\
.$ (B.3)

Integration over a control volume $\omega_{k}$ and application of the
divergence theorem gives the discrete mass balance which reads

$\displaystyle|\omega_{k}|\partial_{t}n_{\alpha,k}+\sum_{\sigma_{kl}\neq\emptyset}|\sigma_{kl}|J^{\textrm{num}}_{\alpha,kl}$
$\displaystyle=0\leavevmode\nobreak\ ,\quad\alpha\in\\{a,b\\}$ (B.4)
$\displaystyle|\omega_{k}|\partial_{t}\left(T_{k}\left(c^{V}_{a}n_{a,k}+c^{V}_{b}n_{b,k}+c^{V}_{D}n_{D}\right)\right)+\sum_{\sigma_{kl}\neq\emptyset}|\sigma_{kl}|J^{\textrm{num}}_{T,kl}$
$\displaystyle=0\leavevmode\nobreak\ ,$ (B.5)

where $J^{\textrm{num}}_{\alpha,kl}$ is the numerical flux of density
$n_{\alpha}$ across the edge
$\sigma_{kl}=\partial\omega_{k}\cap\partial\omega_{k}$;
$J^{\textrm{num}}_{T,kl}$ is the numerical flux The discrete heat flux is
defined as

$\displaystyle
J^{\textrm{num}}_{T,kl}=\frac{\lambda}{h_{kl}}\left[T_{l}B(Q_{T,kl})-T_{k}B(-Q_{T,kl})\right]$
(B.6) $\displaystyle\text{where}\quad
Q_{T,kl}=\frac{1}{h_{kl}\lambda}\sum_{\alpha=a,b}\gamma_{\alpha}c^{V}_{\alpha}J^{\textrm{num}}_{\alpha,kl}\leavevmode\nobreak\
.$

Here $B(x)=\frac{x}{e^{x}-1}$ is the Bernoulli function and
$h_{kl}=|\mathbf{x}_{k}-\mathbf{x}_{l}|$ denotes the distance of the
collocation points.

$\displaystyle F^{\textrm{num}}_{\alpha,kl}$
$\displaystyle=\frac{1}{h_{kl}}\left(n_{\alpha,l}B(Q_{kl})-n_{\alpha,k}B(-Q_{kl})\right)\leavevmode\nobreak\
,\quad\alpha\in\\{a,b\\}\leavevmode\nobreak\ ,$ (B.7)
$\displaystyle\text{where}\quad Q_{kl}$ $\displaystyle=-\frac{1}{h_{kl}}(\log
T_{l}-\log T_{k})\leavevmode\nobreak\ ,\quad\text{for}\leavevmode\nobreak\
k=1,\dots,N\leavevmode\nobreak\ ,\leavevmode\nobreak\ .$

We assume the following discrete version of the mobility matrix $\mathbf{M}$

$\displaystyle
M^{\textrm{num}}_{kl}=\frac{\tau}{\varepsilon}\begin{bmatrix}-\left(\frac{\bar{n}_{b,kl}}{\bar{n}_{kl}}\frac{1}{D_{ab}}+\frac{1}{D_{aD}}\right)&\frac{\bar{n}_{a,kl}}{\bar{n}_{kl}}\frac{1}{D_{ab}}\\\
\\\
\frac{\bar{n}_{b,kl}}{\bar{n}_{kl}}\frac{1}{D_{ab}}&-\left(\frac{\bar{n}_{a,kl}}{\bar{n}_{kl}}\frac{1}{D_{ab}}+\frac{1}{D_{bD}}\right)\\\
\end{bmatrix}$ (B.8)

where
$\bar{n}_{\alpha,kl}:=\tfrac{1}{2}\left(n_{\alpha,k}+n_{\alpha,l}\right)$ and
$\bar{n}_{kl}:=\bar{n}_{a,kl}+\bar{n}_{b,kl}$. The linear system,

$\displaystyle
F^{\textrm{num}}_{kl}=M^{\textrm{num}}_{kl}J^{\textrm{num}}_{kl}\leavevmode\nobreak\
,$ (B.9)

gives the discrete mass fluxes.

###  B.1 Numerical solution in 1D

We wrote a Pluto.jl notebook which shows a solution of the discrete system
(B.3) endowed with Dirichlet boundary conditions on an interval. The solution
of the discrete system itself is solved by the finite volume package
Voronoi.jl [12]. Additionally, we coded a parallel solution with analytically
inverted $M$, that is
$J^{\textrm{num}}_{kl}=\left(M^{\textrm{num}}_{kl}\right)^{-1}F^{\textrm{num}}_{kl}$,
to check the implicitly calculated fluxes in equation (B.9). The two
approaches agree, the notebook can be downloaded here [45].

Figure 2: Relaxation of system (B.1) after temperature step of the boundary
condition ($\Delta T=130\textrm{K}$). The Dirichlet boundary conditions for
the pressure $p$ and the molar fraction $x_{a}$ are constant and equal for the
both boundaries.