# Structural electroneutrality in Onsager–Stefan–Maxwell transport with
charged species

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

###### Abstract

We present a method to embed local electroneutrality within
Onsager–Stefan–Maxwell electrolytic-transport models, circumventing their
formulation as differential systems with an algebraic constraint. Flux-
explicit transport laws are formulated for general multicomponent
electrolytes, in which the conductivity, component diffusivities, and
transference numbers relate to Stefan–Maxwell coefficients through invertible
matrix calculations. A construction we call a ‘salt–charge basis’ implements
Guggenheim’s transformation of species electrochemical potentials into
combinations describing a minimal set of neutral components, leaving a unique
combination associated with electricity. Defining conjugate component
concentrations and fluxes that preserve the structures of the Gibbs function
and energy dissipation retains symmetric Onsager reciprocal relations. The
framework reproduces Newman’s constitutive laws for binary electrolytes and
the Pollard–Newman laws for molten salts; we also propose laws for salt
solutions in two-solvent blends, such as lithium-ion-battery electrolytes.
Finally, we simulate a potentiostatic Hull cell containing a non-ideal binary
electrolyte with concentration-dependent properties.

##  1 Introduction

The growing impetus for improved electrochemical energy storage makes it vital
to understand ion transport. Concentrated solution theory [1, 2] is widely
used to model mass and heat transport in electrolytic materials on the
continuum scale. One of the most popular implementations is Doyle, Fuller and
Newman’s landmark physics-based model of the dual-insertion lithium-ion
battery, which incorporates a concentrated-solution model of a binary
electrolyte, comprising a simple salt dissolved in an uncharged solvent [3,
4].

Concentrated solution theory is grounded in both classical thermodynamics and
linear irreversible thermodynamics, with a model structure that sharply
distinguishes between equilibrium properties and transport properties. It
derives from the Onsager–Stefan–Maxwell transport laws [5, 1], written for $n$
species in an isotropic, isothermal, isobaric common phase as

$-\nabla\mu_{i}=\sum_{j=1}^{n}M_{ij}\vec{N}_{j},$ (1.1)

where $\mu_{i}$ signifies the electrochemical potential of species
$i\in\left\\{1,...,n\right\\}$ and $\vec{N}_{i}$, its total molar flux. The
transport coefficients make up an $n\times n$ matrix M whose entries are

$M_{ij}=\begin{cases}-\dfrac{RT}{c_{\text{T}}\mathscr{D}_{ij}}\>\>\>\text{if}\>\>i\neq
j\vspace{6pt}\\\ \dfrac{RT}{c_{\text{T}}}\displaystyle{\sum_{k\neq
i}^{n}}\dfrac{c_{k}}{\mathscr{D}_{ik}c_{j}}\>\>\>\text{if}\>\>i=j,\end{cases}$
(1.2)

in which $R$ represents the gas constant, $T$, the absolute temperature, and
$c_{\text{T}}$, the total species molarity; $c_{i}$ is the molarity of species
$i$, and $\mathscr{D}_{ij}$, the Stefan–Maxwell diffusivity of species $i$
through species $j$. It has been shown for solutions of uncharged constituents
that Onsager’s reciprocal relations and the second law of thermodynamics
require M to be symmetric positive semidefinite [6]. Also, M affords a single
null eigenvector $\textbf{c}=[c_{1},c_{2},...,c_{n}]^{\top}$, a further
structure which ensures that convection of material at the bulk velocity does
not dissipate energy [7]. (In this work bold lower-case letters always
indicate column matrices.) For fluids wholly made up of uncharged species, the
Onsager–Stefan–Maxwell equations, when combined with an equation that
establishes a convective velocity, have been shown to constitute a closed
model, neither underdetermined nor overdetermined numerically [8].

An important practicality emerges when considering transport in liquids such
as electrolytic solutions, which conduct electricity via the motion of ions.
When species that carry charge are mobile, their thermal motion screens local
charge imbalances on a length scale comparable to the Debye length (typically
of the order of 1 nm) [9]. Thus, if volume elements in a simulation have
characteristic size much larger than the Debye length, the local excess charge
density $\rho_{\textrm{e}}$ within them can be assumed to vanish locally, that
is,

$\rho_{\textrm{e}}=F\sum_{j=1}^{n}z_{j}c_{j}=0,$ (1.3)

where $F$ is Faraday’s constant and $z_{i}$ is the equivalent charge of
species $i$. Local electroneutrality is foundational to concentrated solution
theory and has been assumed nearly ubiquitously in its applications.

Adopting the local electroneutrality approximation simplifies transport models
by allowing one charged-species concentration to be eliminated from the
governing equations, both reducing the complexity of thermodynamic
parametrization and streamlining calculations. Electroneutrality makes
concentrated solution theory hard to reconcile with the Onsager–Stefan–Maxwell
theory for nonelectrolytes, however, because the additional algebraic
constraint imposed by equation (1.3) appears to induce an overdetermination
within equations (1.1). In particular, it is not clear whether the symmetry
and rank of the transport matrix are preserved within electrolytic systems
after electroneutrality restricts the available space of species compositions.

Newman and colleagues carefully circumvented the issues raised by
electroneutrality when formulating models of binary electrolytes and salt
melts [2, 10, 11]. These analyses emphasized the importance of separating
terms into neutral and non-neutral contributions: individual ion
concentrations were lumped into salt concentrations; species fluxes were
expressed in terms of salt chemical potentials and current density; and terms
associated with charge density and charge convection were explicitly
neglected.

Possible applications of electroneutral Onsager–Stefan–Maxwell theory readily
exceed the scope of the few cases that have been considered to date. For
example, aqueous vanadium flow batteries—a transport system comprising at
least eight distinct species—have recently been modelled [12, 13], but it is
not immediately clear how such a model could be extended to account for
additional components in a consistent way. At the same time, the increased
demand for battery simulations that account for electrolyte additives or
contaminants, solution-phase reaction intermediates, or products of side
reactions raises the need for a general, numerically robust model formulation.

This paper generalizes Newman’s particular implementations of electroneutral
concentrated solution theory for binary solutions and salt melts. We implement
Guggenheim’s thermodynamic programme [14] by creating a structure-preserving,
invertible linear transformation that maps the free-energy changes incurred by
individual charged species into changes involving uncharged combinations of
species, assembling the species electrochemical potentials into a basis set of
component chemical potentials and a single quantity that can be interpreted as
an electric potential. The same transformation regroups charged-species fluxes
into component fluxes and current density. As well as reconciling the
electroneutrality approximation with the fundamental equations of linear
irreversible thermodynamics in a structure-preserving way, the theory produces
a general set of transference numbers and an ionic conductivity that match
definitions given for binary electrolytes by Newman, as well as leading to a
naturally symmetric set of component diffusivities.

Our development enables the extension of Newman’s flux-explicit transport laws
for binary electrolytes to general multicomponent electrolytes containing any
number of charged or uncharged species, while maintaining the crucial benefits
of thermodynamic consistency, convection-invariance of diffusion, and
symmetric reciprocal relations inherited from the Onsager–Stefan–Maxwell
equations. A numerical simulation of a Hull cell that leverages a solver for
Onsager–Stefan–Maxwell transport illustrates the practical utility of this
structure-based approach to local electroneutrality within concentrated
solution theory.

##  2 Construction of a salt–charge basis

The key mathematical concepts employed below are linear transformations and
inner products, typically represented by matrix operations. These tools are
used to develop a change of composition basis from one that lists charged-
species molarities into one that lists neutral-salt molarities alongside
charge density. Instead of beginning with species concentrations, however, the
new composition basis is identified by following Guggenheim’s process of
grouping species electrochemical potentials into new quantities that represent
electrically neutral combinations of species [14]. For any given electrolytic
material, one can establish a minimal set of such combinations by positing a
set of independent hypothetical chemical reactions that we call _simple
association equilibria_ , which leave uncharged species alone or bring charged
species together to form neutral salts. We refer to products of the simple
association equilibria as _components_ , rather than species.

Every electrochemical system is globally electroneutral at equilibrium and
every electrochemical reaction necessarily balances both charge and atoms
between its reactants and products. Moreover, any physically realizable
electrolytic solution must contain at least two countercharged species. It
follows that the $n$ species comprising any electrolyte combine to form at
most $n-1$ distinct uncharged components. Thus, association equilibria can be
used to combine ions into a set of uncharged products of association
reactions, whose energetics in solution is described unambiguously in terms of
classical chemical potentials. This idea dates back to the work of Guggenheim
almost a century ago [14]; the process we outline in this paper is merely
intended to systematize Guggenheim’s framework, and consequently, to enable an
analysis of electrolyte thermodynamics with the tools of modern linear
algebra. We will call the elementary structure that establishes the species
charges and a minimal set of hypothetical reactions that neutralize them a
salt–charge basis.

It is helpful to illustrate with an example before formalizing the process for
defining a salt–charge basis. Consider a solution of $\ch{Na+}$, $\ch{Cl-}$,
$\ch{Mg}^{2+}$ and $\ch{SO_{4}^{2-}}$ in water ($\ch{H2O}$). Among these five
species, one can write at most $5-1=4$ independent reactions to form
electrically neutral entities. For instance, one could write a trivial
equilibrium between the naturally uncharged water molecule and itself, and
then write three equilibria that combine distinct countercharged pairs of ions
to neutralize their charges:

$\displaystyle\ch{H2O}\leftrightharpoons\ch{H2O}$
$\displaystyle\ch{Na+}+\ch{Cl-}\leftrightharpoons\ch{NaCl}$
$\displaystyle\ch{Mg}^{2+}+2\ch{Cl-}\leftrightharpoons\ch{MgCl2}$
$\displaystyle 2\ch{Na+}+\ch{SO_{4}^{2-}}\leftrightharpoons\ch{Na2SO4}.$

Along with the list of species charges $\left\\{0,1,-1,2,-2\right\\}$, these
simple association equilibria define a salt–charge basis for this five-species
aqueous solution, in this case involving the uncharged components $\ch{H2O}$,
$\ch{NaCl}$, $\ch{MgCl2}$, and $\ch{Na2SO4}$.

The choice of reactions to form a salt–charge basis is not generally unique,
but any neutral component not defined can always be recovered by writing
additional reactions involving products formed by the selected simple
association equilibria. In the example system, observe that one can create the
fourth possible simple binary salt, $\ch{MgSO_{4}}$, via the recombination
reaction

$\ch{MgCl2}+\ch{Na2SO4}\leftrightharpoons 2\ch{NaCl}+\ch{MgSO4}.$ (2.1)

Through this equilibrium, any thermodynamic property of $\ch{MgSO4}$ within
this system—such as its chemical potential, partial molar volume, or partial
molar entropy—is determined by the properties of dissolved $\ch{NaCl}$,
$\ch{Na2SO4}$, and $\ch{MgCl2}$. The energetics of any multi-ion salt with
higher formula-unit stoichiometry also follows from the simple association
equilibria. Consider

$\ch{MgCl2}+\ch{Na2SO4}\leftrightharpoons\ch{NaCl}+\ch{MgNaClSO4},$

which forms a more complex salt, $\ch{MgNaClSO4}$, from all four available
ions.

The formal procedure for creating a salt–charge basis generalizes the
preceding example. Given a set of species that comprise an electrolytic
solution, start by placing them in an ordered list, wherein each species’
position assigns its numerical index. Let $n_{\textrm{c}}$ denote the number
of charged species. List all $(n-n_{\textrm{c}})$ uncharged species first, and
charged species afterward. For consistency with prior approaches, ensure that
the $n$th and $(n-1)$th species have opposite charges [15, 16]. Then let
$\textbf{z}=[z_{1},...,z_{n}]^{\top}$ be a column matrix comprising the
species equivalent charges in their assigned order. The ordering ensures that
entries $z_{1},...,z_{n-n_{\textrm{c}}}$ of z equal zero, entries
$z_{n-n_{\textrm{c}}+1},...,z_{n}$ of z are nonzero, and $z_{n}/z_{n-1}<0$.

Next, write a set of $n-1$ simple association equilibria involving the
species, beginning with $n-n_{\textrm{c}}$ trivial reactions for the naturally
uncharged species, and following with $n_{\textrm{c}}-1$ linearly independent
association reactions, within which countercharged pairs of charged species
combine to form simple binary salts. (Here ‘simple’ means that the reactant
stoichiometric coefficients are coprime positive integers.) Let column matrix
$\bm{\nu}_{i}=[\nu_{i1},\nu_{i2},...,\nu_{in}]^{\top}$ contain the
stoichiometric coefficients $\nu_{ij}$ of reactant species $j$ in the $i$th
simple association equilibrium. For uncharged species
$k=1,...,n-n_{\textrm{c}}$, stoichiometric column $\bm{\nu}_{k}$ will be the
$k$th column of the identity matrix, notated henceforth as $\textbf{i}_{k}$.
For $k=n-n_{\textrm{c}}+1,...,n-1$, $\bm{\nu}_{k}$ has two positive-integer
entries, corresponding to the stoichiometric coefficients of species in a
formula unit of the simple binary salt they form. The restriction to simple
salts ensures that the total number of ions in each formula unit is minimal.

Guggenheim’s condition for association equilibria requires that the ion
stoichiometry within a formula unit of salt always serves to neutralize its
charge [17]. In terms of reactant stoichiometry in the $i$th simple
association equilibrium, the condition is expressed as

$\bm{\nu}_{i}^{\top}\textbf{z}=0,$ (2.2)

where superscript $\top$ indicates the matrix transpose. Note that this
condition also holds for every trivial simple association equilibrium based on
an uncharged reactant.

It is useful to view equation (2.2) as a statement about a matrix inner
product, namely: the Guggenheim condition requires that each column of
stoichiometric coefficients $\bm{\nu}_{i}$ is orthogonal to the charges z.
Since the set of stoichiometric columns is also linearly independent by
construction, the reactant stoichiometries from the $n-1$ simple association
equilibria establish an $(n-1)$-dimensional linear space orthogonal to the
span of z. Moreover, the $n$-tuple
$\left\\{\bm{\nu}_{1},...,\bm{\nu}_{n-1},\textbf{z}\right\\}$ is linearly
independent, and consequently forms a basis for the entire $n$-dimensional
composition space. This set of column matrices formally represents the chosen
salt–charge basis.

##  3 Chemical potential in a salt–charge basis

Selecting a salt–charge basis enables one to construct a linear transformation
that maps the $n_{\textrm{c}}$ charged-species electrochemical potentials into
$n_{\textrm{c}}-1$ chemical potentials of neutral salts, leaving a unique
additional combination of electrochemical potentials that can be used to
describe charge. The transformation is represented by an $n\times n$ matrix Z,

$\textbf{Z}=\left[\begin{array}[]{c}\bm{\nu}_{1}^{\top}\\\ \vdots\\\
\bm{\nu}_{n-1}^{\top}\\\
{\textbf{z}^{\top}}/{\|\textbf{z}\|}\end{array}\right],$ (3.1)

in which $\|\textbf{v}\|=\sqrt{\textbf{v}^{\top}\textbf{v}}$ for any column v.
Its rows are linearly independent, so Z is invertible. Each of the first $n-1$
rows is orthogonal to the last ($n$th) row, and each of the first
$n-n_{\textrm{c}}$ rows, which correspond to uncharged species, is orthogonal
to every other row. Bear in mind that internal rows corresponding to salts,
that is, rows $i=n-n_{\textrm{c}}+1,...,n-1$ of Z, are not necessarily
mutually orthogonal.

The transformation Z combines species electrochemical potentials into
component chemical potentials as follows. Let the column matrix
$\bm{\mu}=[\mu_{1},...,\mu_{n}]^{\top}$ contain the species electrochemical
potentials in order. Then multiplication by Z maps $\bm{\mu}$ into new
coordinates $\bm{\mu}_{Z}$ that stand over the chosen salt–charge basis. In
the form of a matrix equation, $\bm{\mu}_{Z}$ is defined as

$\bm{\mu}_{Z}=\textbf{Z}\bm{\mu}.$ (3.2)

The first $n-1$ rows of $\bm{\mu}_{Z}$ so formed comprise the chemical
potentials of a minimal number of neutral components, with the $i$th row
representing the chemical potential of the product formed by the $i$th simple
association equilibrium used to define the salt–charge basis. The last row
quantifies all of the information in the set of species electrochemical
potentials associated with reversible processes that do not conserve charge.
It is convenient to subdivide $\bm{\mu}_{Z}$, as

$\bm{\mu}_{Z}=\left[\begin{array}[]{c}\bm{\mu}_{\nu}\\\
\mu_{z}\end{array}\right],$ (3.3)

where the ($n-1$)-dimensional column $\bm{\mu}_{\nu}$ contains the $n-1$
component chemical potentials. More explicitly, if $\bm{\nu}_{k}$ has two
nonzero entries,
$\mu_{\nu,k}=(\bm{\mu}_{\nu})_{k}=\bm{\nu}_{k}^{\top}\bm{\mu}$ is a salt
chemical potential formed by the same process used by Newman and Thomas–Alyea
[1]; otherwise (if $\bm{\nu}_{k}=\textbf{i}_{k}$), $\mu_{\nu,k}$ reproduces
the chemical potential expected for a naturally uncharged species. As we shall
see later, the remaining quantity, $\mu_{z}$, leads naturally to a definition
of solution voltage.

##  4 Concentration and flux in a salt–charge basis

The transformation to a salt–charge basis embodied by Z can be exploited to
form quantities which describe the amounts and motion of uncharged components.
Natural definitions of component concentrations and component fluxes arise
heuristically from an invariance argument.

We assert as a physical principle that the elementary laws from thermodynamics
and irreversible thermodynamics should retain their structures when expressed
in terms of components, rather than species. The essential thermodynamic
relation among species is the Euler equation,

$\tilde{G}=\sum_{i=1}^{n}\mu_{i}c_{i},$ (4.1)

which defines volumetric Gibbs free energy $\tilde{G}$ in terms of species
molarities. The essential relation from irreversible thermodynamics is the
dissipation function, written at a point within an isothermal, isobaric single
phase as

$T\dot{s}=\sum_{i=1}^{n}\vec{N}_{i}\cdot\left(-\nabla\mu_{i}\right),$ (4.2)

which relates local energy loss $T\dot{s}$ to the species fluxes.

To identify a natural set of composition variables over a given salt–charge
basis, require $\tilde{G}$ to retain its structure whether written in terms of
species or components. Then the matrix relationship

$\tilde{G}=\bm{\mu}^{\top}\textbf{c}=\bm{\mu}_{Z}^{\top}\textbf{c}_{Z}$ (4.3)

expresses $\tilde{G}$ either in terms of the species concentrations c or a
column matrix $\textbf{c}_{Z}$ that lists component concentrations in the
chosen salt–charge basis. Since Z is nonsingular, it follows from the second
equality in equation (4.3) and equation (3.2) that

$\textbf{c}_{Z}=\textbf{Z}^{-\top}\textbf{c},$ (4.4)

where superscript ${-\top}$ indicates the inverse transpose.

Introducing Faraday’s law helps to show how the column matrix $\textbf{c}_{Z}$
parses electrical information. Faraday’s law for charge combines with equation
(4.4) to show that excess charge density can be expressed in two equivalent
ways,

$\rho_{\textrm{e}}=F\textbf{z}^{\top}\textbf{c}~{}~{}\textrm{or}~{}~{}\rho_{\textrm{e}}=F\left(\textbf{Z}\textbf{z}\right)^{\top}\textbf{c}_{Z},$
(4.5)

