Onsager Diffusion v8.4

1 From Fick to Onsager: Multi-Component Diffusion

In a single-component system, diffusion is driven solely by the concentration gradient of that species. Fick's First Law describes the flux as proportional to this gradient:

J = -D \, \nabla C

Here, $D$ is the diffusion coefficient (units: m²/s), which measures how rapidly atoms move down their concentration gradient. The negative sign ensures that flux flows from high to low concentration.

What is the Diffusion Coefficient $D$?

The diffusion coefficient quantifies atomic mobility. It depends on temperature through the Arrhenius equation:

D(T) = D_0 \exp\left( -\frac{Q}{k_B T} \right)

$D_0$ (pre-exponential factor) captures the attempt frequency and lattice geometry.
$Q$ (activation energy) is the energy barrier for an atom to jump to a neighboring site.
$k_B$ is Boltzmann's constant (8.617×10⁻⁵ eV/K).
$T$ is absolute temperature (K).

Higher temperatures exponentially increase $D$, allowing atoms to overcome barriers more easily.

Fick's Second Law: Time Evolution

Combining Fick's First Law with mass conservation ($\partial C/\partial t = -\nabla \cdot J$) yields Fick's Second Law, which governs how concentrations change over time:

\frac{\partial C}{\partial t} = \nabla \cdot \left( D \, \nabla C \right) = D \, \nabla^2 C

(assuming $D$ is constant in space). This is the classic diffusion equation: concentrations smooth out over time, with sharper gradients diffusing faster.

Multi-Component Systems: The Onsager Framework

When multiple species diffuse simultaneously (e.g., dopants Ni²⁺, Co²⁺, Fe²⁺ in MgO), they interact. In vacancy-mediated diffusion, all species compete for the same vacancies. This creates cross-coupling: a gradient in species $j$ can drive a flux in species $i$, even if $i$ has no gradient.

The Onsager formalism generalizes Fick's law to account for these interactions. Instead of a single coefficient $D$, we have a diffusion matrix $D_{ij}$:

J_i = -\sum_j D_{ij}(C,T) \, \nabla C_j

Diagonal terms ($D_{ii}$): Self-diffusion of species $i$ driven by its own gradient.
Off-diagonal terms ($D_{ij}$, $i \neq j$): Cross-diffusion, where species $i$ responds to gradients in species $j$.

What is the Onsager Approximation?

Onsager's reciprocal relations state that the diffusion matrix must satisfy certain symmetry constraints based on microscopic reversibility. In the context of vacancy-mediated diffusion, the Onsager framework provides a rigorous way to calculate $D_{ij}$ from the individual tracer diffusivities and local concentrations.

The key insight: all species share the same vacancy mechanism, so their motions are fundamentally coupled. The Manning–Darken model (Section 2) gives explicit expressions for $D_{ij}$ consistent with Onsager's principles.

Time Evolution in Multi-Component Systems

Just as in the single-component case, we combine the generalized flux law with mass conservation to obtain:

\frac{\partial C_i}{\partial t} = \nabla \cdot \left( \sum_{j} D_{ij} \, \nabla C_j \right)

This is the governing equation solved by this tool. Unlike Fick's Second Law, the matrix $D_{ij}$ depends on the local concentrations of all species, so it must be recalculated at every spatial point and every timestep. This dynamic coupling produces surprising effects:

Uphill diffusion: Species can move against their own gradient.
Kirkendall effect: Asymmetric diffusion shifts the original interface.
Vacancy wind: Fast-diffusing species "drag" slower species along.

1.1 Tracer Diffusivities $D_i^*$

Each species (including the matrix cations, e.g., Mg²⁺ in MgO) has an intrinsic tracer diffusivity $D_i^*(T)$, which is the diffusion coefficient measured in a dilute limit when the species does not interact with others:

D_i^*(T) = D_{0,i} \, \exp\left( -\frac{Q_i}{k_B T} \right)

These $D_i^*$ values are the inputs to the Onsager diffusion matrix. They capture how easily each species jumps through the lattice in isolation. The matrix $D_{ij}$ is constructed from these tracer diffusivities plus concentration-dependent corrections.

2 The Manning–Darken Model: Building the Diffusion Matrix

In crystalline materials like MgO, cation diffusion occurs primarily through the vacancy mechanism: an atom jumps into a neighboring vacant lattice site. All cation species (Mg²⁺, Ni²⁺, Co²⁺, etc.) share the same pool of cation vacancies. This shared mechanism creates the cross-coupling described by the Onsager matrix.

The Manning–Darken model provides explicit formulas for the diffusion matrix elements $D_{ij}$ in terms of the tracer diffusivities $D_i^*$ and local concentrations $C_i$. These expressions respect Onsager's reciprocal relations and capture the physics of vacancy-mediated transport.

2.1 Diagonal Elements (Self-Diffusion)

The diagonal elements $D_{ii}$ describe how species $i$ responds to its own concentration gradient:

D_{ii} = D_i^* + C_i \left( D_{\text{matrix}}^* - D_i^* \right)

Physical interpretation: At low concentration ($C_i \to 0$), $D_{ii} \approx D_i^*$ recovers the tracer limit. As $C_i$ increases, the effective diffusivity blends toward the matrix diffusivity $D_{\text{matrix}}^*$ (the tracer diffusivity of the host cations, e.g., Mg²⁺ in MgO). This reflects the fact that in a concentrated alloy, the environment becomes more like the matrix.

2.2 Off-Diagonal Elements (Cross-Diffusion)

The off-diagonal elements $D_{ij}$ ($i \neq j$) describe how species $i$ responds to gradients in species $j$:

D_{ij} = C_i \left( D_{\text{matrix}}^* - D_j^* \right)

