Multicomponent Maxwell–Stefan Diffusion

• Multicomponent Maxwell–Stefan diffusion is a continuum model characterizing cross-diffusion in non-dilute mixtures by rigorously incorporating interspecies friction.
• It employs a nonlinear flux–gradient relation derived from non-equilibrium thermodynamics, resulting in a strongly coupled parabolic PDE system.
• Analytical and numerical methods use Perron–Frobenius theory and entropy dissipation to ensure physical constraints and exponential decay to equilibrium.

Multicomponent Maxwell–Stefan (MS) diffusion is the fundamental continuum model for cross-diffusion in non-dilute mixtures, especially gases, that rigorously accounts for interspecies friction and the resulting nonlinear coupling of diffusive fluxes. Unlike Fickian diffusion, which prescribes uncoupled flux–gradient relations, the Maxwell–Stefan formalism originates from non-equilibrium thermodynamics and statistical physics, producing constrained systems where the constitutive relation for each species involves the relative velocities of all others. The resulting PDE system forms a nonlinear, cross-diffusive, generally non-diagonal, strongly coupled parabolic evolution, which is essential for capturing phenomena such as reverse (uphill) diffusion and ensuring thermodynamic consistency. Analytical studies focus on entropy methods and matrix Perron–Frobenius theory to resolve the degeneracy and quasi-positivity of the MS diffusion operator, while numerical simulations must invert nonlinear, singular cross-diffusion matrices at each step. Recent developments establish global-in-time existence and exponential decay to equilibrium for weak solutions in physically relevant geometries under isothermal, isobaric constraints (Jüngel et al., 2012).

1. Mathematical Formulation and Governing Equations

Let $\Omega\subset\mathbb{R}^d$ ($d\leq3$) be a bounded $C^{1,1}$ domain. Consider $N+1$ chemical species with molar concentrations $c_i(x,t)\geq0$, $i\in\{1,\dots,N+1\}$, normalized by $\sum_{i=1}^{N+1}c_i(x,t)=1$. The isothermal, isobaric multicomponent mass balances read

$$\partial_t c_i + \nabla\cdot J_i = r_i(c), \quad i=1,\dots,N+1,$$

for $(x,t)\in\Omega\times(0,\infty)$, where $J_i=c_i u_i$ is the molar diffusive flux, and $r_i(c)$ is a species production rate satisfying $\sum_i r_i(c)=0$ for conservation. Initial data and homogeneous Neumann boundary conditions,

$$\nabla c_i\cdot n=0 \quad\text{on }\partial\Omega,\,\, c_i(\cdot,0)=c_i^0(\cdot),$$

preserve total mass.

The flux–gradient relation is given by the Maxwell–Stefan system, derived from thermodynamic force balance: $$d_i = f_i, \quad d_i=-RT\nabla \ln c_i, \quad f_i=-\sum_{j\neq i} c_i c_j D_{ij} (u_i - u_j),$$ where $D_{ij}>0$ are the binary MS diffusivities ($D_{ij}=D_{ji}$).

In matrix form, this gives

$\nabla c = A(c) J,$

where $A(c)$ is the $(N+1)\times(N+1)$ quasi-positive, irreducible, concentration-dependent diffusion matrix: $$a_{ij} = d_{ij} c_i\, (i\neq j), \quad a_{ii}=-\sum_{j\neq i} d_{ij} c_j, \quad d_{ij}=1/D_{ij}>0, \ i\neq j.$$

The algebraic constraint $\sum_{i=1}^{N+1} J_i=0$ arises from the total mass balance and is encoded in the singularity of $A(c)$: $\ker A(c)=\mathrm{span}\{c\}$, $0\in\sigma(A)$ is a simple eigenvalue.

2. Diffusion Matrix Structure and Entropy Methods

The non-symmetric, singular, cross-diffusion operator $A(c)$ is quasi-positive and irreducible when all $c_i>0$. By Perron–Frobenius theory, its spectrum satisfies $\sigma(-A)\subset\{0\}\cup [\delta,\Lambda)$, $0<\delta\leq\Lambda<\infty$, and $A$ is invertible on its image. This structure is crucial for mathematical analysis: invertibility on the orthogonal complement of $\ker A$ ensures the solvability of the flux–gradient system up to the natural conservation constraints.