in which the distributive property of the matrix transpose was used to
rearrange the last equality. Since the first $n-1$ rows of Z are orthogonal to
$\textbf{z}^{\top}$,
$\textbf{Z}\textbf{z}=\left\|\textbf{z}\right\|\textbf{i}_{n}$ by equation
(3.1), and consequently equation (4.5) implies that

$\left(\textbf{c}_{Z}\right)_{n}=\frac{\rho_{\textrm{e}}}{F\left\|\textbf{z}\right\|},$
(4.6)

i.e., the last entry of $\textbf{c}_{Z}$ quantifies excess charge density.
Thus $\textbf{c}_{Z}$ naturally partitions into a column that represents an
ordered set of $n-1$ neutral component concentrations, written as
$\textbf{c}_{\nu}$, followed by an entry expressing the excess charge in molar
units,

$\textbf{c}_{Z}=\left[\begin{array}[]{c}\textbf{c}_{\nu}\\\
\dfrac{\rho_{\textrm{e}}}{F\left\|\textbf{z}\right\|}\end{array}\right],$
(4.7)

which is zero when the local electroneutrality approximation holds. Entry
$c_{\nu,i}=\left(\textbf{c}_{\nu}\right)_{i}$ of the component concentrations
$\textbf{c}_{\nu}$ represents the molarity of the product formed by the $i$th
simple association equilibrium.

Natural flux descriptors over a salt–charge basis can be created by a similar
route, by demanding that the structure of the energy dissipation $T\dot{s}$ is
preserved. Assemble the species fluxes into a column matrix
$\vec{\textbf{n}}=[\vec{N}_{1},...,\vec{N}_{n}]^{\top}$.111The overarrow
indicates that each entry in the column $\vec{\textbf{n}}$ is a vector
quantity. Whenever the matrix multiplication makes sense,
$\textbf{A}\vec{\textbf{B}}$ forms a matrix $\vec{\textbf{C}}$ with vector
entries $\vec{C}_{ij}=\sum_{k}A_{ik}\vec{B}_{kj}$, and
$\vec{\textbf{A}}\cdot\vec{\textbf{B}}$ forms C, with scalar entries
$C_{ij}=\sum_{k}\vec{A}_{ik}\cdot\vec{B}_{kj}$. Structure preservation of
equation (4.2) requires that

$T\dot{s}=\vec{\textbf{n}}^{\top}\cdot\left(-\nabla\bm{\mu}\right)=\vec{\textbf{n}}_{Z}^{\top}\cdot\left(-\nabla\bm{\mu}_{Z}\right),$
(4.8)

where the column $\vec{\textbf{n}}_{Z}$ lists fluxes over the chosen
salt–charge basis Z. (The notation $\nabla\textbf{v}$ indicates a column
matrix whose $i$th entry is the vector $\nabla v_{i}$, where $v_{i}$ is the
$i$th entry of v.) The transformation from $\vec{\textbf{n}}$ into
$\vec{\textbf{n}}_{Z}$ therefore must be

$\vec{\textbf{n}}_{Z}=\textbf{Z}^{-\top}\vec{\textbf{n}},$ (4.9)

identical to the mapping that sends species concentrations into component
concentrations.

Faraday’s law—this time, for current—also provides detail about the new fluxes
$\vec{\textbf{n}}_{Z}$. In matrix form, Faraday’s law casts the current
density $\vec{i}$ in terms of species fluxes as

$\vec{i}=F\textbf{z}^{\top}\vec{\textbf{n}}=F\left(\textbf{Z}\textbf{z}\right)^{\top}\vec{\textbf{n}}_{Z}.$
(4.10)

Again, the fact that
$\textbf{Z}\textbf{z}=\left\|\textbf{z}\right\|\textbf{i}_{n}$ comes into
play, yielding

$\left(\vec{\textbf{n}}_{Z}\right)_{n}=\frac{\vec{i}}{F\left\|\textbf{z}\right\|}.$
(4.11)

Thus the last entry of $\vec{\textbf{n}}_{Z}$ represents a renormalized
current density. It is again natural to partition $\vec{\textbf{n}}_{Z}$, as

$\vec{\textbf{n}}_{Z}=\left[\begin{array}[]{c}\vec{\textbf{n}}_{\nu}\\\
\dfrac{\vec{i}}{F\left\|\textbf{z}\right\|}\end{array}\right],$ (4.12)

in which $\vec{\textbf{n}}_{\nu}$ is an $(n-1)$-dimensional column matrix
representing the set of _total component fluxes_.

Other fundamental thermodynamic laws also readily transform. For example,
prior definitions imply that

$\sum_{i=1}^{n}c_{i}\nabla\mu_{i}=\textbf{c}^{\top}\nabla\bm{\mu}=\textbf{c}_{Z}^{\top}\nabla\bm{\mu}_{Z}=\vec{0},$
(4.13)

so the change to a salt–charge basis retains a structurally identical
isothermal, isobaric Gibbs–Duhem equation.

One can also cast the total species concentration $c_{\textrm{T}}$ in terms of
the component concentrations. Letting 1 stand for a column matrix whose
entries are all ‘1’, one has that

$c_{\textrm{T}}=\sum_{i=1}^{n}c_{i}=\textbf{1}^{\top}\textbf{c}=\bm{\nu}_{Z}^{\top}\textbf{c}_{Z},$
(4.14)

in which the column matrix $\bm{\nu}_{Z}$, defined as

$\bm{\nu}_{Z}=\textbf{Z}\textbf{1}=\left[\begin{array}[]{c}\bm{\nu}\\\
\dfrac{\textbf{z}^{\top}\textbf{1}}{\left\|\textbf{z}\right\|}\end{array}\right],$
(4.15)

summarizes information about net reactant stoichiometry in the simple
association equilibria. The second equality here indicates that the
stoichiometric column $\bm{\nu}_{Z}$ also partitions into $(n-1)$ entries
$\bm{\nu}$ that describe components and a final entry associated with charge.
When component $i<n$ derives from a naturally uncharged species, the
corresponding entry of $\bm{\nu}$, that is, $\nu_{i}$, is 1; when component
$i<n$ is a salt, $\nu_{i}$ is the total number of ions in its formula unit.

Molar concentrations are convenient for stating the thermodynamic laws, but
they make characterization efforts error-prone because molarity varies with
temperature and pressure at fixed composition. It is useful to introduce two
alternative composition descriptors, called the _species fractions_ y and
_component fractions_ $\textbf{y}_{Z}$, and defined as

$\textbf{y}=\frac{1}{c_{\textrm{T}}}\textbf{c}~{}~{}~{}~{}\textrm{and}~{}~{}~{}~{}\textbf{y}_{Z}=\frac{1}{c_{\textrm{T}}}\textbf{c}_{Z},$
(4.16)

respectively. These afford the fundamental properties that

$\textbf{1}^{\top}\textbf{y}=1~{}~{}~{}~{}\textrm{and}~{}~{}~{}~{}\bm{\nu}_{Z}^{\top}\textbf{y}_{Z}=1,$
(4.17)

and are independent of temperature or pressure.

##  5 Governing equations over a salt–charge basis

The move to a salt–charge basis preserves the structure of material balances.
In the absence of homogenous reactions, all the species continuity equations
are summarized in matrix form as

$\frac{\partial\textbf{c}}{\partial t}=-\nabla\cdot\vec{\textbf{n}}.$ (5.1)

(Given any column $\vec{\textbf{v}}$, the operation
$\nabla\cdot\vec{\textbf{v}}$ forms a column matrix whose $i$th entry is the
scalar $\nabla\cdot\vec{v}_{i}$.) Multiplication by $\textbf{Z}^{-\top}$ and
use of equations (4.4) and (4.9) yields

$\frac{\partial\textbf{c}_{Z}}{\partial t}=-\nabla\cdot\vec{\textbf{n}}_{Z},$
(5.2)

confirming structure preservation. Bringing in the decompositions from
equations (4.7) and (4.12), one can alternatively express equation (5.2) as a
set of $n-1$ component mole balances,

$\frac{\partial\textbf{c}_{\nu}}{\partial
t}=-\nabla\cdot\vec{\textbf{n}}_{\nu},$ (5.3)

along with the single equation

$\frac{\partial\rho_{\textrm{e}}}{\partial t}=-\nabla\cdot\vec{i},$ (5.4)

which establishes charge continuity. Note that any model governed by species
material balances which also incorporates Faraday’s law is compatible with
Maxwellian electrodynamics, because Maxwell’s equations imply a charge-
continuity relationship identical to equation (5.4).

Because the transformations that form $\textbf{c}_{Z}$ and
$\vec{\textbf{n}}_{Z}$ preserve the structure of the dissipation function,
they also preserve the structure of the Onsager–Stefan–Maxwell constitutive
laws. Equation (1.1) takes the matrix form

$-\nabla\bm{\mu}=\textbf{M}\vec{\textbf{n}}.$ (5.5)

Equations (3.2) and (4.9) show that this transforms to

$-\nabla\bm{\mu}_{Z}=\textbf{M}_{Z}\vec{\textbf{n}}_{Z}$ (5.6)

over a salt–charge basis. The transport matrix $\textbf{M}_{Z}$ in the
salt–charge representation of the Onsager–Stefan–Maxwell laws relates to the
matrix M from equation (1.2) through

$\textbf{M}_{Z}=\textbf{Z}\textbf{M}\textbf{Z}^{\top},$ (5.7)

a congruence transformation. It follows from this congruence that M is
symmetric if and only if $\textbf{M}_{Z}$ is symmetric, so that the
Onsager–Stefan–Maxwell transport coefficients retain symmetry when expressed
over a salt–charge basis. Sylvester’s law of inertia implies further that the
signatures and ranks of M and $\textbf{M}_{Z}$ equate. As mentioned before,
the non-dissipative nature of convection implies that M affords c as its sole
null eigenvector; in light of congruence relation (5.7) and Sylvester’s law,
this implies that $\textbf{c}_{Z}$ is the sole null eigenvector of
$\textbf{M}_{Z}$.

Thus, transformation to a salt–charge basis preserves the structure of the
differential balance equations as well as all spectral properties of the
Onsager–Stefan–Maxwell equations. Both M and $\textbf{M}_{Z}$ are symmetric
matrices, and each affords a single null eigenvector. Although the null spaces
of M and $\textbf{M}_{Z}$ generally differ, each matrix affords just one null
eigenvector (c and $\textbf{c}_{Z}$, respectively). The nonzero eigenvalues of
M generally differ in magnitude from those of $\textbf{M}_{Z}$, but both sets
are all positive.

##  6 The salt–charge potential

It is prudent to introduce a thermodynamic potential with units of voltage to
quantify electrical energy and power. To put electricity in a thermodynamic
context, however, one must acknowledge an essential limitation on how
electrochemical equilibria are parametrized.

When proposing the electrochemical potential concept, Guggenheim introduced
constitutive laws

$\mu_{i}=RT\ln a_{i}+Fz_{i}\Phi,$ (6.1)

in which $a_{i}$ represents the chemical activity of species $i$ and $\Phi$ is
an electric potential in the solution phase [14, 17, 1]. Observing that the
equilibrium compositions of practical macroscopic systems cannot be adjusted
along paths that violate global electroneutrality, Guggenheim noted that the
partitioning of electrochemical potentials into chemical and electrical parts
is necessarily ambiguous. Therefore a model of electrolyte energetics must
never depend on the chemical activity of a single charged species in isolation
[17]. To cast this notion formally, define a non-neutral free-energy
contribution $\mu_{z}$, written in terms of Guggenheim’s species activities
and solution potential as

$\mu_{z}=\frac{1}{\left\|\textbf{z}\right\|}\sum\limits_{i=1}^{n}z_{i}\mu_{i}=\frac{RT}{\left\|\textbf{z}\right\|}\sum\limits_{i=1}^{n}z_{i}\ln
a_{i}+F\left\|\textbf{z}\right\|\Phi.$ (6.2)

Guggenheim’s principle states that no macroscopic experiment can discern the
term that involves the solution potential $\Phi$ in this expression from those
involving the species activities $a_{i}$ [14].

Rather than using the solution potential $\Phi$, we instead employ a quantity
we call the _salt–charge potential_ $\Phi_{z}$, defined as

$\Phi_{z}=\frac{\mu_{z}}{F\left\|\textbf{z}\right\|}.$ (6.3)

Guggenheim’s principle requires that no experiment can distinguish chemical
and electrical contributions to $\Phi_{z}$. Despite this ambiguity, one can
still understand the salt–charge potential $\Phi_{z}$ as a type of electrical
potential. Inserting decompositions (3.3) and (4.7), then introducing
definition (6.3) recasts the Euler equation for $\tilde{G}$ as

$\tilde{G}=\bm{\mu}_{\nu}^{\top}\textbf{c}_{\nu}+\Phi_{z}\rho_{\textrm{e}}.$
(6.4)

Here, the $\Phi_{z}\rho_{\textrm{e}}$ term apparently quantifies excess
coulomb energy, that is, the available energy associated with imbalanced
charge. Similarly, inserting decompositions (3.3) and (4.12) restates the
dissipation as

$T\dot{s}=\vec{\textbf{n}}_{\nu}^{\top}\cdot\left(-\nabla\bm{\mu}_{\nu}\right)+\vec{i}\cdot\left(-\nabla\Phi_{z}\right).$
(6.5)

The $-\vec{i}\cdot\nabla\Phi_{z}$ term is the local ohmic loss—the rate of
energy dissipation associated with electrical current flow. Observe that local
electroneutrality makes the coulomb term disappear from $\tilde{G}$, but
leaves the ohmic term in $T\dot{s}$.

Although the species material balances summarized in equation (5.1) imply
charge-continuity equation (5.4), and are therefore compatible with Maxwell’s
equations, equilibrium thermodynamics cannot provide insight into the
Maxwellian electric field, because Guggenheim’s principle mandates that it is
impossible to determine experimentally whether or not $-\nabla\Phi_{z}$ is a
purely electrical quantity.

Still, the energy density and dissipative power associated with salt–charge
potential (and its gradient) have unambiguous meanings, because the coulomb
term in $\tilde{G}$ and ohmic term in $T\dot{s}$ do not vary with the choice
of components. Whereas one can only speak of ‘a’ salt–charge basis since Z is
not always unique, one may refer to ‘the’ salt–charge potential, which is
unique for a given z.

##  7 Thermodynamic factors in a salt–charge basis

To implement a simulation of any transport system, one requires constitutive
relationships that express how thermodynamic potentials depend on the
concentrations of extensive properties. The particular mapping that sends
composition gradients into electrochemical-potential gradients in isothermal,
isobaric situations is expressed by a matrix of equilibrium properties
commonly referred to as _thermodynamic factors_.

For an isothermal, isobaric, $n$-species mass-transport system, the
thermodynamic factors depend parametrically on $\tfrac{1}{2}n\left(n-1\right)$
composition-dependent state functions—one state function for every pair of
species [6, 16, 18]. In their discussion of Onsager reciprocal relations for
multicomponent diffusion, Monroe and Newman introduced an
$\left(n-1\right)\times\left(n-1\right)$ matrix of dimensionless Darken
factors Q [19], defined such that

$\operatorname{diag}\left(\textbf{y}\right)\nabla\bm{\mu}=RT\left[\begin{array}[]{c}\textbf{I}_{n-1}\\\
-\textbf{1}^{\top}\end{array}\right]\textbf{Q}\left[\begin{array}[]{cc}\textbf{I}_{n-1}&\textbf{o}\end{array}\right]^{\top}\nabla\textbf{y},$
(7.1)

where $\textbf{I}_{n-1}$ indicates the
$\left(n-1\right)\times\left(n-1\right)$ identity matrix, o represents a
column of zeroes, and the linear operator
$\operatorname{diag}\left(\textbf{v}\right)$ forms column matrix v into a
diagonal square matrix whose $i$th diagonal entry is $v_{i}$.222 Generally
$\operatorname{diag}\left(\textbf{a}\right)\textbf{b}=\operatorname{diag}\left(\textbf{b}\right)\textbf{a}$;
both products form a column c with entries $c_{i}=a_{i}b_{i}$. Also
$\textbf{1}^{\top}\operatorname{diag}(\textbf{a})=\textbf{a}^{\top}$ and
$\operatorname{diag}(\textbf{a})\textbf{1}=\textbf{a}$. The notation
$\operatorname{diag}\left(\textbf{a}\right)^{-1}$ indicates the diagonal
matrix whose $i$th entry is $1/a_{i}$, such that
$\operatorname{diag}\left(\textbf{a}\right)^{-1}\textbf{a}=\textbf{1}$.  Each
entry $Q_{ij}$ of Q represents the thermodynamic derivative of the (log)
activity of species $i$ with respect to the particle fraction of species
$j<n$, leaving temperature, pressure, and all particle fractions save those of
species $j$ and $n$ fixed. (Experiments to measure $Q_{ij}$ adopt a convention
that when one varies $y_{j}$, one specifically adjusts $y_{n}$ to maintain the
constraint that $\textbf{1}^{\top}\textbf{y}=1$.) In terms of Guggenheim’s
species activities, $Q_{ij}$ is defined as

$\displaystyle Q_{ij}=y_{i}\left(\frac{\partial\ln a_{i}}{\partial
y_{j}}\right)_{T,p,y_{k\neq
j,n}}=\delta_{ij}+y_{i}\left(\frac{\partial\ln\lambda_{i}}{\partial
y_{j}}\right)_{T,p,y_{k\neq j,n}},$ (7.2)

in which $\delta_{ij}$ represents the Kronecker delta. The expression on the
right of equation (7.2) introduces the _activity coefficient_ of species $i$,
$\lambda_{i}$, defined such that $a_{i}=\lambda_{i}y_{i}$. Introducing
$\lambda_{i}$ emphasizes the expectation that $\textbf{Q}=\textbf{I}_{n-1}$
within solutions whose mixing free energies are thermodynamically ideal.

Although the simplicity of Q in the ideal case is useful, Maxwell relations
among these Darken factors in non-ideal situations are obscure. It is
sometimes convenient to work instead with the $(n-1)\times(n-1)$ Hessian
matrix K,

$K_{ij}=\frac{\partial^{2}}{\partial y_{i}\partial
y_{j}}\left(\frac{\tilde{G}}{RTc_{\textrm{T}}}\right)_{T,p,y_{k\neq i,j,n}},$
(7.3)

for which Maxwell relations imply that $\textbf{K}=\textbf{K}^{\top}$ [6].
Thermodynamic stability demands that K is positive semidefinite [16]. The
Darken factors Q depend on K through

$\textbf{Q}=\textbf{Y}^{-1}\textbf{K},$ (7.4)

where the definition

$\textbf{Y}^{-1}=\left[\begin{array}[]{cc}\textbf{I}_{n-1}&\textbf{o}\end{array}\right]\left(\operatorname{diag}\left(\textbf{y}\right)-\textbf{y}\textbf{y}^{\top}\right)\left[\begin{array}[]{c}\textbf{I}_{n-1}\\\
\textbf{o}^{\top}\end{array}\right]$ (7.5)

puts in matrix form the $(n-1)\times(n-1)$ inverse composition matrix
$\textbf{Y}^{-1}$ employed by Monroe and Newman [20]. Substituting equations
(7.4) and (7.5) into equation (7.1) and algebraically simplifying the result
gives