Physical interpretation: This term arises from the vacancy wind effect. When species $j$ diffuses (creating a flux of $j$-atoms), it generates a counter-flow of vacancies. Species $i$ can "surf" on this vacancy current, causing a flux of $i$-atoms even when $i$ has no concentration gradient.

If $D_j^* \ll D_{\text{matrix}}^*$ (slow dopant), then $D_{ij} > 0$: a gradient in $j$ pulls $i$ along.
If $D_j^* \gg D_{\text{matrix}}^*$ (fast dopant), then $D_{ij} < 0$: species $i$ moves against the $j$-gradient.

2.3 Why the Matrix Diffusivity Appears

The quantity $D_{\text{matrix}}^*$ is the tracer diffusivity of the matrix cations (e.g., Mg²⁺ in MgO). It sets the reference scale for vacancy motion. In a pure MgO crystal, vacancies would move at a rate determined by Mg²⁺ jumps. When dopants are added, their different jump rates create local imbalances in the vacancy concentration, which couples all species together.

The Manning–Darken expressions naturally interpolate between two limits:
(1) Dilute limit ($C_i \to 0$): Each dopant sees only the matrix, so $D_{ii} \to D_i^*$ and cross-terms vanish.
(2) Concentrated limit: The environment becomes dominated by the alloy composition, and all diffusivities blend toward a common effective value.

2.4 Final Flux Law with Thermodynamic Corrections

The Manning–Darken expressions give the kinetic diffusion matrix. In non-ideal solutions, we must also account for thermodynamic driving forces (chemical potential gradients). The solver includes these effects through the thermodynamic factor $\Phi_i$:

J_i(x) = \Phi_i(x) \left[ -\sum_j D_{ij}(x) \, \nabla C_j(x) \right]

For a regular solution (Section 5), $\Phi_i$ depends on the interaction parameter $\Omega$:

\Phi_i = 1 - \frac{2 \Omega \, C_i(1 - C_i)}{k_B T}

When $\Phi_i < 0$, the system enters the spinodal region and phase separates (uphill diffusion occurs naturally). When $\Phi_i > 1$, thermodynamic repulsion accelerates diffusion. This thermodynamic correction multiplies the kinetic flux, capturing effects like interface sharpening and demixing.

3 Numerical Algorithm & Stability

The solver uses an explicit finite-difference scheme. Because the system is strongly coupled, numerical stability depends on the largest eigenvalue of the diffusion matrix.

\Delta t \le \frac{\Delta x^2}{2\,\lambda_{\max}}

Instead of computing eigenvalues explicitly, the code uses Gershgorin bounds to estimate $\lambda_{\max}$ locally at every grid point. This ensures stability even for systems with strong coupling or highly unequal diffusivities.

4 Viewing the Results

Concentration (2D)

Standard diffusion curves.
• Asymmetry: Kirkendall shift.
• Local bumps: Vacancy drag.
• Flat zones: Flux cancellation.

Flux Surface (3D)

A time–position surface plot of $J(x,t)$.
• Reversal: Flux changing sign.
• Ridges: Tilted flow paths.

Time Traces (3D)

Evolution of concentration markers.
• Backwards movement: Uphill diffusion.
• Crossing: Species overtaking.

5 Non-Ideal Thermodynamics ($\Omega$)

Deviations from ideal mixing are included through the regular-solution parameter $\Omega$. It modifies the chemical potential curvature:

\Phi = 1 - \frac{2\Omega\, C(1-C)}{k_B T}

Spinodal Decomposition

If $\Phi < 0$, the system is in the spinodal region. Diffusion drives demixing (phase separation) instead of smoothing.

• $\Omega = 0$: Ideal solution.
• $\Omega > 0$: Repulsive (miscibility gaps possible).
• Large $\Omega$: Sharp interfaces.

6 Reference Values ($\Omega$)

Experimental thermodynamic assessments ($1 \text{ kJ/mol} \approx 0.0104 \text{ eV}$).

System (Matrix-Dopant)	$\Omega$ (eV)	Type	Physical Reason
MgO Matrix
MgO - NiO	0.00	Ideal	Complete solid solution (Similar radii).
MgO - CoO	0.05	Near Ideal	Very slight repulsion.
MgO - FeO	0.04	Near Ideal	Mixes well at high Temp.
MgO - ZnO	0.18	Repulsive	Structural mismatch (Rock Salt vs Wurtzite).
MgO - CuO	0.21	Repulsive	Jahn-Teller distortion.
MgO - CaO	0.55	Miscibility Gap	Ca²⁺ is large. Phase separates < 2000K.
MgO - Cr₂O₃	> 0.80	Precipitation	Solubility limit simulation.
Metallic / Test
Cu - Zn (Brass)	-0.10	Ordering	Prefer A-B neighbors (Intermetallics).
Au - Cu	-0.08	Ordering	Forms ordered compounds.
Simulation Test	-0.15	Strong Ordering	Use to see "box-like" profiles.

7 Advanced Topics & Limits

This solver implements vacancy-mediated Onsager diffusion using a simplified physical model. Limitations include:

1D geometry: The system is a one-dimensional slab.
Constant prefactors: $D_0$ and $Q$ are concentration-independent.
No defect chemistry: Charged vacancies or Frenkel pairs are not explicitly tracked.
Single Crystal: No grain boundaries or dislocations.

Even with these simplifications, the solver reproduces Kirkendall shifts, uphill diffusion, vacancy wind, and interface sharpening.

8 Notes & References

Usage: Intended for teaching, research, and rapid prototyping.

Citation: If used in academic work, please cite "Onsager Multicomponent Diffusion Solver v8.4, 2025" and the Manning-Darken coupling expressions.

Tech Stack: WebWorker for physics engine, Chart.js (2D), Plotly (3D).