A key technique is the introduction of the entropy (Gibbs) functional: $$h(c) = \sum_{i=1}^{N+1} c_i(\ln c_i - 1), \qquad H[c]=\int_\Omega h(c)\,dx.$$ For nonnegative, mass-conserving reaction rates $r_i(c)$ with $\sum_i r_i(c)\ln c_i \leq 0$, one obtains the entropy–dissipation estimate: $$\frac{d}{dt} H[c] + K\sum_{i=1}^{N+1} \int_\Omega |\nabla \sqrt{c_i}|^2\,dx \leq 0,$$ where $K>0$ depends only on $\{D_{ij}\}$. The entropy is non-increasing and provides control over the system's long-time behavior and regularity.

By introducing entropy variables

$$w_i = \frac{\partial h}{\partial c_i} = \ln\frac{c_i}{c_{N+1}}, \quad i=1,\dots,N,$$

the PDE becomes, for $c'=(c_1,\dots,c_N)$,

$\partial_t c' - \nabla\cdot(B(w)\nabla w) = r'(c),$

with $B(w)=A_0(c)^{-1} H(c)^{-1}$ symmetric positive-definite, and the inverse change-of-variables mapping $w\mapsto c$ ensures positivity of concentrations for bounded $w$.

3. Existence Theory and Regularity of Solutions

Under the assumptions:

• $\Omega\subset\mathbb{R}^d$, $d\leq3$, $C^{1,1}$,
• $D_{ij}>0$ constants,
• $0\leq c_i^0\leq 1$ a.e., $\sum_i c_i^0=1$, $H[c^0]<\infty$,
• $r_i\in C^0([0,1]^{N+1})$, $\sum_i r_i=0$, $\sum_i r_i\ln c_i\leq0$,

the global-in-time existence of bounded weak solutions is established: for all $i$,

$$0\leq c_i\leq 1,\quad c_i\in L^2_{\text{loc}}(0,\infty; H^1(\Omega)),\quad \partial_t c_i\in L^2_{\text{loc}}(0,\infty; (H^2(\Omega))').$$

The PDE system holds in the weak sense, and the invariances ensure the physical constraints of non-negativity and conservation. Entropy variables yield a reformulation with symmetric, positive-definite operators, ensuring coercivity and allowing a priori $H^1$ bounds on $c_i$.

The analytical proof combines:

• The Perron–Frobenius spectral theory for quasi-positive matrices,
• The entropy–dissipation method for parabolic regularity,
• A time-discretization (implicit Euler), spatial regularization (fourth-order $H^2$ term), and passage of limits using compactness,
• Uniform $L^\infty$ and positivity bounds via the structure of the entropy variable transformation.

4. Reduction, Cross-diffusion, and Reaction Coupling

Reduction to $N$ components exploits $\sum_i c_i=1$ and $\sum_i J_i=0$, yielding an $N\times N$ reduced diffusion matrix $A_0(c')$. The general system is thus a quasilinear, strongly coupled, non-diagonal cross-diffusion PDE,

$$\partial_t c - \nabla\cdot\big(-A(c)^{-1} \nabla c\big)=r(c),$$

or, in entropy variables and reduced form, a symmetric cross-diffusion operator.

Cross-diffusion means that the flux of each species depends not only on its own concentration gradient but also on those of all other species, mediated by $A(c)^{-1}$. This leads to phenomena impossible in Fickian models, such as uphill or reverse diffusion.

5. Long-Time Asymptotics and Decay to Equilibrium

For conservative mixtures ($r_i\equiv 0$) and homogeneous Neumann boundary conditions, the unique steady state is spatially uniform: $$\bar c_i = \frac{1}{|\Omega|}\int_\Omega c_i^0(x)\,dx,\quad \sum\bar c_i=1.$$ Defining the relative entropy

$$H^*[c]=\int_\Omega\sum_{i=1}^{N+1}c_i\ln\frac{c_i}{\bar c_i}\,dx,$$

combining the entropy dissipation with logarithmic Sobolev and Csiszár–Kullback inequalities yields

$$H^*[c(t)] \leq e^{-\lambda t} H^*[c(0)], \quad \sum_{i=1}^{N+1}\|c_i(t)-\bar c_i\|_{L^1(\Omega)}\leq C e^{-\mu t},$$

for constants $\lambda,\mu>0$ depending only on the domain and diffusivities. Thus, solutions converge exponentially fast to the homogeneous equilibrium.

6. Implications for Modeling and Numerical Computation

The MS system strictly refines the classical Fickian approach, capturing multicomponent cross-effects and conforming to the constraints of nonequilibrium thermodynamics. The singularity and cross-diffusive structure demand careful treatment in analysis and numerics: inversion of $A(c)$ is only valid on the constrained subspace, and entropy structure must be preserved to maintain physicality (e.g., $0<c_i<1$).

Numerical algorithms must solve, at each time step and potentially grid point, a nonlinear algebraic system (often reduced to $N$ species) with