$\nabla\bm{\mu}=RT\left(\textbf{I}-\textbf{1}\textbf{y}^{\top}\right)\left[\begin{array}[]{cc}\textbf{K}&\textbf{o}\\\
\textbf{o}^{\top}&0\end{array}\right]\nabla\textbf{y},$ (7.6)

which shows how the isothermal, isobaric composition dependence of
$\nabla\bm{\mu}$ is parametrized by the $\tfrac{1}{2}n\left(n-1\right)$
independent entries of K. Observe that the block matrix here has a nullspace
spanned by the null eigenvector $\textbf{i}_{n}$, which reflects the fact that
gradients $\nabla y_{n}$ have no independent effect on species electrochemical
potentials.

Some of the benefits of the Darken matrix Q can be brought into the analysis
of equations expressed in terms of K. The composition Hessian separates into
two parts,

$\textbf{K}=\textbf{Y}+\Delta\textbf{K},$ (7.7)

such that Y accounts for derivatives of ideal mixing free energy and
$\Delta\textbf{K}$ expresses the Hessian of excess free energy. The
composition matrix Y here is written explicitly as

$\textbf{Y}=\left[\begin{array}[]{cc}\textbf{I}_{n-1}&\textbf{o}\end{array}\right]\left(\operatorname{diag}\left(\textbf{y}\right)^{-1}+\frac{\textbf{1}\textbf{1}^{\top}}{\textbf{i}_{n}^{\top}\textbf{y}}\right)\left[\begin{array}[]{c}\textbf{I}_{n-1}\\\
\textbf{o}^{\top}\end{array}\right].$ (7.8)

Upon substitution into equation (7.6), algebraic simplification yields an
equivalent alternative form,

$\displaystyle\nabla\bm{\mu}=\,$ $\displaystyle
RT\operatorname{diag}\left(\textbf{y}\right)^{-1}\nabla\textbf{y}$
$\displaystyle+RT\left(\textbf{I}-\textbf{1}\textbf{y}^{\top}\right)\left[\begin{array}[]{cc}\Delta\textbf{K}&\textbf{o}\\\
\textbf{o}^{\top}&0\end{array}\right]\nabla\textbf{y},$ (7.11)

in which the $(n-1)\times(n-1)$ matrix $\Delta\textbf{K}$ is symmetric, and
vanishes when mixing is ideal.

The nonideal term in equation (7) translates readily into a salt–charge basis.
Use the congruence relationship

$\Delta\textbf{K}_{Z}=\textbf{Z}\left[\begin{array}[]{cc}\Delta\textbf{K}&\textbf{o}\\\
\textbf{o}^{\top}&0\end{array}\right]\textbf{Z}^{\top}$ (7.12)

to define a transformed excess Hessian $\Delta\textbf{K}_{Z}$. All $n^{2}$
entries of this symmetric matrix are generally nonzero. The nullspace of
$\Delta\textbf{K}_{Z}$ has a minimum dimension of $1$, however, because
$\textbf{Z}^{-\top}\textbf{i}_{n}$ is always a null eigenvector. Multiplying
equation (7) through by Z and inserting $\Delta\textbf{K}_{Z}$, one finds

$\displaystyle\nabla\bm{\mu}_{Z}=$ $\displaystyle
RT\textbf{Z}\operatorname{diag}\left(\textbf{y}\right)^{-1}\nabla\textbf{y}$
$\displaystyle+RT\left(\textbf{I}-\bm{\nu}_{Z}\textbf{y}_{Z}^{\top}\right)\Delta\textbf{K}_{Z}\nabla\textbf{y}_{Z}$
(7.13)

after using equations (3.2), (4.4), and (4.17) to simplify. The assertion that
species potentials within an electrolyte at constant temperature and pressure
depend on all the independent species fractions necessitates that the
differential component potentials have this structure.

It is noteworthy that the ideal part of the component chemical-potential
gradients—the first term on the right of equation (7)—involves entries of the
salt–charge transformation Z as prefactors. Thus the ideal mixing free energy
of a dissolved salt scales linearly with its ion stoichiometry, in addition to
the ideal logarithmic dependence on its mole fraction. Note also that
multiplication from the left by
$\textbf{c}_{Z}^{\top}=c_{\textrm{T}}\textbf{y}_{Z}^{\top}$ verifies that the
Gibbs–Duhem equation is satisfied by constitutive laws (7) for any choice of
the symmetric excess Hessian matrix $\Delta\textbf{K}_{Z}$.

Although equation (7) is formally correct, the parameters within
$\Delta\textbf{K}_{Z}$ are not all measurable. In particular, as discussed
after equation (6.2), adherence to Guggenheim’s principle demands that the
coefficients of $\nabla\textbf{y}_{Z}$ on the right side of the equation for
$\nabla\mu_{z}$ should remain ambiguous because they cannot be disentangled
with equilibrium experiments. Therefore it is appropriate to discard this
component of the expansion in equation (7).

The notation is simplified by defining an $n\times(n-1)$ matrix N, whose
columns represent the reactant stoichiometry in the simple association
equilibria:

$\textbf{N}=\left[\begin{array}[]{ccc}\bm{\nu}_{1}&\dots&\bm{\nu}_{n-1}\end{array}\right],~{}~{}~{}\textrm{so}~{}~{}~{}\textbf{Z}=\left[\begin{array}[]{c}\textbf{N}^{\top}\\\
\dfrac{\textbf{z}^{\top}}{\left\|\textbf{z}\right\|}\end{array}\right].$
(7.14)

Further observing that $\bm{\nu}=\textbf{N}^{\top}\textbf{1}$, the general
expression

$\displaystyle\nabla\bm{\mu}_{\nu}=$ $\displaystyle
RT\textbf{N}^{\top}\operatorname{diag}\left(\textbf{y}\right)^{-1}\nabla\textbf{y}$
$\displaystyle+RT\left(\left[\begin{array}[]{cc}\textbf{I}_{n-1}&\textbf{o}\end{array}\right]-\bm{\nu}\textbf{y}_{Z}^{\top}\right)\Delta\textbf{K}_{Z}\nabla\textbf{y}_{Z}$
(7.16)

expresses the portion of the expansion that remains after striking the last
row of matrix equation (7) to accommodate Guggenheim’s principle. In the next
section, this form of the component chemical-potential gradients will
facilitate analyzing constraints on properties under the electroneutrality
approximation.

##  8 Structural implications of electroneutrality

Imposing local electroneutrality on the Onsager–Stefan–Maxwell transport model
impacts both the variables involved and the parameters that quantify material
properties. The approximation separately impacts the thermodynamic and
dynamical aspects of the general theory.

###  8.1 Electroneutral composition

The simplest consequence of electroneutrality is that it limits the available
composition space, because one charged-species concentration becomes linearly
dependent on the others. It is useful to establish some relationships that
demonstrate the effect of this constraint.

When excess charge density vanishes everywhere, the $n$-dimensional column of
species concentrations c becomes limited to a set of locally electroneutral
species concentrations, which we notate as $\textbf{c}^{0}$. Necessarily, any
solution compositions in the space that $\textbf{c}^{0}$ occupies must be
instantiable by assembling uncharged components in the appropriate proportions
and then dissociating them into their constituent species. This idea can be
put in mathematical terms by leveraging a salt–charge basis. Over such a basis
$\textbf{c}^{0}$ is determined by the mapping

$\textbf{c}^{0}=\textbf{Z}^{\top}\left[\begin{array}[]{c}\textbf{c}_{\nu}\\\
0\end{array}\right]=\textbf{N}\textbf{c}_{\nu},$ (8.1)

which results from equations (4.4) and (4.7) when $\rho_{\textrm{e}}=0$. Thus
the first $(n-1)$ columns of $\textbf{Z}^{\top}$ (that is, the columns of N)
serve as a basis for the subspace occupied by $\textbf{c}^{0}$, and the
$(n-1)$-dimensional column of electroneutral component concentrations
$\textbf{c}_{\nu}$ expresses the coordinates of a given $\textbf{c}^{0}$ over
this basis. Electroneutrality constrains the total concentration’s domain in a
similar way. The electroneutral total concentration, $c_{\textrm{T}}^{0}$,
depends on $\textbf{c}_{\nu}$ as

$c_{\textrm{T}}^{0}=\textbf{1}^{\top}\textbf{c}^{0}=\bm{\nu}_{Z}^{\top}\left[\begin{array}[]{c}\textbf{c}_{\nu}\\\
0\end{array}\right]=\bm{\nu}^{\top}\textbf{c}_{\nu},$ (8.2)

simplified with the definition of $\bm{\nu}$ from equation (4.15).

Some additional physical principles outside thermodynamics restrict the range
of neutral compositions. Since negative species molarities are meaningless,
every entry of $\textbf{c}^{0}$ must be non-negative. Because mass transport
cannot occur in the absence of species, $c_{\textrm{T}}^{0}$ is strictly
positive. These inequalities place bounds on the entries of the component-
composition column $\textbf{c}_{\nu}$.

The restrictions on $\textbf{c}^{0}$ and $c_{\textrm{T}}^{0}$ also imply that
every entry of the electroneutral species fractions, $\textbf{y}^{0}$, defined
as

$\textbf{y}^{0}=\frac{1}{c_{\textrm{T}}^{0}}\textbf{c}^{0}=\textbf{Z}^{\top}\left[\begin{array}[]{c}\textbf{y}_{\nu}\\\
0\end{array}\right]=\textbf{N}\textbf{y}_{\nu},$ (8.3)

must be strictly non-negative. Under electroneutrality this object retains the
property that $\textbf{1}^{\top}\textbf{y}^{0}=1$, established by equation
(4.17). Furthermore, the condition

$\bm{\nu}^{\top}\textbf{y}_{\nu}=1~{}~{}~{}~{}~{}~{}~{}~{}\textrm{(electroneutral)}$
(8.4)

constrains component fractions under electroneutrality.

Note that requiring $\textbf{c}^{0}$ or $\textbf{y}^{0}$ to have non-negative
entries does not ensure that every component concentration within
$\textbf{c}_{\nu}$ (or $\textbf{y}_{\nu}$) is non-negative. Over a salt–charge
basis, negative component concentrations can occur when the choice of products
in the simple association equilibria is not unique. Consider the solution of
$\ch{Na+}$, $\ch{Cl-}$, $\ch{Mg}^{2+}$ and $\ch{SO_{4}^{2-}}$ in $\ch{H2O}$
from section 2 by way of example. The electroneutral composition space for
this solution reaches extremes corresponding to four distinct binary
electrolytes: aqueous $\ch{NaCl}$ with a trace amount of $\ch{MgSO4}$; aqueous
$\ch{MgCl2}$ with trace $\ch{Na2SO4}$; aqueous $\ch{Na2SO4}$ with trace
$\ch{MgCl2}$; and aqueous $\ch{MgSO4}$ with trace $\ch{NaCl}$. Suppose the
neutral salts for a salt–charge basis are chosen to be $\ch{NaCl}$,
$\ch{MgCl2}$, and $\ch{Na2SO4}$, as in section 2. In the laboratory, these
salts cannot be mixed with water to produce, say, 1 M aqueous $\ch{MgSO4}$. To
formulate 1 M $\ch{MgSO4}$ with the given precursor salts, one would have to
add $\ch{Na2SO4}$ to supply the sulfate, add $\ch{MgCl2}$ to introduce
magnesium, and then utilize a separation process (e.g., precipitation,
heterogeneous extraction, etc.) to remove $\ch{NaCl}$. Thus, aqueous 1 M
$\ch{MgSO4}$ is represented by a $\textbf{c}_{\nu}$ column with entries for 1
M $\ch{Na2SO4}$, 1 M $\ch{MgCl2}$, and $-2$ M $\ch{NaCl}$. Although the
$\ch{NaCl}$ molarity is negative, that of every ion remains non-negative.
Since $\ch{MgSO4}$ can be formed by recombination reaction (2.1), its chemical
potential depends unambiguously on those of $\ch{Na2SO4}$, $\ch{MgCl2}$, and
$\ch{NaCl}$. The negative concentration just indicates that $\ch{NaCl}$ is a
product of a recombination reaction involving the salt–charge basis’
components, rather than a reactant.

###  8.2 Component activity coefficients

Electroneutrality has ramifications for Guggenheim’s constitutive framework.
Definition (6.1) can be rewritten in terms of species activity coefficients
$\lambda_{i}^{0}$ as

$\displaystyle\mu_{i}^{0}$
$\displaystyle=\mu_{i}^{\ominus}+RT\ln\left(\lambda_{i}^{0}y_{i}^{0}\right)+z_{i}F\Phi,~{}~{}~{}~{}\textrm{or}$
$\displaystyle\bm{\mu}^{0}$
$\displaystyle=\bm{\mu}^{\ominus}+RT\ln\left(\textbf{y}^{0}\right)+RT\ln\left(\bm{\lambda}^{0}\right)+F\Phi\textbf{z},$
(8.5)

where the composition-independent parameters $\mu_{i}^{\ominus}$ that make up
the $n$-dimensional column $\bm{\mu}^{\ominus}$ quantify species
electrochemical potentials in a secondary reference state, and
$\ln\left(\textbf{v}\right)$ represents a column matrix wherein each entry is
the natural logarithm of the corresponding entry of v. Thermodynamic
consistency mandates that every activity coefficient comprising
$\bm{\lambda}^{0}$ is positive. Since $(\textbf{c}^{0})^{\top}\textbf{z}=0$ by
design, these constitutive laws meet the requirement that

$\tilde{G}=\left(\textbf{c}^{0}\right)^{\top}\bm{\mu}^{0}=\textbf{c}_{\nu}^{\top}\bm{\mu}_{\nu}^{0}$
(8.6)

as a consequence of equations (3.2) and (8.1).

Constitutive laws for component chemical potentials in accord with
Guggenheim’s formulation are identified by using the second equality in
equation (8.6), which implies

$\left(\textbf{c}^{0}\right)^{\top}\textbf{Z}^{-1}\left(\textbf{Z}\bm{\mu}^{0}-\bm{\mu}_{Z}^{0}\right)=0.$
(8.7)

Insertion of equations (7.14) and (8.2) into equation (8.7), followed by
partitioning of $\bm{\mu}_{Z}^{0}$ with equation (3.3), application of the
relation $\textbf{Z}\textbf{z}=\left\|\textbf{z}\right\|\textbf{i}_{n}$, and
substitution of $\Phi$ in favor of $\Phi_{z}^{0}$ with equations (6.2) and
(6.3), reveals that Guggenheim’s laws are consistent with component chemical
potentials constituted by equations of the form

$\mu_{\nu,k}^{0}=\mu_{\nu,k}^{\ominus}+RT\sum\limits_{m=1}^{n}\nu_{km}\ln
y_{m}^{0}+RT\nu_{k}\ln\lambda_{\nu,k}^{0},$ (8.8)

wherein $\mu_{\nu,k}^{\ominus}\left(T,p\right)$ is the chemical potential of
component $k$ in the reference state and $\nu_{k}$ is the $k$th entry of the
column $\bm{\nu}$ that quantifies total component stoichiometry. The parameter
$\lambda_{\nu,k}^{0}$ is called the _mean molar activity coefficient_ of
component $k$.

All of the component chemical-potential constitutive laws can therefore be
summarized by a matrix equation,

$\bm{\mu}_{\nu}^{0}=\bm{\mu}_{\nu}^{\ominus}+\bm{\mu}_{\nu}^{0,\textrm{ideal}}+RT\operatorname{diag}\left(\bm{\nu}\right)\ln\bm{\lambda}_{\nu}^{0},$
(8.9)

in which the term

$\bm{\mu}_{\nu}^{0,\textrm{ideal}}=RT\textbf{N}^{\top}\ln\textbf{y}^{0}$
(8.10)

accounts for the part of a component’s chemical potential that is attributable
to ideal mixing of its constituent species, and the reference chemical
potentials and mean molar activity coefficients of each component make up the
$(n-1)$-dimensional columns $\bm{\mu}_{\nu}^{\ominus}$ and
$\bm{\lambda}_{\nu}^{0}$, respectively: the expression

$\bm{\mu}_{\nu}^{\ominus}=\textbf{N}^{\top}\bm{\mu}^{\ominus},~{}~{}~{}\textrm{or}~{}~{}~{}\mu_{\nu,k}^{\ominus}=\sum_{m=1}^{n}\nu_{km}\mu_{m}^{\ominus},$
(8.11)

where $k\in\left\\{1,...,n-1\right\\}$, shows how reactant stoichiometries in
the simple association equilibria relate the reference potentials of
components and species, and

$\displaystyle\ln\bm{\lambda}_{\nu}^{0}$
$\displaystyle=\operatorname{diag}\left(\bm{\nu}\right)^{-1}\textbf{N}^{\top}\ln\bm{\lambda}^{0},~{}~{}~{}\textrm{or}$
$\displaystyle\ln\lambda_{\nu,k}^{0}$
$\displaystyle=\sum_{m=1}^{n}\frac{\nu_{km}}{\nu_{k}}\ln\lambda_{m}^{0},$
(8.12)

puts the mean molar component activity coefficients in terms of species
activity coefficients. Note that the versions of equations (8.11) and (8.2)
written in terms of summations identify with the definitions stated by both
Guggenheim [17] and Newman [1].

###  8.3 Electroneutral thermodynamic factors

Applying electroneutrality to equations (4.13) and (6.4) reveals equilibrium
energetics to be independent of $\Phi_{z}^{0}$. The electroneutral Euler and
Gibbs–Duhem equations,

$\tilde{G}=\left(\bm{\mu}_{\nu}^{0}\right)^{\top}\textbf{c}_{\nu}$ (8.13)

and

$\textbf{c}_{\nu}^{\top}\nabla\bm{\mu}_{\nu}^{0}=\vec{0},$ (8.14)

respectively, imply that isothermal, isobaric component chemical potentials
have the functionality
$\bm{\mu}_{\nu}^{0}=\bm{\mu}_{\nu}^{0}\left(\textbf{c}_{\nu}\right)$, or,
equivalently, $\bm{\mu}_{\nu}^{0}\left(\textbf{y}_{\nu}\right)$.333The
implication
$\bm{\mu}_{\nu}^{0}\left(\textbf{c}_{\nu}\right)\implies\bm{\mu}_{\nu}^{0}\left(\textbf{y}_{\nu}\right)$
relies on the fact that total concentration depends only on the species
particle fractions at fixed temperature and pressure. This is justified by
Euler equation (4.1) because the molar Gibbs energy
$\overline{G}=\tilde{G}/c_{\textrm{T}}$ must produce the molar volume
$1/c_{\textrm{T}}$ through the derivative $(\partial\overline{G}/\partial
p)_{T,y_{i}}=1/c_{\textrm{T}}$. The mean molar activity coefficients in
electroneutral constitutive laws (8.9) thus depend only on $\textbf{y}_{\nu}$
at constant temperature and pressure.

These facts come together with the discussion of Darken factors from section 7
to establish Maxwell relations and component thermodynamic factors for the
electroneutral Onsager–Stefan–Maxwell diffusion driving forces. Let

$\Lambda_{\nu,ij}^{0}=\nu_{i}\left(\frac{\partial\ln\lambda_{\nu,i}^{0}}{\partial
y_{\nu,j}}\right)_{T,p,y_{\nu,k\neq j,n-1}}$ (8.15)

represent the partial derivative of the mean molar activity coefficient of
component $i$ with respect to the fraction of component $j$, leaving all of
the $n-1$ component fractions save the $j$th and $(n-1)$th fixed. In light of
Gibbs–Duhem equation (8.14), the component chemical-potential gradients can be
expanded as

$\displaystyle\nabla\bm{\mu}_{\nu}^{0}=$
$\displaystyle\nabla\bm{\mu}_{\nu}^{0,\textrm{ideal}}$
$\displaystyle+RT\left(\textbf{I}_{n-1}-\bm{\nu}\textbf{y}_{\nu}^{\top}\right)\left[\begin{array}[]{cc}\bm{\Lambda}_{\nu}^{0}&\textbf{o}\\\
\textbf{o}^{\top}&0\end{array}\right]\nabla\textbf{y}_{\nu},$ (8.18)

in which $\bm{\Lambda}_{\nu}^{0}$ is an $(n-2)\times(n-2)$ matrix of
electroneutral component Darken factors, and

$\displaystyle\nabla\bm{\mu}_{\nu}^{0,\textrm{ideal}}$
$\displaystyle=RT\textbf{N}^{\top}\operatorname{diag}\left(\textbf{y}^{0}\right)^{-1}\nabla\textbf{y}^{0}$
$\displaystyle=RT\textbf{N}^{\top}\operatorname{diag}\left(\textbf{N}\textbf{y}_{\nu}\right)^{-1}\textbf{N}\nabla\textbf{y}_{\nu}$
(8.19)

expresses the gradient of the ideal contribution to every component chemical
potential.

Equation (8.3) contains $(n-2)^{2}$ excess component Darken factors
$\Lambda^{0}_{\nu,ij}$, but these are not all independent because of Maxwell
relations, whose structure can be understood by revisiting the electroneutral
form of equation (7). After insertion of equation (8.3) and simplification
with equation (8.3), equation (7) becomes

$\nabla\bm{\mu}_{\nu}^{0}=\nabla\bm{\mu}_{\nu}^{0,\textrm{ideal}}+RT\left(\textbf{I}_{n-1}-\bm{\nu}\textbf{y}_{\nu}^{\top}\right)\Delta\textbf{K}_{\nu}^{0}\nabla\textbf{y}_{\nu},$
(8.20)

in which the $(n-1)$ by $(n-1)$ matrix $\Delta\textbf{K}_{\nu}^{0}$ is formed
by striking the $n$th row and $n$th column of $\Delta\textbf{K}_{Z}^{0}$,

$\Delta\textbf{K}_{\nu}^{0}=\left[\begin{array}[]{cc}\textbf{I}_{n-1}&\textbf{o}\end{array}\right]\Delta\textbf{K}_{Z}^{0}\left[\begin{array}[]{c}\textbf{I}_{n-1}\\\
\textbf{o}^{\top}\end{array}\right],$ (8.21)

and is consequently symmetric. Equating the right sides of equations (8.3) and
(8.20) shows that

$\left(\textbf{I}_{n-1}-\bm{\nu}\textbf{y}_{\nu}^{\top}\right)\left[\begin{array}[]{cc}\bm{\Lambda}_{\nu}^{0}&\textbf{o}\\\
\textbf{o}^{\top}&0\end{array}\right]=\left(\textbf{I}_{n-1}-\bm{\nu}\textbf{y}_{\nu}^{\top}\right)\Delta\textbf{K}_{\nu}^{0},$
(8.22)

which relates the thermodynamic factors that make up $\bm{\Lambda}_{\nu}^{0}$
to the truncated excess Hessian $\Delta\textbf{K}_{\nu}^{0}$. Compatibility of
the species and component perspectives thus demands that the last column of
$\Delta\textbf{K}_{\nu}^{0}$ lies in the nullspace of
$\left(\textbf{I}_{n-1}-\bm{\nu}\textbf{y}_{\nu}^{\top}\right)$, that is, it
must be proportional to $\bm{\nu}$. (Guggenheim’s principle demands that the
constant of proportionality cannot be determined, however.) The last row of
$\Delta\textbf{K}_{\nu}^{0}$ is therefore proportional to $\bm{\nu}^{\top}$,
as required by the Hessian’s symmetry. Multiplication through by
$\textbf{y}_{\nu}^{\top}$ additionally verifies that the last row of equation
(8.22) is linearly dependent, since Gibbs–Duhem relations constrain both sides
of the equality in the same way.

In short, the last row of equation (8.22) is redundant, and the last column,
trivial. Discarding these, isolating the truncated Hessian that remains by
exploiting a Sherman–Morrison formula, and further rearranging demonstrate
that for all $i\neq j$, the Maxwell relations

$\nu_{i}\left(\frac{\partial\ln\lambda_{\nu,i}^{0}}{\partial
y_{\nu,j}}\right)_{T,p,y_{\nu,k\neq
j,n-1}}=\nu_{j}\left(\frac{\partial\ln\lambda_{\nu,j}^{0}}{\partial
y_{\nu,i}}\right)_{T,p,y_{\nu,k\neq i,n-1}}$ (8.23)

hold; the matrix $\bm{\Lambda}^{0}_{\nu}$ is symmetric. This substantiates the
claim that $\tfrac{1}{2}\left(n-1\right)\left(n-2\right)$ properties, which
quantify isothermal, isobaric gradients of $n-1$ component chemical
potentials, parametrize the mixing free energy of any $n$-ary electroneutral
electrolyte.

Henceforth, we will write the electroneutral gradients of component chemical
potentials in the more compact form

$\nabla\bm{\mu}_{\nu}^{0}=RT\textbf{X}_{\nu}^{0}\nabla\textbf{y}_{\nu},$
(8.24)

in which the $(n-1)\times(n-1)$ matrix of _electroneutral thermodynamic
factors_ $\textbf{X}_{\nu}^{0}$ is defined as

$\displaystyle\textbf{X}_{\nu}^{0}=$
$\displaystyle\textbf{N}^{\top}\operatorname{diag}\left(\textbf{N}\textbf{y}_{\nu}\right)^{-1}\textbf{N}$
$\displaystyle+\left(\textbf{I}_{n-1}-\bm{\nu}\textbf{y}_{\nu}^{\top}\right)\left[\begin{array}[]{cc}\bm{\Lambda}_{\nu}^{0}&\textbf{o}\\\
\textbf{o}^{\top}&0\end{array}\right],$ (8.27)

wherein $\bm{\Lambda}^{0}_{\nu}=\left(\bm{\Lambda}^{0}_{\nu}\right)^{\top}$,
and equation (8.15) establishes how independent entries in
$\bm{\Lambda}_{\nu}^{0}$ can be measured experimentally.

###  8.4 Electroneutrality and dynamics

Application of the electroneutrality approximation to the dynamical governing
equations is simpler than electroneutral thermodynamics. When excess charge
density vanishes, there is no change to the $(n-1)$-dimensional component
material balance equation (5.3) that was produced by moving to a salt–charge
basis. Charge continuity, embodied by equation (5.4), does simplify slightly,
to

$\nabla\cdot\vec{i}=0,~{}~{}~{}~{}~{}~{}~{}~{}\textrm{(electroneutral)}$
(8.28)

which expresses Kirchhoff’s law of the node.

Entries in the Onsager drag matrix from equation (1.2) change their structure
under electroneutrality. The transport matrix $\textbf{M}^{0}$ that relates
species fluxes to species driving forces in the electroneutral case has
entries

$M_{ij}^{0}=\begin{cases}-\dfrac{RT}{c_{\text{T}}^{0}\mathscr{D}_{ij}}\>\>\>\text{if}\>\>i\neq
j\vspace{6pt}\\\ \dfrac{RT}{c_{\text{T}}^{0}}\displaystyle{\sum_{k\neq
i}^{n}}\dfrac{c_{k}^{0}}{\mathscr{D}_{ik}c_{j}^{0}}\>\>\>\text{if}\>\>i=j,\end{cases}$
(8.29)

where the neutral composition descriptors $\textbf{c}^{0}$ and
$c_{\textrm{T}}^{0}$ are given in terms of the component concentrations
$\textbf{c}_{\nu}$ by equations (8.1) and (8.2), respectively. Although the
$n$ species concentrations here depend only on $n-1$ component concentrations
from a salt–charge basis, the effect is relatively minor because the
constraint does not affect the rank of $\textbf{M}^{0}$.

When the whole set of Onsager–Stefan–Maxwell equations is put over a
salt–charge basis, the transport coefficients that appear after applying
electroneutrality follow from congruence relation (5.7):

$\textbf{M}_{Z}^{0}=\textbf{Z}\textbf{M}^{0}\textbf{Z}^{\top}=\left[\begin{array}[]{cc}\textbf{M}_{\nu}&\textbf{m}_{z}\\\
\textbf{m}_{z}^{\top}&M_{zz}\end{array}\right].$ (8.30)

The second equality here defines sub-blocks that will be useful for expressing
other macroscopic transport properties in terms of the electroneutral Onsager
drag coefficients over a salt–charge basis, $\textbf{M}_{Z}^{0}$. Sub-block
$\textbf{M}_{\nu}$ is an $\left(n-1\right)\times\left(n-1\right)$ square
matrix, $\textbf{m}_{z}$, an $(n-1)$-dimensional column matrix, and $M_{zz}$,
a scalar.

Recall from section 5 that the change to a salt–charge basis implies
$\textbf{c}_{Z}$ is the sole null eigenvector of $\textbf{M}_{Z}$. After
electroneutrality is adopted, the column $\textbf{c}_{Z}^{0}$, defined as

$\textbf{c}_{Z}^{0}=\left[\begin{array}[]{c}\textbf{c}_{\nu}\\\
0\end{array}\right],$ (8.31)

is a null eigenvector of $\textbf{M}_{Z}^{0}$. This in turn implies that

$\textbf{M}_{\nu}\textbf{c}_{\nu}=\textbf{o},$ (8.32)

which is to say, the column of component concentrations is a null eigenvector
of the symmetric sub-block $\textbf{M}_{\nu}$, and also that

$\textbf{m}_{z}^{\top}\textbf{c}_{\nu}=0,$ (8.33)

so that the column $\textbf{m}_{z}$ is orthogonal to the component
concentrations.

Once the transport matrix from equation (8.30) has been incorporated, the
Onsager–Stefan–Maxwell equations over a salt–charge basis from equation (5.6)
become

$-\left[\begin{array}[]{c}\nabla\bm{\mu}_{\nu}^{0}\\\
F\left\|\textbf{z}\right\|\nabla\Phi_{z}^{0}\end{array}\right]=\left[\begin{array}[]{cc}\textbf{M}_{\nu}&\textbf{m}_{z}\\\
\textbf{m}_{z}^{\top}&M_{zz}\end{array}\right]\left[\begin{array}[]{c}\vec{\textbf{n}}_{\nu}\\\
\dfrac{\vec{i}}{F\left\|\textbf{z}\right\|}\end{array}\right]$ (8.34)

under local electroneutrality. The transport matrix here is symmetric positive
semidefinite and affords the single null eigenvector $\textbf{c}_{Z}^{0}$,
preserving the spectral structure of the species flux laws. Bear in mind that
the component chemical potential gradients $\nabla\bm{\mu}_{\nu}^{0}$ are
parametrized by a set of thermodynamic factors $\textbf{X}_{\nu}^{0}$
dependent on $\textbf{y}_{\nu}$ and a matrix of material parameters through
equations (8.3), in turn constrained by the symmetry of the excess property
matrix $\bm{\Lambda}_{\nu}^{0}$. Henceforth we call equation (8.34), with
$\nabla\bm{\mu}_{\nu}^{0}$ parametrized by equations (8.24) and (8.3), the
_electroneutral Onsager–Stefan–Maxwell equations_.

##  9 Mass continuity and excess fluxes

Compatibility of the material balances with the familiar mass continuity
equation is ensured by identifying the mass-average (barycentric) velocity as
a combination of the species fluxes. Developing excess fluxes relative to the
mass-average velocity allows an understanding of the non-dissipative
phenomenon of convection, which is distinct from the dissipative diffusion
process.

To incorporate mass into the governing framework, let $\overline{m}_{i}$
represent the molar mass of species $i$, and assemble the set of these masses
into an $n$-dimensional column $\overline{\textbf{m}}$. The mass density
$\rho$ of a multi-species solution is generally given by

$\rho=\sum_{i=1}^{n}\overline{m}_{i}c_{i}=\overline{\textbf{m}}^{\top}\textbf{c}.$
(9.1)

Further let $\psi_{i}$ represent a weighting factor with units of inverse
molarity (or, equivalently, molar volume),

$\psi_{i}=\frac{\overline{m}_{i}}{\rho},~{}~{}~{}~{}~{}\textrm{or}~{}~{}~{}~{}~{}\bm{\psi}=\frac{\overline{\textbf{m}}}{\rho},$
(9.2)

where the second relationship summarizes the definition for all $i$ by
introducing an $n$-dimensional column matrix $\bm{\psi}$ with entries
$\psi_{i}$. Note that equation (9.1) implies a general constraint among the
weighting factors,

$\bm{\psi}^{\top}\textbf{c}=1,$ (9.3)

that is, the species mass fractions $\psi_{i}c_{i}$ sum to unity. The mass-
average velocity $\vec{v}$ is generally defined as

$\vec{v}=\sum_{i=1}^{n}\psi_{i}\vec{N}_{i}=\bm{\psi}^{\top}\vec{\textbf{n}}.$
(9.4)

With equations (9.1), (9.3), and (9.4) in hand, one can multiply equation
(5.1) through from the left by the constant row matrix $\rho\bm{\psi}^{\top}$
($=\overline{\textbf{m}}^{\top}$) to get

$\frac{\partial\rho}{\partial t}=-\nabla\cdot\left(\rho\vec{v}\right),$ (9.5)

showing consistency of the material balances with the common understanding of
mass continuity.

Defining a bulk velocity such as $\vec{v}$ distinguishes convective mass
transport from diffusion as follows. Let the excess molar flux of species $i$
relative to the mass-average velocity, $\vec{J}_{i}$, be defined as

$\vec{J}_{i}=\vec{N}_{i}-c_{i}\vec{v},~{}~{}~{}~{}~{}\textrm{or}~{}~{}~{}~{}~{}\vec{\textbf{j}}=\vec{\textbf{n}}-\vec{v}\,\textbf{c},$
(9.6)

where the latter equation introduces the $n$-dimensional column
$\vec{\textbf{j}}$ whose entries are each of the $\vec{J}_{i}$ in sequence.
Insertion of definition (9.6) into equation (5.5) gives

$-\nabla\bm{\mu}=\textbf{M}\left(\vec{\textbf{j}}+\vec{v}\,\textbf{c}\right)=\textbf{M}\vec{\textbf{j}},$
(9.7)

which is to say, because c is a null eigenvector of M, the physical content of
the Onsager–Stefan–Maxwell equations is independent of convection.

Importantly, the very presence of a bulk velocity constrains excess fluxes. In
light of equations (9.3) and (9.4), multiplying equation (9.6) through by
$\bm{\psi}^{\top}$ reveals that

$\bm{\psi}^{\top}\vec{\textbf{j}}=\vec{0},$ (9.8)

establishing a kinematic relation among the set of excess molar species fluxes
relative to the barycentric velocity.

Electroneutrality affects density and the kinematic relation in different
ways. Application of electroneutrality to equation (9.1) and insertion of
equation (8.1) show that density is a function of $n-1$ component
concentrations,

$\rho^{0}=\overline{\textbf{m}}^{\top}\textbf{c}^{0}=\overline{\textbf{m}}_{\nu}^{\top}\textbf{c}_{\nu},$
(9.9)

in which $\overline{\textbf{m}}_{\nu}=\textbf{N}^{\top}\overline{\textbf{m}}$
is an $(n-1)$-dimensional column comprising the component molar masses. Excess
species fluxes then transform into the salt–charge basis as

$\vec{\textbf{j}}_{Z}=\textbf{Z}^{-\top}\vec{\textbf{j}}=\left[\begin{array}[]{c}\vec{\textbf{j}}_{\nu}\\\
\dfrac{\vec{i}-\rho_{\textrm{e}}\vec{v}}{F\left\|\textbf{z}\right\|}\end{array}\right],$
(9.10)

because $\rho_{\textrm{e}}=F\textbf{z}^{\top}\textbf{c}$. Imposing local
electroneutrality on $\vec{\textbf{j}}_{Z}$ here shows that the current
density, embodied by the last entry of $\vec{\textbf{j}}_{Z}^{\,0}$, becomes
independent of convection. This entry does not vanish entirely, however, and
consequently electroneutrality does not reduce the dimensionality of the space
of excess fluxes in the same way that it reduces the domain of $\rho^{0}$.

Defining $\bm{\psi}_{Z}=\textbf{Z}\bm{\psi}$, one can change basis to write
kinematic relation (9.8) as

$\bm{\psi}_{Z}^{\top}\vec{\textbf{j}}_{Z}=\vec{0},$ (9.11)

showing that the column of excess component fluxes is orthogonal to
$\bm{\psi}_{Z}^{\top}$. Electroneutral species fluxes $\vec{\textbf{j}}^{\,0}$
depend on the component excess fluxes over a salt–charge basis,
$\vec{\textbf{j}}_{\nu}$, as

$\vec{\textbf{j}}^{\,0}=\textbf{Z}^{\top}\vec{\textbf{j}}_{Z}^{\,0}=\textbf{N}\vec{\textbf{j}}_{\nu}+\dfrac{\vec{i}}{F\left\|\textbf{z}\right\|}\cdot\dfrac{\textbf{z}}{\left\|\textbf{z}\right\|},$
(9.12)

a form that notably includes the current, as well as the excess component
fluxes.

##  10 Flux-explicit transport laws

Typical constitutive laws for electrolytic transport, including the standard
concentrated solution theory and the Nernst–Planck equations, are written in
flux-explicit forms. Laws giving component flux and current density in terms
of the gradients of composition and voltage can be derived from force-explicit
Onsager–Stefan–Maxwell equations by exploiting the inversion formulas reported
previously by the authors [18]. Generally, flux-explicit laws for diffusion
must be expressed using excess species fluxes in order to exclude the effects
of convection; transport-law inversion therefore necessitates the selection of
a representative bulk velocity.

Section 9 set out the preliminary definitions needed to develop flux-explicit
transport laws relative to the mass-average velocity by the inversion method
of Van-Brunt et al. [18]. This method begins with a statement of the
electroneutral Onsager–Stefan–Maxwell equations (8.34) in terms of excess
fluxes. Through equation (9.10), they become

$-\left[\begin{array}[]{c}\nabla\bm{\mu}_{\nu}^{0}\\\
F\left\|\textbf{z}\right\|\nabla\Phi_{z}^{0}\end{array}\right]=\left[\begin{array}[]{cc}\textbf{M}_{\nu}&\textbf{m}_{z}\\\
\textbf{m}_{z}^{\top}&M_{zz}\end{array}\right]\left[\begin{array}[]{c}\vec{\textbf{j}}_{\nu}\\\
\dfrac{\vec{i}}{F\left\|\textbf{z}\right\|}\end{array}\right]$ (10.1)

in terms of $\vec{\textbf{j}}_{\nu}$. To implement the inversion, we first
observe that within the flux-explicit form of these laws, the electroneutral
excess component fluxes $\textbf{j}_{\nu}$ and $\vec{i}$ must be constrained
such that ${\bm{\psi}_{Z}^{0}}^{\top}\textbf{j}_{Z}^{0}=0$, wherein the
kinematic relation over the salt–charge basis, $\bm{\psi}_{Z}^{0}$, is defined
such that $\bm{\psi}_{Z}^{0}=\textbf{Z}\bm{\psi}^{0}$. (Here $\bm{\psi}^{0}$
expresses the basic kinematic relation that constrains excess species fluxes
relative to the barycentric velocity under local electroneutrality, identified
by replacing $\rho$ with $\rho^{0}$ in equation (9.2).) Therefore, the flux-
explicit transport laws must involve a coefficient matrix $\textbf{L}_{Z}^{0}$
with null eigenvector $\bm{\psi}_{Z}^{0}$. This matrix is formed uniquely from
$\textbf{M}_{Z}^{0}$ through the limit process

$\textbf{L}_{Z}^{0}=\lim_{\gamma\to
0}\left(\textbf{M}_{Z}^{0}+\frac{\bm{\psi}_{Z}^{0}{\bm{\psi}_{Z}^{0}}^{\top}}{\gamma}\right)^{-1},$
(10.2)

or, equivalently, from the algebraic equation

$\textbf{L}_{Z}^{0}=\left(\textbf{M}_{Z}^{0}+\gamma\bm{\psi}_{Z}^{0}{\bm{\psi}_{Z}^{0}}^{\top}\right)^{-1}-\frac{\textbf{c}_{Z}^{0}\left(\textbf{c}_{Z}^{0}\right)^{\top}}{\gamma},$
(10.3)

whose result is identical for any nonzero value of the augmentation parameter
$\gamma$ [18]. After partitioning $\textbf{L}_{Z}^{0}$ into the block form

$\textbf{L}_{Z}^{0}=\left[\begin{array}[]{cc}\textbf{L}_{\nu}&\textbf{l}_{z}\\\
\textbf{l}_{z}^{\top}&L_{zz}\end{array}\right],$ (10.4)

one finds

$\left[\begin{array}[]{c}\vec{\textbf{j}}_{\nu}\\\
\dfrac{\vec{i}}{F\left\|\textbf{z}\right\|}\end{array}\right]=-\left[\begin{array}[]{cc}\textbf{L}_{\nu}&\textbf{l}_{z}\\\
\textbf{l}_{z}^{\top}&L_{zz}\end{array}\right]\left[\begin{array}[]{c}\nabla\bm{\mu}_{\nu}^{0}\\\
F\left\|\textbf{z}\right\|\nabla\Phi_{z}^{0}\end{array}\right],$ (10.5)

the desired flux-explicit form of the electroneutral Onsager transport laws
with respect to a salt–charge basis.

##  11 Distinguishing ionic conductivity, migration, and molecular diffusion

One can leverage the flux-explicit electroneutral Onsager transport laws to
generalize the process by which concentrated solution theory distinguishes the
various dissipative phenomena that drive electrochemical mass transport [1].
First and foremost, the last row of equation (10.5) establishes the MacInnes
equation—a current–voltage relation that can be viewed as a modified form of
Ohm’s law accounting for concentration overpotential [21]. Identifying the
_ionic conductivity_ $\kappa^{0}$ as

$\kappa^{0}=F^{2}\left\|\textbf{z}\right\|^{2}L_{zz},$ (11.1)

and further letting an $(n-1)$-entry column $\bm{\xi}$ represent a set of what
we will call _component migration coefficients_ , defined by

$\bm{\xi}=\frac{1}{L_{zz}}\textbf{l}_{z},$ (11.2)

the MacInnes equation can be written as

$\vec{i}=-\kappa^{0}\nabla\Phi_{z}^{0}-\frac{\kappa^{0}\bm{\xi}^{\top}\nabla\bm{\mu}_{\nu}^{0}}{F\left\|\textbf{z}\right\|}$
(11.3)

in terms of the salt–charge potential and the component chemical potentials
with respect to the chosen salt–charge basis. The first term on the right
describes charge flow by ionic conduction, i.e., current driven by electric-
potential gradients, in the absence of composition gradients.

Migration is defined as the excess mass flux driven by electric current flow
in the absence of composition gradients. To isolate this effect, one writes
the species flux laws in a partially inverted form, with current density
replacing the electric-potential gradient. Newman does this by eliminating one
of the species fluxes with Faraday’s law (4.10) before performing a
pseudoinversion of the electroneutral Onsager–Stefan–Maxwell equations [1].
Equivalently, one can rearrange the Onsager flux laws from equation (10.5), by
using the MacInnes equation—the last row of equation (10.5)—to replace
$\nabla\Phi_{z}^{0}$ with $\vec{i}/F\left\|\textbf{z}\right\|$ in the laws for
$\vec{\textbf{j}}_{\nu}$ expressed by the remaining rows. Let the symmetric
matrix

$\overline{\textbf{L}}_{\nu}=\frac{RT}{c_{\textrm{T}}^{0}}\left(\textbf{L}_{\nu}-\frac{\kappa^{0}\bm{\xi}\bm{\xi}^{\top}}{F^{2}\left\|\textbf{z}\right\|^{2}}\right)$
(11.4)

define the set of _Onsager component diffusivities_ relative to the mass-
average velocity, $\overline{\textbf{L}}_{\nu}$. The excess component fluxes
become

$\vec{\textbf{j}}_{\nu}=-\frac{c_{\textrm{T}}^{0}\overline{\textbf{L}}_{\nu}}{RT}\nabla\bm{\mu}_{\nu}^{0}+\bm{\xi}\frac{\vec{i}}{F\left\|\textbf{z}\right\|},$
(11.5)

in which $\overline{\textbf{L}}_{\nu}$ has the standard units of area per
time. Here, the first term on the right describes the component flux
associated with diffusion, and the second, migration.

Note that in the absence of current density, equation (11.5) reduces to a
generalized form of Fick’s law that expresses component excess fluxes solely
in terms of the symmetric component-diffusivity matrix
$\overline{\textbf{L}}_{\nu}$. Miller deployed transport laws with this
structure to describe solutions containing multiple dissolved salts in some of
his experimental tests of the Onsager reciprocal relations [22], which
confirmed the symmetry of $\overline{\textbf{L}}_{\nu}$.

When quantifying migration coefficients, one should bear in mind that their
values must be constrained to ensure that $\bm{\psi}_{Z}^{0}$ is a null
eigenvector of $\textbf{L}_{Z}^{0}$. By partitioning $\bm{\psi}_{Z}^{0}$ into
two parts $\bm{\psi}_{\nu}$ and $\psi_{z}$ such that

$\bm{\psi}_{Z}^{0}=\left[\begin{array}[]{c}\textbf{N}^{\top}\bm{\psi}^{0}\\\
\dfrac{\textbf{z}^{\top}\bm{\psi}^{0}}{\left\|\textbf{z}\right\|}\end{array}\right]=\left[\begin{array}[]{c}\bm{\psi}_{\nu}\\\
\psi_{z}\end{array}\right],$ (11.6)

one can write

$\bm{\psi}_{\nu}^{\top}\bm{\xi}=-\psi_{z}$ (11.7)

to phrase this constraint entirely in terms of species charges and molar
masses.

Within concentrated solution theory, transport laws are generally written for
species fluxes, rather than component fluxes [1]. By leveraging equation
(9.12), the species fluxes can be separated into terms associated with
diffusion and migration using the properties defined above. Inserting
equations (11.5) yields

$\vec{\textbf{j}}^{\,0}=-\frac{\textbf{N}\overline{\textbf{L}}_{\nu}c_{\textrm{T}}^{0}}{RT}\nabla\bm{\mu}_{\nu}^{0}+\left(\textbf{N}\bm{\xi}+\frac{\textbf{z}}{\left\|\textbf{z}\right\|}\right)\frac{\vec{i}}{F\left\|\textbf{z}\right\|}.$
(11.8)

As well as distinguishing diffusion from migration at the species level, this
expression helps to elucidate the spectral structure of the diffusivity matrix
$\overline{\textbf{L}}_{\nu}$. Multiplying through from the left by
$(\bm{\psi}^{0})^{\top}$, applying the electroneutral form of equation (9.8),
introducing $\bm{\psi}_{\nu}$ and $\psi_{z}$ with equation (11.6), and
applying equation (11.7) show that

$-\frac{c_{\textrm{T}}^{0}\bm{\psi}_{\nu}^{\top}\overline{\textbf{L}}_{\nu}}{RT}\nabla\bm{\mu}_{\nu}^{0}=\vec{0}$
(11.9)

for any $\nabla\bm{\mu}_{\nu}^{0}$. Thus, because
$\overline{\textbf{L}}_{\nu}$ is symmetric, it must afford $\bm{\psi}_{\nu}$
as a null eigenvector:
$\overline{\textbf{L}}_{\nu}\bm{\psi}_{\nu}=\textbf{o}$.444If
$\nabla\bm{\mu}_{\nu}^{0}$ is restricted to a space that satisfies the
Gibbs–Duhem relation from equation (8.14), equation (11.9) permits that
$\overline{\textbf{L}}_{\nu}\bm{\psi}_{\nu}=k\bm{c}_{\nu}$ for some $k$. We
choose $k=0$ for consistency with the model formed by treating component
transport in the absence of current as transport within a non-electrolytic
system containing $n-1$ uncharged species, whose flux-explicit form matches
equation (11.5) with $\bm{\xi}=\textbf{o}$.

A property count is helpful here. Symmetry of the $(n-1)\times(n-1)$ matrix
$\overline{\textbf{L}}_{\nu}$ reduces its number of independent entries by
$\tfrac{1}{2}\left(n-1\right)\left(n-2\right)$; the fact that it affords
$\bm{\psi}_{\nu}$ as a null eigenvector adds $(n-1)$ additional constraints.
Thus $\tfrac{1}{2}\left(n-1\right)\left(n-2\right)$ of the entries in
$\overline{\textbf{L}}_{\nu}$ are independently specifiable. The
$(n-1)$-dimensional column $\bm{\xi}$ is constrained through species molar
masses and equivalent charges by equation (11.7), leaving $(n-2)$ independent
migration coefficients. Finally, equation (11.1) defines a single ionic
conductivity. Summing up, the symmetric, positive-semidefinite flux-explicit
Onsager matrix with respect to the salt–charge basis breaks down as

$\textbf{L}_{Z}^{0}=\frac{c_{\textrm{T}}^{0}}{RT}\left[\begin{array}[]{cc}\overline{\textbf{L}}_{\nu}&\textbf{o}\\\
\textbf{o}^{\top}&0\end{array}\right]+\frac{\kappa^{0}}{F^{2}\left\|\textbf{z}\right\|^{2}}\left[\begin{array}[]{cc}\bm{\xi}\bm{\xi}^{\top}&\bm{\xi}\\\
\bm{\xi}^{\top}&1\end{array}\right],$ (11.10)

in which $\overline{\textbf{L}}_{\nu}=\overline{\textbf{L}}_{\nu}^{\top}$,
$\overline{\textbf{L}}_{\nu}\bm{\psi}_{\nu}=\textbf{o}$, and
$\bm{\psi}_{\nu}^{\top}\bm{\xi}+\psi_{z}=0$, with $\bm{\psi}_{\nu}$ and
$\psi_{z}$ defined in terms of stoichiometry in the simple association
equilibria, species molar masses, species charges, and the electroneutral
density $\rho^{0}$ through equations (9.2) and (11.6). Consequently,
$\textbf{L}_{Z}^{0}$ depends on $\tfrac{1}{2}n(n-1)$ independent transport
coefficients. This equals the number of Stefan–Maxwell diffusivities that
underpin the original transport matrix $\textbf{M}^{0}$ in equation (8.29).

Recently molecular-dynamics techniques have been deployed to measure the naive
Onsager diffusion matrix $\textbf{L}^{0}$, which sits relative to the mass-
average velocity and appears in the inverted form of equation (9.7) [23,
24],555A molecular dynamics simulation can be described with equation (11.11)
if the simulated control volume conserves charge and has characteristic
dimensions much larger than the Debye length.

$\vec{\textbf{j}}^{\,0}=-\textbf{L}^{0}\nabla\bm{\mu}^{0}.$ (11.11)

The transport properties that make up $\textbf{L}^{0}$ relate to
$\textbf{M}^{0}$ directly through

$\textbf{M}^{0}=\lim_{\gamma\to
0}\left[\textbf{L}^{0}+\frac{\textbf{c}^{0}\left(\textbf{c}^{0}\right)^{\top}}{\gamma}\right]^{-1},$
(11.12)

where we have used the limit form of the inversion described by Van-Brunt et
al. [18]. Through equation (11.10), the congruence transformation

$\textbf{L}_{Z}^{0}=\textbf{Z}^{-\top}\textbf{L}^{0}\textbf{Z}^{-1}$ (11.13)

shows how the component diffusivities, migration coefficients, and ionic
conductivity over a salt–charge basis contribute to the naive Onsager matrix
$\textbf{L}^{0}$.

##  12 Alternative convective velocities

Up to now we have adopted a convention of defining $\textbf{L}_{Z}^{0}$ and
$\vec{\textbf{j}}^{\,0}$ relative to the mass-average velocity, but it is
straightforward to change the kinematic constraint that defines the convective
reference velocity in equation (9.4) to any other $\bm{\psi}^{\prime}$
satisfying equation (9.3). (For example, letting
$\bm{\psi}^{\prime}=\textbf{1}/c_{\textrm{T}}$ defines the mole-average
velocity.) The conversion of the set of excess molar fluxes relative to the
mass-average velocity, composed of $\vec{J}_{i}$, to excess fluxes relative to
another velocity, composed of $\vec{J}_{i}^{\,\psi^{\prime}}$, is accomplished
by the projection operation

$\vec{\textbf{j}}^{\psi^{\prime}}=\left(\textbf{I}-\textbf{c}\bm{\psi}^{\prime\top}\right)\vec{\textbf{j}}.$
(12.1)

Observe that a similar projection can be used to move from excess fluxes that
satisfy kinematic relation $\bm{\psi}^{\prime}$ back to fluxes in excess of
convection at the mass-average reference velocity:

$\vec{\textbf{j}}=\left(\textbf{I}-\textbf{c}\bm{\psi}^{\top}\right)\vec{\textbf{j}}^{\psi^{\prime}}.$
(12.2)

Transformation (12.1) can be expressed in terms of excess fluxes over a
salt–charge basis as

$\vec{\textbf{j}}^{\,\psi^{\prime}}_{Z}=\left[\textbf{I}-\textbf{c}_{Z}(\bm{\psi}^{\prime}_{Z})^{\top}\right]\vec{\textbf{j}}_{Z},$
(12.3)

where $\bm{\psi}^{\prime}_{Z}=\textbf{Z}\bm{\psi}^{\prime}$.

Under electroneutrality, the kinematic relation $\bm{\psi}_{Z}$ is replaced by
$\bm{\psi}_{Z}^{0}$. Given any other kinematic relation $\bm{\psi}^{\prime}$
defined under the electroneutrality approximation, one can exploit equation
(10.2) or (10.3), with $\bm{\psi}^{\prime}_{Z}$ in place of
$\bm{\psi}_{Z}^{0}$, to identify the flux-explicit electroneutral transport
relations relative to $\bm{\psi}^{\prime}$, as

$\vec{\textbf{j}}^{\,0,\psi^{\prime}}_{Z}=-\textbf{L}^{0,\psi^{\prime}}_{Z}\nabla\bm{\mu}^{0}_{Z},$
(12.4)

in which

$\textbf{L}^{0,\psi^{\prime}}_{Z}=\left[\textbf{I}-\textbf{c}^{0}_{Z}(\bm{\psi}^{\prime}_{Z})^{\top}\right]\textbf{L}^{0}_{Z}\left[\textbf{I}-\bm{\psi}^{\prime}_{Z}(\textbf{c}^{0}_{Z})^{\top}\right]$
(12.5)

establishes a congruence between $\textbf{L}^{0,\psi^{\prime}}_{Z}$ and
$\textbf{L}^{0}_{Z}$. One can also compute $\textbf{L}^{0,\psi^{\prime}}_{Z}$
directly from $\textbf{M}_{Z}^{0}$, using

$\textbf{L}^{0,\psi^{\prime}}_{Z}=\lim_{\gamma\to
0}\left(\textbf{M}_{Z}^{0}+\frac{\bm{\psi}^{\prime}_{Z}{\bm{\psi}^{\prime}_{Z}}^{\top}}{\gamma}\right)^{-1},$
(12.6)

which is analogous to equation (10.2). In the new reference frame for the
excess fluxes established by $\bm{\psi}^{\prime}$, diffusivities and migration
coefficients can be identified by partitioning
$\textbf{L}^{0,\psi^{\prime}}_{Z}$ as

$\displaystyle\textbf{L}_{Z}^{0,\psi^{\prime}}=$
$\displaystyle\frac{c_{\textrm{T}}^{0}}{RT}\left[\begin{array}[]{cc}\overline{\textbf{L}}_{\nu}^{\psi^{\prime}}&\textbf{o}\\\
\textbf{o}^{\top}&0\end{array}\right]$ (12.9)
$\displaystyle+\frac{\kappa^{0}}{F^{2}\left\|\textbf{z}\right\|^{2}}\left[\begin{array}[]{cc}\bm{\xi}^{\psi^{\prime}}\big{(}\bm{\xi}^{\psi^{\prime}}\big{)}^{\top}&\bm{\xi}^{\psi^{\prime}}\\\
\big{(}\bm{\xi}^{\psi^{\prime}}\big{)}^{\top}&1\end{array}\right],$ (12.12)

similar to how $\textbf{L}_{Z}^{0}$ is partitioned in equation (11.10). Note
that the second matrix on the right of equation (12.9) is proportional to the
outer product of $[\bm{\xi}^{\psi^{\prime}},1]^{\top}$ with itself, and is
consequently positive semidefinite in general.

One can write excess component fluxes relative to the alternative convective
velocity defined by $\bm{\psi}^{\prime}$, and with respect to the salt–charge
basis, as

$\vec{\textbf{j}}_{\nu}^{\,\psi^{\prime}}=-\frac{c_{\textrm{T}}^{0}\overline{\textbf{L}}_{\nu}^{\psi^{\prime}}}{RT}\nabla\bm{\mu}_{\nu}^{0}+\bm{\xi}^{\,\psi^{\prime}}\frac{\vec{i}}{F\left\|\textbf{z}\right\|},$
(12.13)

which derives from equation (11.5). Here, the component diffusivities relative
to the average velocity determined by $\bm{\psi}^{\prime}$ are defined as

$\overline{\textbf{L}}_{\nu}^{\psi^{\prime}}=\left[\textbf{I}_{n-1}-\textbf{c}_{\nu}(\bm{\psi}^{\prime}_{\nu})^{\top}\right]\overline{\textbf{L}}_{\nu}\left[\textbf{I}_{n-1}-\bm{\psi}^{\prime}_{\nu}(\textbf{c}_{\nu})^{\top}\right],$
(12.14)

and

$\bm{\xi}^{\psi^{\prime}}=\left[\textbf{I}_{n-1}-\textbf{c}_{\nu}(\bm{\psi}^{\prime}_{\nu})^{\top}\right]\bm{\xi}-\psi_{z}^{\prime}\textbf{c}_{\nu}$
(12.15)

defines migration coefficients in the new frame. Similar to the projection in
equation (12.2), one can create transformations that change the reference
frame from velocity ${\bm{\psi}^{\prime}}^{\top}\vec{\textbf{n}}$ back to the
mass-average velocity by replacing $\bm{\psi}^{\prime}_{\nu}$ with
$\bm{\psi}_{\nu}$ and $\psi^{\prime}_{z}$ with $\psi_{z}$ in equations (12.14)
and (12.15), as well as swapping $\overline{\textbf{L}}_{\nu}^{\psi^{\prime}}$
with $\overline{\textbf{L}}_{\nu}$ and $\bm{\xi}^{\psi^{\prime}}$ with
$\bm{\xi}$.

In Newman’s implementations of concentrated solution theory [1], one of the
species velocities—say, that of species $m$—is used as the convective
reference. Thus the coefficients that Newman identifies in the flux-explicit
equations follow from taking
$\bm{\psi}^{\prime}_{Z}=\textbf{Z}\textbf{i}_{m}/c^{0}_{m}$ in equation
(12.5). This produces species flux laws in the form

$\vec{\textbf{j}}^{\,0,\psi^{\prime}}=-\frac{c_{\textrm{T}}^{0}\textbf{N}\overline{\textbf{L}}_{\nu}^{\psi^{\prime}}}{RT}\nabla\bm{\mu}_{\nu}^{0}+\left(\textbf{N}\bm{\xi}^{\psi^{\prime}}+\frac{\textbf{z}}{\left\|\textbf{z}\right\|}\right)\frac{\vec{i}}{F\left\|\textbf{z}\right\|}.$
(12.16)

Newman’s development also produces a slightly different MacInnes equation,

$\vec{i}=-\kappa^{0}\nabla\Phi_{z}^{0}-\frac{\kappa^{0}}{F\left\|\textbf{z}\right\|}\big{(}\bm{\xi}^{\psi^{\prime}}\big{)}^{\top}\nabla\bm{\mu}_{\nu}^{0},$
(12.17)

wherein the migration coefficients are constrained by kinematic relation
$\bm{\psi}^{\prime}$, but $\kappa^{0}$ remains unaffected because its value is
independent of the choice of convective velocity.

Note that the embedding of an arbitrary reference velocity in flux-explicit
transport matrices (e.g., $\textbf{L}^{0}$) makes properties derived from them
somewhat ambiguous. When tabulating transport parameters, it may be preferable
to work with the electroneutral transport matrix $\textbf{M}^{0}_{Z}$, which
conveys frame-invariant information. The computation of $\textbf{M}_{\nu}$ can
be accomplished readily given a flux-explicit form; in fact, one can use the
procedure of Van-Brunt et al. [18] and a block inversion identity to show for
any $\bm{\psi}^{\prime}$ that

$\textbf{M}_{\nu}=\frac{RT}{c^{0}_{\textrm{T}}}\lim_{\gamma\rightarrow
0}\left(\overline{\textbf{L}}^{\,\psi^{\prime}}_{\nu}+\frac{\textbf{c}_{\nu}\textbf{c}_{\nu}^{\top}}{\gamma}\right)^{-1},$
(12.18)

and also that

$\textbf{m}_{z}=\textbf{M}_{\nu}\bm{\xi}^{\psi^{\prime}}~{}~{}\textrm{and}~{}~{}M_{zz}=\frac{F^{2}\|\textbf{z}\|^{2}}{\kappa^{0}}+(\bm{\xi}^{\psi^{\prime}})^{\top}\textbf{M}_{\nu}\bm{\xi}^{\psi^{\prime}},$
(12.19)

elucidating the structure of $\textbf{M}_{Z}^{0}$ through equation (8.30).
Despite its computation from a transport matrix
$\textbf{L}_{Z}^{0,\psi^{\prime}}$, it should be emphasized that
$\textbf{M}_{Z}^{0}$ enjoys uniqueness and is independent of the convective
velocity. Notably, the submatrix $\textbf{M}_{\nu}$ depends only on the
Onsager component diffusivities, and the only matrix inversion that needs to
be performed is the one in equation (12.18). The entries of $\textbf{m}_{z}$
are zero in the absence of migration, and $M_{zz}$ matches the ionic
resistivity plus a quadratic correction that also vanishes when migration
does.

##  13 Transference numbers

Equations (11.3) and (11.8) provide a pair of constitutive equations similar,
but not identical, to the familiar equations from concentrated solution theory
and Nernst–Planck dilute-solution theory. To achieve complete agreement, one
must consider how migration coefficients relate to the transference numbers
more familiar from electrochemistry [21, 1, 25]. Let the transference number
of species $i$ relative to the mass-average velocity be $t_{i}$. These can be
formed into an $n$-dimensional column matrix t, related to the migration
coefficients as

$\textbf{t}=\frac{\operatorname{diag}\left(\textbf{z}\right)}{\left\|\textbf{z}\right\|}\left(\textbf{N}\bm{\xi}+\frac{\textbf{z}}{\left\|\textbf{z}\right\|}\right).$
(13.1)

Multiplication of this result through by $\textbf{1}^{\top}$ shows that
$\textbf{1}^{\top}\textbf{t}=1$, that is, transference numbers sum to unity.
Note that if migration is neglected entirely ($\bm{\xi}=\textbf{o}$), the
right side of equation (13.1) remains non-zero; each transference number
becomes a partition of unity weighted by the corresponding squared species
charge. The migration coefficients we introduced here deviate from this
standard.

Sometimes it is desirable to change the reference velocity relative to which
transference is expressed. Incorporating the alternative migration coefficient
from equation (12.15) into equation (13.1) shows that

$\textbf{t}^{\psi^{\prime}}=\frac{\operatorname{diag}\left(\textbf{z}\right)}{\left\|\textbf{z}\right\|}\left(\textbf{N}\bm{\xi}^{\psi^{\prime}}+\frac{\textbf{z}}{\left\|\textbf{z}\right\|}\right)$
(13.2)

yields the set of species transference numbers relative to the reference
velocity associated with the kinematic relation represented by
$\bm{\psi}^{\prime}$. These also afford the property that
$\textbf{1}^{\top}\textbf{t}^{\psi^{\prime}}=1$.

Because $\operatorname{diag}\left(\textbf{z}\right)$ appears on the right of
equation (13.1), $t_{i}$ vanishes for all $i$ such that $z_{i}=0$. Thus,
transference numbers as traditionally defined only have meaning for charged
species. Importantly, specifying a set of independent transference numbers
does not generally suffice to specify the Onsager transport matrix. Newman
addresses this point by using transference numbers for uncharged species
expressed as ratios $t_{i}/z_{i}$, noting that the ratio is finite for any
species $i$ with $z_{i}=0$ [1]. Indeed, one can see with equation (13.1) that
$t_{i}/z_{i}$ can be written under the assumption that $t_{i}\propto z_{i}$ as

$\frac{t_{i}}{z_{i}}=\frac{1}{\left\|\textbf{z}\right\|}\left(\textbf{N}\bm{\xi}+\frac{\textbf{z}}{\left\|\textbf{z}\right\|}\right)_{i},$
(13.3)

in which the right side is generally nonzero for any $i$.

Rather than taxing ourselves by ignoring apparent divisions by zero when
transference numbers are involved, we propose that it is simpler to think of
transport fundamentally in terms of migration coefficients, because they are
defined unambiguously in terms of Onsager coefficients for species with and
without charge. The nontrivial transference numbers that make up t are
subsidiary quantities that afford a useful physical interpretation, but the
set as a whole does not have as much essential physical content. The
distinction between t and $\bm{\xi}$ becomes particularly important for
electrolytes containing multiple uncharged species, such as the cosolvent
blends commonly used in lithium-ion-battery electrolytes [26, 27], fuel-cell
membranes susceptible to gas crossover [28], or electrolytes for lithium–air
batteries [29]. Within these systems, applied currents can drive migration of
one neutral species relative to another—a dissipative phenomenon distinct from
diffusion that cannot be accounted for by traditional transference numbers.

##  14 Component diffusion coefficients

The entries of the Onsager component-diffusivity matrices
$\overline{\textbf{L}}_{\nu}$ and
$\overline{\textbf{L}}_{\nu}^{\,\psi^{\prime}}$ have the expected units, and
each entry is paired with a particular chemical-potential gradient.
Nevertheless, the dependences of these transport matrices on the convective
velocity obfuscates their connection to the energy dissipation. Newman
addresses this issue by identifying combinations of the Stefan–Maxwell
coefficients between species that quantify the frame-invariant part of the
dissipation arising from component interdiffusion [1, 30]. The process Newman
uses to identify thermodynamic component diffusivities within binary solutions
can be generalized to solutions comprising more than two components as
follows.

Instead of using entries of Onsager’s diffusivity matrix
$\overline{\textbf{L}}_{\nu}$, a frame-invariant set of thermodynamic
component diffusivities is defined by revisiting the Stefan–Maxwell form.
Physically, one expects that in situations where electroneutrality is
maintained ($\rho_{\textrm{e}}=0$) and the current density vanishes uniformly
($\vec{i}=\vec{0}$), the mass-transport model describing an $n$-species
electrolytic phase should be no different than the model describing a
corresponding $(n-1)$-species nonelectrolytic phase in which every uncharged
species is taken to be a component of the original $n$-ary electrolyte. This
equivalence suggests that $\textbf{M}_{\nu}$ possesses the structure

$\frac{c_{\text{T}}^{0}}{RT}\left(\textbf{M}_{\nu}\right)_{ij}=\begin{cases}-\dfrac{\nu_{i}\nu_{j}}{\mathfrak{D}_{ij}}\>\>\>\text{if}\>\>i\neq
j\vspace{6pt}\\\ \displaystyle{\sum_{k\neq
i}^{n-1}}\dfrac{\nu_{i}\nu_{k}c_{\nu,k}}{\mathfrak{D}_{ik}c_{\nu,j}}\>\>\>\text{if}\>\>i=j\end{cases},$
(14.1)

which defines an alternative set of binary diffusion coefficients
$\mathfrak{D}_{ij}$ that describe drag interactions between components $i$ and
$j$ as they move relative to each other. It should be borne in mind that the
indices in $\mathfrak{D}_{ij}$ refer to components in a given salt–charge
basis.

In equation (14.1), the diffusivities that parametrize current-free,
electroneutral component–component interactions are defined analogously to how
species–species interactions are defined in equation (1.2), with a minor
stoichiometric correction to account for the fact that component $i$ within a
salt–charge basis is actually made up of $\nu_{i}$ species through the simple
association equilibria. We shall distinguish component diffusion coefficients
from Stefan–Maxwell coefficients by writing them in Fraktur font, rather than
a script font.

Intrinsic properties of the Onsager matrix $\textbf{M}_{Z}$ demand that the
coefficients defined by equation (14.1) are symmetric in their indices,
$\mathfrak{D}_{ij}=\mathfrak{D}_{ji}$, and that $\mathfrak{D}_{ii}$ is
undefined for all $i$. Hence equation (14.1) defines $\tfrac{1}{2}(n-1)(n-2)$
independent diffusion coefficients—one for each pair of components. The matrix
$\textbf{M}_{\nu}$ must be positive semidefinite, a property guaranteed if
each $\mathfrak{D}_{ij}$ is positive, and by construction it affords
$\textbf{c}_{\nu}$ as a null eigenvector. Furthermore, a block inversion
identity shows that equations (12.6) and (12.9) imply

$\overline{\textbf{L}}_{\nu}^{\psi^{\prime}}=\lim_{\gamma\to
0}\left(\frac{c_{\text{T}}^{0}}{RT}\textbf{M}_{\nu}+\frac{\bm{\psi}^{\prime}_{\nu}{\bm{\psi}^{\prime}_{\nu}}^{\top}}{\gamma}\right)^{-1},$
(14.2)

an inverted form of equation (12.18). Thus
$\overline{\textbf{L}}_{\nu}^{\psi^{\prime}}$, the portion of the
electroneutral Onsager component-diffusivity matrix that describes situations
in which current density vanishes uniformly, relates directly to
$\textbf{M}_{\nu}$.

When casting $\overline{\textbf{L}}_{\nu}$ in terms of component diffusivities
$\mathfrak{D}_{ij}$ with definition (14.1) and inversion formula (14.2), it is
necessary to assume that every entry of $\textbf{c}_{\nu}$ is nonzero almost
everywhere. This must be borne in mind when considering systems where products
of the simple association equilibria are not unique. The situation described
in section 2, where two binary salts without a common ion are dissolved into
an uncharged solvent, provides a simple example. That system can pass locally
through compositions where a particular component concentration vanishes, but
all of the ions retain non-negative concentrations, as described in the
discussion at the end of section 8.1. Such a state nominally induces singular
behavior in one of the diagonal entries defined in equation (14.1). Observe,
however, that the component diffusivities $\mathfrak{D}_{ij}$ are defined
solely in terms of the off-diagonal entries of $\textbf{M}_{\nu}$. Therefore,
although one must take care when applying inversion formula (14.2), the
diffusion coefficients $\mathfrak{D}_{ij}$ remain well posed across the
physically available range of $\textbf{c}_{\nu}$.

##  15 Reference electrodes

Introducing the concept of a reference electrode allows the salt–charge
potential to be related unambiguously to an experimentally measurable voltage,
despite the degree of freedom in $\Phi_{z}^{0}$ demanded by Guggenheim’s
principle. A heterogeneous half-reaction suitable for use within a reference
electrode is generally written as

$\sum_{i=1}^{n}s_{i}\textrm{X}_{i}^{z_{i}}+\sum_{j}s^{\prime}_{j}\textrm{Y}_{j}\leftrightharpoons
n_{\textrm{e}^{-}}\textrm{e}^{-},$ (15.1)

in which $\textrm{X}_{i}^{z_{i}}$ is the chemical symbol of species $i$,
$\textrm{Y}_{j}$ is the symbol of uncharged entity $j$ that forms in an
immediately adjacent phase, $\textrm{e}^{-}$ denotes electrons, the integer
$s_{i}$ is the stoichiometric coefficient of species $i$ in the half-reaction,
integer $s^{\prime}_{j}$ is the stoichiometric coefficient of
$\textrm{Y}_{j}$, and the positive integer $n_{\textrm{e}^{-}}$ is the number
of electrons involved. Since half-reactions balance charge, the set
$\\{s_{i}\\}_{i=1}^{n}$ must generally satisfy

$\sum_{i=1}^{n}s_{i}z_{i}=\textbf{s}^{\top}\textbf{z}=n_{\textrm{e}^{-}}z_{\textrm{e}^{-}},$
(15.2)

where $z_{\textrm{e}^{-}}$ is the equivalent charge of electrons. (If a half-
reaction is written as a reduction, then $s_{i}$ or $s^{\prime}_{i}$ is
negative when species $i$ is a reactant and positive when it is a product.)
Note that the matrix form in the middle of this expression indicates that by
definition, the column s representing half-reaction stoichiometry is not
orthogonal to the column z of species charges.

The electrode potential $\Phi^{\ominus}$ associated with reference-electrode
reaction (15.1) is defined such that [1, 17]

$n_{\textrm{e}^{-}}Fz_{\textrm{e}^{-}}\nabla\Phi^{\ominus}=\sum_{i=1}^{n}s_{i}\nabla\mu_{i}.$
(15.3)

After defining
$\textbf{s}_{Z}=\textbf{Z}^{-\top}\textbf{s}=[\textbf{s}_{\nu},s_{z}]^{\top}$
and assuming an electroneutral state, this is written equivalently as

$n_{\textrm{e}^{-}}Fz_{\textrm{e}^{-}}\nabla\Phi^{\ominus}=\textbf{s}_{Z}^{\top}\nabla\bm{\mu}_{Z}^{0}=\textbf{s}_{\nu}^{\top}\nabla\bm{\mu}_{\nu}^{0}+s_{z}F\left\|\textbf{z}\right\|\nabla\Phi_{z}^{0},$
(15.4)

and since

$n_{\textrm{e}^{-}}z_{\textrm{e}^{-}}=\textbf{s}^{\top}\textbf{z}=\textbf{s}_{Z}^{\top}\textbf{Z}\textbf{z}=s_{z}\left\|\textbf{z}\right\|,$
(15.5)

it follows that

$\nabla\Phi^{\ominus}=\frac{\textbf{s}_{\nu}^{\top}\nabla\bm{\mu}_{\nu}^{0}}{n_{\textrm{e}^{-}}Fz_{\textrm{e}^{-}}}+\nabla\Phi_{z}^{0},$
(15.6)

a relation that allows the voltage of a given reference electrode to be
inferred after one has computed the component and salt-charge potentials. Note
that at uniform temperature, pressure, and composition,

$\left(\nabla\Phi^{\ominus}\right)_{T,p,\textbf{y}_{\nu}}=\left(\nabla\Phi_{z}^{0}\right)_{T,p,\textbf{y}_{\nu}}$
(15.7)

for any reference-electrode half-reaction, showing the intrinsic utility of
the salt–charge potential.

##  16 Examples

###  16.1 The binary electrolyte

The model of a binary electrolytic solution—a liquid comprising a single
simple salt dissolved in a single uncharged solvent—is the simplest and most
common implementation of concentrated solution theory. Newman’s identification
of transport laws for this system begins by adopting Onsager–Stefan–Maxwell
equations (1.1) with the transport coefficients defined in equation (1.2). He
uses these in combination with the electroneutrality approximation (1.3),
Guggenheim condition (2.2) for the salt, and Faraday’s law (4.10) to derive
species flux laws in the form of equation (12.16) and a MacInnes equation in
the form of equation (12.17), casting both in terms of properties relative to
the solvent velocity. The development involves apparently ad hoc manipulation
of the transport equations to obtain compact forms of the flux-explicit laws.

We now demonstrate how Newman’s process can be formalized with a salt–charge
basis under electroneutrality. Let $0$, $+$, and $-$ distinguish the solvent,
cation, and anion species, respectively, and subscripts $0,\text{e}$ denote
the solvent and salt components. Set the stoichiometric coefficients for the
reactants in the simple association equilibrium forming the salt as

$\nu_{+}=-\frac{\nu z_{-}}{z_{+}-z_{-}},\>\>\>\>\>\>\nu_{-}=\frac{\nu
z_{+}}{z_{+}-z_{-}},$ (16.1)

and note that $\nu=\nu_{+}+\nu_{-}$. This framing demonstrates that
$z_{+}\nu_{+}+z_{-}\nu_{-}=0$ more directly, and will produce results more
obviously consistent with Newman [31]. The products of the simple association
equilibria for this system are unique, so the salt–charge basis is unique up
to ordering of the ions. We let the salt–charge transformation Z take the form

$\textbf{Z}=\begin{bmatrix}1&0&0\\\ 0&-\frac{\nu z_{-}}{z_{+}-z_{-}}&\frac{\nu
z_{+}}{z_{+}-z_{-}}\\\
0&\frac{z_{+}}{\left\|\textbf{z}\right\|}&\frac{z_{-}}{\left\|\textbf{z}\right\|}\end{bmatrix},$
(16.2)

in which $\left\|\textbf{z}\right\|=\sqrt{z_{+}^{2}+z_{-}^{2}}$ because
$z_{0}=0$. Given this, one can write the set of species concentrations in the
neutral state, $\textbf{c}^{0}$, as a parametric function of the component
concentrations $c_{0}$ and $c_{\textrm{e}}$ through equation (8.1),

$\left[\begin{array}[]{c}c_{0}^{0}\\\ c_{+}^{0}\\\
c_{-}^{0}\end{array}\right]=\textbf{Z}^{\top}\left[\begin{array}[]{c}c_{0}\\\
c_{\textrm{e}}\\\ 0\end{array}\right]=\left[\begin{array}[]{c}c_{0}\\\
-\frac{z_{-}\nu}{z_{+}-z_{-}}\cdot c_{\textrm{e}}\\\
\frac{z_{+}\nu}{z_{+}-z_{-}}\cdot c_{\textrm{e}}\end{array}\right].$ (16.3)

The electroneutral component chemical potentials can be written using
equations (8.24) and (8.3), which together with the relation $y_{0}=1-\nu
y_{\textrm{e}}$ show that

$\displaystyle\nabla\left[\begin{array}[]{c}\mu_{0}^{0}\\\
\mu_{\textrm{e}}^{0}\end{array}\right]$
$\displaystyle=RT\left[\begin{array}[]{cc}X_{00}^{0}&X_{0\textrm{e}}^{0}\\\
X_{\textrm{e}0}^{0}&X_{\textrm{ee}}^{0}\end{array}\right]\nabla\left[\begin{array}[]{c}y_{0}\\\
y_{\textrm{e}}\end{array}\right]$ (16.10)
$\displaystyle=\frac{RT\left(1+\Lambda_{00}\nu
y_{\textrm{e}}y_{0}\right)}{y_{\textrm{e}}y_{0}}\left[\begin{array}[]{c}-y_{\textrm{e}}\\\
y_{0}\end{array}\right]\nu\nabla y_{\textrm{e}},$ (16.13)

parametrized by the single excess activity-coefficient derivative function
$\Lambda_{00}\left(y_{\textrm{e}}\right)$.

Using the salt–charge transformation from equation (16.2) in equation (8.30)
shows that

$\textbf{M}_{\nu}=\frac{\nu
RT\left(z_{+}\mathscr{D}_{0+}-z_{-}\mathscr{D}_{0-}\right)}{c_{\textrm{T}}\left(z_{+}-z_{-}\right)\mathscr{D}_{0+}\mathscr{D}_{0-}}\left[\begin{array}[]{cc}\dfrac{c_{\textrm{e}}}{c_{0}}&-1\\\
-1&\dfrac{c_{0}}{c_{\textrm{e}}}\end{array}\right],$ (16.14)
$\textbf{m}_{z}=\frac{RT\left(z_{+}\mathscr{D}_{0-}+z_{-}\mathscr{D}_{0+}\right)}{c_{\textrm{T}}\left\|\textbf{z}\right\|^{2}\mathscr{D}_{0+}\mathscr{D}_{0-}}\left[\begin{array}[]{c}-1\\\
\dfrac{c_{0}}{c_{\textrm{e}}}\end{array}\right],~{}~{}~{}~{}\textrm{and}$
(16.15)
$M_{zz}=-\tfrac{RT\left\|\textbf{z}\right\|^{2}}{c_{\textrm{T}}z_{+}z_{-}}\Big{[}\tfrac{1}{\mathscr{D}_{+-}}+\tfrac{\left(z_{+}-z_{-}\right)\left(z_{+}^{3}\mathscr{D}_{0-}+z_{-}^{3}\mathscr{D}_{0+}\right)c_{0}}{\left\|\textbf{z}\right\|^{4}\mathscr{D}_{0+}\mathscr{D}_{0-}\nu
c_{\textrm{e}}}\Big{]}.$ (16.16)

Since the $\textbf{M}_{Z}^{0}$ matrix is invariant with respect to the choice
of convective velocity, these expressions can be used in equation (8.34) or
(10.1), or even alternative versions of equation (10.1) in which the set of
excess component fluxes relative to the mass-average velocity,
$\vec{\textbf{j}}_{Z}^{\,0}$, is replaced with a set relative to some other
velocity, $\vec{\textbf{j}}_{Z}^{\,0,\psi^{\prime}}$.

Concentrated solution theory uses flux-explicit forms of the transport laws.
Following Newman, we refer the excess fluxes to the solvent velocity (i.e.,
$\vec{N_{0}}/c_{0}$) by taking

$\bm{\psi}^{\prime}_{Z}=\frac{\textbf{Z}}{c_{0}^{0}}\left[\begin{array}[]{c}1\\\
0\\\ 0\end{array}\right]=\frac{1}{c_{0}}\left[\begin{array}[]{c}1\\\ 0\\\
0\end{array}\right]$ (16.17)

in equation (12.3). This changes the dynamical variables to a set of
electroneutral excess species fluxes

$\vec{\textbf{j}}_{Z}^{\,0,0}=\left[\begin{array}[]{c}0\\\
\vec{J}_{\textrm{e}}^{\,0,0}\\\
\dfrac{\vec{i}}{F\left\|\textbf{z}\right\|}\end{array}\right]$ (16.18)

whose first member is trivial, and whose second describes the excess salt flux
relative to the solvent velocity.

Since matrix $\textbf{M}_{Z}^{0}$ is built up of the sub-blocks defined by
equations (16.14) through (16.16), the transport matrix $\textbf{L}_{Z}^{0,0}$
for the flux-explicit constitutive formulation relative to the solvent
velocity can be computed directly with equation (12.6), using
$\bm{\psi}^{\prime}_{Z}$ from equation (16.17). The ionic conductivity is
identified as the entry in the last row and column of $\textbf{L}_{Z}^{0,0}$
through equation (12.9), yielding

$\kappa^{0}=-\frac{c_{\textrm{T}}F^{2}z_{+}z_{-}}{RT}\left[\frac{1}{\mathscr{D}_{+-}}+\frac{c_{0}\left(z_{+}-z_{-}\right)}{\nu
c_{\textrm{e}}\left(z_{+}\mathscr{D}_{0+}-z_{-}\mathscr{D}_{0-}\right)}\right]^{-1},$
(16.19)

whence the last column (or last row) of $\textbf{L}_{Z}^{0,0}$ reveals that

$\bm{\xi}^{0}=\left[\begin{array}[]{c}\xi_{0}^{0}\\\
\xi_{\textrm{e}}^{0}\end{array}\right]=\left[\begin{array}[]{c}0\\\
-\dfrac{\left(z_{+}-z_{-}\right)\left(z_{-}\mathscr{D}_{0+}+z_{+}\mathscr{D}_{0-}\right)}{\left\|\textbf{z}\right\|\nu\left(z_{+}\mathscr{D}_{0+}-z_{-}\mathscr{D}_{0-}\right)}\end{array}\right]$
(16.20)

are the component migration coefficients relative to the solvent velocity.
Migration coefficients can be used to compute transference numbers with
equation (13.2), yielding

$\textbf{t}^{0}=\left[\begin{array}[]{c}t_{0}^{0}\\\ t_{+}^{0}\\\
t_{-}^{0}\end{array}\right]=\left[\begin{array}[]{c}0\\\
\frac{z_{+}\mathscr{D}_{0+}}{z_{+}\mathscr{D}_{0+}-z_{-}\mathscr{D}_{0-}}\\\
-\frac{z_{-}\mathscr{D}_{0-}}{z_{+}\mathscr{D}_{0+}-z_{-}\mathscr{D}_{0-}}\end{array}\right].$
(16.21)

The transference number $t_{0}^{0}$ vanishes necessarily because the solvent
species is uncharged, conveniently matching the solvent component’s null
migration coefficient. Thus, for binary electrolytes, no information is lost
by using transference numbers.

A third flux-explicit transport property, $\mathfrak{D}$, is related to the
Stefan–Maxwell coefficients by comparing the computed expression (16.14) with
equation (14.1). This shows that

$\mathfrak{D}=\frac{\left(z_{+}-z_{-}\right)\mathscr{D}_{0+}\mathscr{D}_{0-}}{z_{+}\mathscr{D}_{0+}-z_{-}\mathscr{D}_{0-}}$
(16.22)

defines the thermodynamic diffusivity of salt in solvent. Given
$\textbf{M}_{\nu}$ from equation (16.14) and that
$\bm{\psi}^{\prime}_{\nu}=[1/c_{0}~{}0]^{\top}$, Onsager component
diffusivities relative to the solvent velocity,
$\overline{\textbf{L}}_{\nu}^{0}$, are given by equation (14.2), producing

$\overline{\textbf{L}}_{\nu}^{0}=\left[\begin{array}[]{cc}0&0\\\
0&\dfrac{c_{\textrm{e}}\left(z_{+}-z_{-}\right)\mathscr{D}_{0+}\mathscr{D}_{0-}}{\nu
c_{0}\left(z_{+}\mathscr{D}_{0+}-z_{-}\mathscr{D}_{0-}\right)}.\end{array}\right]$
(16.23)

Thus $\overline{\textbf{L}}_{\nu}^{0}$ is seen to be proportional to
$\mathfrak{D}$, through a composition dependent prefactor.

Equations (16.19), (16.20), and (16.22) establish three macroscopic transport
properties that parametrize equations (12.16) and (12.17) in terms of
Stefan–Maxwell diffusivities. Notably, equations (16.19), (16.21), and (16.22)
respectively match the definition of ionic conductivity in Newman’s equation
12.23, cation transference number relative to the solvent velocity in Newman’s
12.11, and thermodynamic diffusivity in Newman’s 12.10 [1].

Identifying flux-explicit species laws and a MacInnes equation involving a
single thermodynamic factor, with properties $\kappa^{0}$, $t_{+}^{0}$, and
$\mathfrak{D}$ defined in terms of Stefan–Maxwell coefficients, completes
Newman’s development of concentrated solution theory. The structural
properties laid out above make possible an opposite process, however, in which
one writes flux laws a priori with the structure that kinematic and
thermodynamic consistency necessitates. Once such laws have been formulated,
they can be analyzed to find how the Onsager–Stefan–Maxwell matrix
$\textbf{M}_{Z}^{0}$ depends on the posited set of flux-explicit properties.
This route more naturally reflects how properties are typically measured, by
fitting data to formulas derived from flux-explicit laws. It also develops the
relationships needed to incorporate standard properties into numerical solvers
that exploit the spectral structure of force-explicit transport formulations
[8, 18].

For the binary electrolyte, $\textbf{i}_{0}$ (the first column of the identity
matrix) sits in the nullspace of the symmetric Onsager component-diffusivity
matrix relative to the solvent velocity, $\overline{\textbf{L}}_{\nu}^{0}$,
which therefore has one nonzero entry (in the second row and second column),

$\overline{\textbf{L}}_{\nu}^{0}=\left[\begin{array}[]{cc}0&0\\\
0&\mathscr{L}\end{array}\right]=\left[\begin{array}[]{cc}0&0\\\
0&\dfrac{c_{\textrm{e}}\mathfrak{D}}{\nu c_{0}}\end{array}\right],$ (16.24)

in which the single Onsager diffusivity of salt relative to the solvent
velocity, $\mathscr{L}$, is expressed in terms of a thermodynamic diffusivity
$\mathfrak{D}$ by applying inversion formula (14.2) to definition (14.1).
Similarly, $\textbf{i}_{0}$ is orthogonal to $\bm{\xi}^{0}$, which therefore
also must have just one nonzero entry, that is,

$\bm{\xi}^{0}=\left[\begin{array}[]{c}0\\\ \xi,\end{array}\right]$ (16.25)

where $\xi$ stands for the salt migration coefficient relative to the solvent
velocity. This property relates to the cation transference number relative to
solvent as

$\xi=-\frac{\left(z_{+}-z_{-}\right)\left[\big{(}1-t_{+}^{0}\big{)}z_{+}^{2}-z_{-}^{2}t_{+}^{0}\right]}{\nu
z_{+}z_{-}\left\|\textbf{z}\right\|}$ (16.26)

through equation (13.2).

The known structures of $\overline{\textbf{L}}_{\nu}^{0}$ and $\bm{\xi}^{0}$
allow species flux laws relative to the solvent velocity to be written
directly in the form of equation (12.16). Bringing in equations (16.1),
(16.21), and (16.24), one finds

$\displaystyle\vec{J}_{0}^{\,0}$ $\displaystyle=0,$ (16.27)
$\displaystyle\vec{J}_{+}^{\,0}$
$\displaystyle=-\frac{\mathfrak{D}c_{\textrm{T}}^{0}}{RT}\frac{\nu_{+}c_{\textrm{e}}}{\nu
c_{0}}\nabla\mu_{\textrm{e}}^{0}+\frac{t_{+}^{0}\vec{i}}{Fz_{+}},$ (16.28)
$\displaystyle\vec{J}_{-}^{\,0}$
$\displaystyle=-\frac{\mathfrak{D}c_{\textrm{T}}^{0}}{RT}\frac{\nu_{-}c_{\textrm{e}}}{\nu
c_{0}}\nabla\mu_{\textrm{e}}^{0}+\frac{\big{(}1-t_{+}^{0}\big{)}\vec{i}}{Fz_{-}}.$
(16.29)

Note that equations (16.28) and (16.29) match Newman’s equations 12.8 and
12.9, respectively [1].

Finally, posit an ionic conductivity $\kappa^{0}$ to write an electroneutral
MacInnes equation in the form of (12.17), as

$\frac{\vec{i}}{\kappa^{0}}=-\nabla\Phi_{z}^{0}-\frac{\left(z_{+}-z_{-}\right)}{\nu
Fz_{+}z_{-}}\left(\frac{z_{+}^{2}}{z_{+}^{2}+z_{-}^{2}}-t_{+}^{0}\right)\nabla\mu_{\textrm{e}}^{0},$
(16.30)

which depends on the salt–charge potential. Introducing a reference-electrode
potential with equation (15.6) gives

$\displaystyle\frac{\vec{i}}{\kappa^{0}}=$
$\displaystyle-\nabla\Phi^{\ominus}+\frac{s_{0}\nabla\mu_{0}^{0}}{n_{\textrm{e}^{-}}Fz_{\textrm{e}^{-}}}$
$\displaystyle+\frac{z_{+}-z_{-}}{n_{\textrm{e}^{-}}Fz_{\textrm{e}^{-}}\nu}\left[t_{+}^{0}\frac{s_{-}}{z_{+}}-\big{(}1-t_{+}^{0}\big{)}\frac{s_{+}}{z_{-}}\right]\nabla\mu_{\textrm{e}}^{0},$
(16.31)

in harmony with Newman’s equation 12.27 [1].

Given $\kappa^{0}$, as well as $\overline{\textbf{L}}_{\nu}^{0}$ and
$\bm{\xi}^{0}$ matrices in the forms from equations (16.24) and (16.25), one
can use relationships (12.18) and (12.19) to construct the
$\textbf{M}_{Z}^{0}$ matrix that would be needed to implement the
electroneutral Onsager–Stefan–Maxwell equations. Performing the calculation,
one finds

$\textbf{M}_{\nu}=\frac{RT}{c_{\textrm{T}}^{0}\mathfrak{D}}\left[\begin{array}[]{cc}\frac{\nu
c_{\textrm{e}}}{c_{0}}&-\nu\\\ -\nu&\frac{\nu
c_{0}}{c_{\textrm{e}}}\end{array}\right],$ (16.32) $\textbf{m}_{z}=\frac{\nu
RT\xi}{c_{\textrm{T}}^{0}\mathfrak{D}}\left[\begin{array}[]{c}1\\\
-\tfrac{c_{0}}{c_{\textrm{e}}}\end{array}\right],$ (16.33)

and

$M_{zz}=\frac{F^{2}\left\|\textbf{z}\right\|^{2}}{\kappa^{0}}+\frac{\nu
RTc_{0}\xi^{2}}{c_{\textrm{T}}^{0}c_{\textrm{e}}\mathfrak{D}}.$ (16.34)

One can also elucidate the Stefan–Maxwell coefficients that describe species
interdiffusion, by inverting the congruence transformation in equation (8.30):
$\textbf{M}^{0}=\textbf{Z}^{-1}\textbf{M}_{Z}^{0}\textbf{Z}^{-\top}$. Off-
diagonal entries of $\textbf{M}^{0}$ relate to $\mathscr{D}_{ij}$ values
through equation (1.2). Once $t_{+}^{0}$ has been introduced in favor of $\xi$
via equation (16.26), one finds

$\displaystyle\frac{1}{\mathscr{D}_{0+}}=$
$\displaystyle-\frac{(1-t_{+}^{0})\left(z_{+}-z_{-}\right)}{z_{-}\mathfrak{D}},$
$\displaystyle\frac{1}{\mathscr{D}_{0-}}=$
$\displaystyle\frac{t_{+}^{0}\left(z_{+}-z_{-}\right)}{z_{+}\mathfrak{D}},$
$\displaystyle\frac{1}{\mathscr{D}_{+-}}=$
$\displaystyle-\frac{F^{2}z_{+}z_{-}c_{\textrm{T}}^{0}}{RT\kappa^{0}}+\frac{t_{+}^{0}(1-t_{+}^{0})\left(z_{+}-z_{-}\right)^{2}c_{0}}{\nu\mathfrak{D}z_{+}z_{-}c_{\textrm{e}}},$
(16.35)

consistent with Newman’s equations 14.1–14.3.

###  16.2 The molten salt

Pollard and Newman applied concentrated solution theory to the other possible
liquid electrolyte comprising three species: a molten solution of two simple
salts that share a common ion [30]. From transport laws based on the
Onsager–Stefan–Maxwell equations, they ultimately developed flux-explicit laws
in the form of equation (12.16) and a MacInnes equation like equation (12.17).

Here we apply the salt–charge framework to obtain transport laws consistent
with Pollard and Newman’s results. The flux laws we formulate will not
precisely identify with the Pollard–Newman equations, however. Whereas they
used the common ion’s velocity as the reference for convection, we instead use
the velocity of one of the salts. This alternative convention creates
component flux laws that involve a single chemical-potential gradient
directly, rather than requiring additional use of the Gibbs–Duhem equation to
eliminate a dependent driving force.

Let the charges of the three species comprising the binary molten salt be
$\textbf{z}=[z_{1},z_{2},z_{3}]^{\top}$. We assume that the third species is
counter-charged to the first two, so that
$\textrm{sgn}\left(z_{1}\right)=\textrm{sgn}\left(z_{2}\right)$ and
$\textrm{sgn}\left(z_{1}\right)=-\textrm{sgn}\left(z_{3}\right)$. Again, the
choice of products of the simple association equilibria is unique; we label
the simple salts formed from these three species with indices corresponding to
the similarly charged ion that distinguishes them: $\text{e}_{1}$ and
$\text{e}_{2}$. The stoichiometric coefficients of the ions within these salts
are written in terms of the species charges and total formula-unit
stoichiometries $\nu_{\text{e}_{1}}$ and $\nu_{\text{e}_{2}}$ as

$\displaystyle\nu_{1}^{\text{e}_{1}}=-\frac{\nu_{\text{e}_{1}}z_{3}}{z_{1}-z_{3}},\>\>\>\>\>\>\nu_{3}^{\text{e}_{1}}=\frac{\nu_{\text{e}_{1}}z_{1}}{z_{1}-z_{3}},$
(16.36)
$\displaystyle\nu_{2}^{\text{e}_{2}}=-\frac{\nu_{\text{e}_{2}}z_{3}}{z_{2}-z_{3}},\>\>\>\>\>\>\nu_{3}^{\text{e}_{2}}=\frac{\nu_{\text{e}_{2}}z_{2}}{z_{2}-z_{3}},$
(16.37)

so that $\nu_{\text{e}_{1}}=\nu_{1}^{\text{e}_{1}}+\nu_{3}^{\text{e}_{1}}$ and
$\nu_{\text{e}_{2}}=\nu_{2}^{\text{e}_{2}}+\nu_{3}^{\text{e}_{2}}$. The matrix
Z that represents the salt–charge basis is unique up to ordering of the
similarly charged ions. We choose to write it as

$\textbf{Z}=\begin{bmatrix}\nu_{1}^{\text{e}_{1}}&0&\nu_{3}^{\text{e}_{1}}\\\
0&\nu_{2}^{\text{e}_{2}}&\nu_{3}^{\text{e}_{2}}\\\
\frac{z_{1}}{\|\textbf{z}\|}&\frac{z_{2}}{\|\textbf{z}\|}&\frac{z_{3}}{\|\textbf{z}\|}\end{bmatrix},$
(16.38)

where $\left\|\textbf{z}\right\|=\sqrt{z_{1}^{2}+z_{2}^{2}+z_{3}^{2}}$,
bearing in mind that the reactant stoichiometric coefficients from the simple
association equilibria satisfy equations (16.36) and (16.37).

With the salt–charge transformation from equation (16.38), the electroneutral
species concentrations $\textbf{c}^{0}$ are cast as parametric functions of
the component concentrations $c_{\textrm{e}_{1}}$ and $c_{\textrm{e}_{2}}$
through equation (8.1), as

$\left[\begin{array}[]{c}c_{1}^{0}\\\ c_{2}^{0}\\\
c_{3}^{0}\end{array}\right]=\textbf{Z}^{\top}\left[\begin{array}[]{c}c_{\textrm{e}_{1}}\\\
c_{\textrm{e}_{2}}\\\
0\end{array}\right]=\left[\begin{array}[]{c}-\frac{\nu_{\text{e}_{1}}z_{3}}{z_{1}-z_{3}}c_{\text{e}_{1}}\\\
-\frac{\nu_{\text{e}_{2}}z_{3}}{z_{2}-z_{3}}c_{\text{e}_{2}}\\\
\nu_{1}^{\textrm{e}_{1}}c_{\text{e}_{1}}+\nu_{2}^{\textrm{e}_{2}}c_{\text{e}_{2}}\end{array}\right].$
(16.39)

The electroneutral component chemical potentials can then be written using
equations (8.24) and (8.3). Together with the relation
$\nu_{\textrm{e}_{2}}y_{\textrm{e}_{2}}=1-\nu_{\textrm{e}_{1}}y_{\textrm{e}_{1}}$,
these show that

$\displaystyle\vec{\nabla}\left[\begin{array}[]{c}\mu_{\textrm{e}_{1}}^{0}\\\
\mu_{\textrm{e}_{2}}^{0}\end{array}\right]=$
$\displaystyle\,RT\left[\begin{array}[]{cc}X_{\textrm{e}_{1}\textrm{e}_{1}}^{0}&X_{\textrm{e}_{1}\textrm{e}_{2}}^{0}\\\
X_{\textrm{e}_{2}\textrm{e}_{2}}^{0}&X_{\textrm{e}_{2}\textrm{e}_{2}}^{0}\end{array}\right]\vec{\nabla}\left[\begin{array}[]{c}y_{\textrm{e}_{1}}\\\
y_{\textrm{e}_{2}}\end{array}\right]$ (16.46) $\displaystyle=$
$\displaystyle\,RT\left[\tfrac{z_{3}\left(z_{1}\nu_{\textrm{e}_{1}}y_{\textrm{e}_{1}}+z_{2}\nu_{\textrm{e}_{2}}y_{\textrm{e}_{2}}\right)}{z_{1}\left(z_{3}-z_{2}\right)\nu_{\textrm{e}_{1}}y_{\textrm{e}_{1}}+z_{2}\left(z_{3}-z_{1}\right)\nu_{\textrm{e}_{2}}y_{\textrm{e}_{2}}}\right.$
$\displaystyle\left.+\tfrac{\nu_{\textrm{e}_{1}}y_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}y_{\textrm{e}_{2}}\Lambda^{0}_{\textrm{e}_{1}\textrm{e}_{1}}}{\nu_{\textrm{e}_{1}}^{2}}\right]\left[\begin{array}[]{c}-\frac{1}{y_{\textrm{e}_{1}}}\\\
\frac{1}{y_{\textrm{e}_{2}}}\end{array}\right]\nu_{\textrm{e}_{2}}\nabla
y_{\textrm{e}_{2}},$ (16.49)

in which a single function
$\Lambda_{\textrm{e}_{1}\textrm{e}_{1}}\left(y_{\textrm{e}_{2}}\right)$
parametrizes the composition gradient of either salt’s mean molar activity
coefficient within the solution. Observe that the ideal contribution to mixing
free energy—the first term in the square brackets—depends in a rather subtle
way on the salt fractions $y_{\textrm{e}_{1}}$ and $y_{\textrm{e}_{2}}$.

Rather than repeating in its entirety the development illustrated for the
binary electrolyte, we leverage the structural knowledge provided by our
framework to abbreviate the process. Because the solution comprises three
species, there are three flux-explicit transport properties $\mathfrak{D}$,
$\xi$ and $\kappa^{0}$. Using the structure posited in equation (14.1), we may
write

$\textbf{M}_{\nu}=\frac{RT\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}}{c_{\text{T}}^{0}\mathfrak{D}}\left[\begin{array}[]{cc}\dfrac{c_{\textrm{e}_{2}}}{c_{\textrm{e}_{1}}}&-1\\\
-1&\dfrac{c_{\textrm{e}_{1}}}{c_{\textrm{e}_{2}}}\end{array}\right].$ (16.50)

In order to eliminate the chemical-potential gradient of salt 1 from the flux-
explicit laws, it is convenient to choose $\bm{\psi}^{\prime}$ such that
$\bm{\psi}^{\prime}_{Z}=\textbf{i}_{\textrm{e}_{1}}/c_{\textrm{e}_{1}}$, that
is, to let the convective velocity be that of the first salt (i.e.,
$\vec{N}_{\textrm{e}_{1}}/c_{\textrm{e}_{1}}$). There is a single non-trivial
migration coefficient $\xi$ relative to this reference velocity:

$\bm{\xi}^{\textrm{e}_{1}}=\left[\begin{array}[]{c}0\\\
\xi\end{array}\right].$ (16.51)

Using equations (12.19) then shows that

$\textbf{m}_{z}=\frac{RT\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}\xi}{c_{\text{T}}^{0}\mathfrak{D}}\left[\begin{array}[]{c}-1\\\
\dfrac{c_{\textrm{e}_{1}}}{c_{\textrm{e}_{2}}}\end{array}\right]~{}~{}~{}~{}\textrm{and}$
(16.52)
$M_{zz}=\frac{F^{2}\|z\|^{2}}{\kappa^{0}}+\frac{RT\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}\xi^{2}}{c_{\text{T}}^{0}\mathfrak{D}}\dfrac{c_{\textrm{e}_{1}}}{c_{\textrm{e}_{2}}}.$
(16.53)

Next, taking the transport matrix (1.2), performing the congruence
transformation (5.7), and comparing the result to equations (16.50), (16.52),
and (16.53) allows identification of the three macroscopic transport
properties in terms of Stefan–Maxwell coefficients. These are

$\displaystyle\mathfrak{D}=$
$\displaystyle\frac{c^{0}_{3}(z_{1}-z_{3})(z_{2}-z_{3})}{\left(\frac{c^{0}_{3}z_{3}^{2}}{\mathscr{D}_{12}}+\frac{c^{0}_{2}z_{2}^{2}}{\mathscr{D}_{13}}+\frac{c^{0}_{1}z_{1}^{2}}{\mathscr{D}_{23}}\right)},$
(16.54)
$\displaystyle\xi\cdot\tfrac{\nu_{\text{e}_{1}}\nu_{\text{e}_{2}}c^{0}_{1}c^{0}_{3}z_{3}^{2}}{\mathfrak{D}\nu_{1}^{\textrm{e}_{1}}}=$
$\displaystyle\tfrac{z_{3}^{2}c_{3}(c^{0}_{1}z_{2}-c^{0}_{2}z_{1})}{\mathscr{D}_{12}\|z\|^{2}}-\tfrac{z_{2}c^{0}_{1}c^{0}_{2}(z_{1}^{2}+z_{3}^{2})+({c^{0}_{2}})^{2}z_{1}z_{2}^{2}}{\mathscr{D}_{13}\|z\|^{2}}$
$\displaystyle+\tfrac{({c^{0}_{1}})^{2}z_{1}^{2}z_{2}+c^{0}_{1}c^{0}_{2}z_{1}(z^{2}_{2}+z_{3}^{2})}{\mathscr{D}_{23}\|z\|^{2}},$
(16.55) $\displaystyle\kappa^{0}=$
$\displaystyle\frac{F^{2}c_{\textrm{T}}^{0}\left(\frac{z_{1}^{2}c_{1}^{0}}{\mathscr{D}_{23}}+\frac{z_{2}^{2}c_{2}^{0}}{\mathscr{D}_{13}}+\frac{z_{3}^{2}c_{3}^{0}}{\mathscr{D}_{12}}\right)}{RT\left(\frac{c_{1}^{0}}{\mathscr{D}_{12}\mathscr{D}_{13}}+\frac{c_{2}^{0}}{\mathscr{D}_{12}\mathscr{D}_{23}}+\frac{c_{3}^{0}}{\mathscr{D}_{13}\mathscr{D}_{23}}\right)}.$
(16.56)

Here the conductivity is identical to Pollard and Newman’s equation 33 [30].
The diffusivity defined in equation (16.54) relates to the diffusion
coefficient identified in their equation 9 [30], which we will label
$\mathscr{D}$, as

$\mathfrak{D}=\mathscr{D}\frac{c_{\text{T}}^{0}}{c_{3}^{0}}\frac{\nu_{\text{e}_{1}}\nu_{\text{e}_{2}}}{\nu^{\text{e}_{1}}_{1}\nu^{\text{e}_{2}}_{1}}.$
(16.57)

Although the prefactor that relates our $\mathfrak{D}$ to Pollard and Newman’s
$\mathscr{D}$ is composition-dependent, the factor only involves the particle
fraction of the common ion, which is always of order unity.

It is equally straightforward to invert this process. Similar to the binary-
electrolyte analysis, entries of $\textbf{M}^{0}$ may be related to the
original Stefan–Maxwell diffusivities $\mathscr{D}_{ij}$ by inverting the
congruence transformation in equation (8.30) and examining the off-diagonal
elements in light of equation (8.29). This comparison shows that the
diffusivities relate to the macroscopic properties as

$\displaystyle\tfrac{1}{\mathscr{D}_{12}}=$
$\displaystyle\tfrac{\left[c_{1}^{0}(z_{1}^{2}+z_{3}^{2})+c_{2}^{0}z_{1}z_{2}\right]\left[c_{2}^{0}(z_{2}^{2}+z_{3}^{2})+c_{1}^{0}z_{1}z_{2}\right]\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}}{\mathfrak{D}c_{1}^{0}c_{2}^{0}\nu_{1}^{\textrm{e}_{1}}\nu_{2}^{\textrm{e}_{2}}\|\textbf{z}\|^{4}}$
$\displaystyle-\tfrac{\xi\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}[c_{1}^{0}z_{1}(z_{1}^{2}-z_{2}^{2}+z_{3}^{2})-c_{2}^{0}z_{2}(-z_{1}^{2}+z_{2}^{2}+z_{3}^{2})]}{\mathfrak{D}\|\textbf{z}\|^{3}c_{2}^{0}\nu^{\textrm{e}_{1}}_{1}}$
$\displaystyle-\left(\tfrac{F^{2}c_{\text{T}}^{0}}{RT\kappa^{0}}+\tfrac{\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}\xi^{2}}{\mathfrak{D}}\tfrac{c_{\textrm{e}_{1}}}{c_{\textrm{e}_{2}}}\right)\tfrac{z_{1}z_{2}}{\|\textbf{z}\|^{2}}$
$\displaystyle\tfrac{1}{\mathscr{D}_{13}}=$
$\displaystyle-\tfrac{z_{3}(c_{1}^{0}z_{2}-c_{2}^{0}z_{1})\left[c_{2}^{0}(z_{2}^{2}+z_{3}^{2})+c_{1}^{0}z_{1}z_{2}\right]\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}}{\mathfrak{D}c_{1}^{0}c_{2}^{0}\nu_{1}^{\textrm{e}_{1}}\nu_{2}^{\textrm{e}_{2}}\|\textbf{z}\|^{4}}$
$\displaystyle+\tfrac{\xi
z_{3}\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}\left[c_{2}^{0}(-z_{1}^{2}+z_{2}^{2}+z_{3}^{2})+2c_{1}^{0}z_{1}z_{2}\right]}{\mathfrak{D}c_{2}^{0}\nu_{1}^{\textrm{e}_{1}}\|\textbf{z}\|^{3}}$
$\displaystyle-\left(\tfrac{F^{2}c_{\text{T}}^{0}}{RT\kappa^{0}}+\tfrac{\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}\xi^{2}}{\mathfrak{D}}\tfrac{c_{\textrm{e}_{1}}}{c_{\textrm{e}_{2}}}\right)\tfrac{z_{1}z_{3}}{\|\textbf{z}\|^{2}}$
$\displaystyle\tfrac{1}{\mathscr{D}_{23}}=$
$\displaystyle-\tfrac{z_{3}(c_{2}^{0}z_{1}-c_{1}^{0}z_{2})\left[c_{2}^{0}(z_{1}^{2}+z_{3}^{2})+c_{2}^{0}z_{1}z_{2}\right]\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}}{\mathfrak{D}c_{1}^{0}c_{2}^{0}\nu_{1}^{\textrm{e}_{1}}\nu_{2}^{\textrm{e}_{2}}\|\textbf{z}\|^{4}}$
$\displaystyle+\tfrac{\xi
z_{3}\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}\left[c_{1}^{0}(z_{1}^{2}-z_{2}^{2}+z_{3}^{2})+2c_{2}^{0}z_{1}z_{2}\right]}{\mathfrak{D}c_{2}^{0}\nu_{1}^{\textrm{e}_{1}}\|\textbf{z}\|^{3}}$
$\displaystyle-\left(\tfrac{F^{2}c_{\text{T}}^{0}}{RT\kappa^{0}}+\tfrac{\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}\xi^{2}}{\mathfrak{D}}\tfrac{c_{\textrm{e}_{1}}}{c_{\textrm{e}_{2}}}\right)\tfrac{z_{2}z_{3}}{\|\textbf{z}\|^{2}}.$
(16.58)

Pollard and Newman do not write down these relationships, presumably because
of their algebraic complexity. Stating explicit forms for the Stefan–Maxwell
coefficients in terms of flux-explicit properties becomes increasingly
cumbersome as the number of species (and the number of charged species) in a
system increases.

With our choice of convective velocity, the matrix
$\overline{\textbf{L}}_{\nu}^{\,\textrm{e}_{1}}$ also affords a single non-
trivial entry,

$\overline{\textbf{L}}_{\nu}^{\,\textrm{e}_{1}}=\left[\begin{array}[]{cc}0&0\\\
0&\frac{\mathfrak{D}}{\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}}\frac{c_{\textrm{e}_{2}}}{c_{\textrm{e}_{1}}}\end{array}\right].$
(16.59)

After deriving transference numbers $t_{1}^{\textrm{e}_{1}}$,
$t_{2}^{\textrm{e}_{1}}$, and $t_{3}^{\textrm{e}_{1}}$ relative to salt 1 with
equation (13.2), the laws governing the species excess fluxes are found to be

$\displaystyle\vec{J}_{1}^{\,\textrm{e}_{1}}$
$\displaystyle=\frac{t_{1}^{\textrm{e}_{1}}\vec{i}}{Fz_{1}},$ (16.60)
$\displaystyle\vec{J}_{2}^{\,\textrm{e}_{1}}$
$\displaystyle=-\frac{c_{\textrm{T}}^{0}\mathfrak{D}}{RT}\frac{\nu^{\textrm{e}_{1}}_{1}c_{\textrm{e}_{2}}}{\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}c_{\textrm{e}_{1}}}\nabla\mu_{\textrm{e}_{2}}^{0}+\frac{t_{2}^{\textrm{e}_{1}}\vec{i}}{Fz_{2}},$
(16.61) $\displaystyle\vec{J}_{3}^{\,\textrm{e}_{1}}$
$\displaystyle=-\frac{c_{\textrm{T}}^{0}\mathfrak{D}}{RT}\frac{\nu^{\textrm{e}_{2}}_{2}c_{\textrm{e}_{2}}}{\nu_{\textrm{e}_{1}}\nu_{\textrm{e}_{2}}c_{\textrm{e}_{1}}}\nabla\mu_{\textrm{e}_{2}}^{0}+\frac{t_{3}^{\textrm{e}_{1}}\vec{i}}{Fz_{3}}.$
(16.62)

Interestingly, the excess flux of ion 1 relative to the velocity of salt 1
does not vanish; it can be driven by migration, but not by diffusion. Equation
(13.2) shows that all three transference numbers here are parametric functions
of the single migration coefficient $\xi$, through

$\left[\begin{array}[]{c}t_{1}^{\textrm{e}_{1}}\\\ t_{2}^{\textrm{e}_{1}}\\\
t_{3}^{\textrm{e}_{1}}\end{array}\right]=\frac{\operatorname{diag}\left(\textbf{z}\right)\textbf{Z}^{\top}}{\left\|\textbf{z}\right\|}\left[\begin{array}[]{c}0\\\
\xi\\\ 1\end{array}\right].$ (16.63)

Thus, as well as satisfying the relation
$t_{1}^{\textrm{e}_{1}}+t_{2}^{\textrm{e}_{1}}+t_{3}^{\textrm{e}_{1}}=1$ that
follows from the last row of this matrix equation, the species transference
numbers relative to salt 1 are also constrained such that
$(z_{2}^{2}+z_{3}^{2})t_{1}^{\textrm{e}_{1}}+z_{1}^{2}(t_{2}^{\textrm{e}_{1}}+t_{3}^{\textrm{e}_{1}})=0$,
which follows from the first row.

The reference velocity under $\bm{\xi}^{\textrm{e}_{1}}$ can be changed to
that of the common ion by moving through the mass-average velocity with
transformation (12.15). Then one can derive the associated transference
numbers, denoted $t_{i}^{3}$, with equation (13.2) to find Pollard and
Newman’s expression

$t_{2}^{3}=\frac{\frac{z_{2}}{\mathscr{D}_{13}}-\frac{z_{3}}{\mathscr{D}_{12}}}{\left(\frac{z_{2}}{\mathscr{D}_{13}}-\frac{z_{3}}{\mathscr{D}_{12}}\right)+\frac{z_{1}c^{0}_{1}}{z_{2}c^{0}_{2}}\left(\frac{z_{1}}{\mathscr{D}_{23}}-\frac{z_{3}}{\mathscr{D}_{12}}\right)},$
(16.64)

which appears in their equation 10 [30]. With a migration coefficient relative