Onsager Diffusion v9.3

1. Overview

This tool simulates multi-component cation diffusion in oxide crystals using the Onsager transport formalism. Unlike simple Fickian diffusion where species diffuse independently, real crystals exhibit correlated motion. Every jump requires a vacancy, and vacancies are shared by all species, producing cross-coupling between fluxes.

2. Physical Models

Initially the deposited film contains 100% selected non-Mg cations, while the substrate is pure MgO. The domain spans film plus slab. Both external boundaries are sealed; Mg can diffuse into the film during annealing. Set time to zero to inspect the initial state. One finite-volume cell crossing the initial interface contains the volume-weighted average of both layers. Film percentages are cation-site fractions, not oxide mass or formula-unit fractions; this constant-site-density model does not resolve separate oxide phases, oxygen transport or density changes.

Oxygen, film composition & density assumptions

What the percentages mean. Film inputs and concentration profiles describe the fractions of cation sites occupied by each element. A film containing four cations at 25% each has 100% of its cation sites assigned to those elements. Oxygen is outside that percentage total; these values are not oxide mass percentages or percentages of all atoms.

Oxygen is implicit. The model does not solve oxygen concentration, oxygen diffusion, exchange with the atmosphere, or the oxygen-to-cation ratio of each oxide. Cation valence enters a phenomenological vacancy correction, but this is not an explicit oxygen or charge-balance calculation.

Film density correction. Each film cation carries its oxide (e.g. Cr₂O₃, TiO₂) with molar mass, bulk density and cations per formula unit. The film's cation density is $n_{film}=1/\sum_i x_i\,M_i/(\rho_i\nu_i)$ (ideal mixing of volume per cation, $\nu_i$ = cations per formula unit). The solver runs in MgO-equivalent depth: the film thickness is multiplied by $n_{film}/n_{MgO}$ ($n_{MgO}$ = 88,825 mol/m³), so every length unit holds the same moles of cations per area as MgO, while the single-crystal substrate keeps its real thickness. Profiles, the interface marker and the CSV use this depth axis, and fluxes in mol/m²/s remain correct. Porosity, thermal expansion and volume change on mixing are not modelled, and the film is still treated as a full cation lattice.

Flux conversion. To display and export fluxes in mol/m²/s, the code multiplies its fraction-based flux by 88,825 mol of cation sites/m³. For MgO, this corresponds to approximately 3.58 g/cm³. This factor converts the flux units; it is valid throughout because the solver depth axis is MgO-equivalent (see film density correction).

Comparison with an experimental film. Converting deposited oxide amounts into a cation inventory requires the oxide stoichiometries and an appropriate density or measured mass per unit area. Quantitatively modeling a film whose site density differs from MgO also requires extending the transport model. Film thickness and cation percentages alone do not supply those properties.

2.1 Onsager Flux Equation

For each species $i$, the flux is governed by the gradients of all concentrations, split into kinetic cross-coupling and thermodynamic driving forces:

J_i = - \sum_{j} \sum_{k} D_{ik}^{kin} \Phi_{kj} \nabla C_j

$D_{ik}^{kin}$: Onsager-Manning kinetic coefficients (Vacancy Wind / Site sharing).
$\Phi_{kj}$: Thermodynamic factor matrix derived from the $\Omega$-matrix.

2.2 Kinetics: The Q-Matrix & Vacancies

The local tracer diffusivity depends on the concentration-weighted activation barrier ($Q$) and the dynamic vacancy fraction ($V_{factor}$):

D_i^*(x) = D_{0,i}\, \exp\left(-\frac{\sum_j C_j(x)\, Q_{ij}}{k_B T}\right) \times V_{factor}(x)

The same $V_{factor}$ multiplies every cation, including the host Mg, because all of them jump through the same cation vacancies. $V_{factor}$ is the local free cation-vacancy fraction divided by its value in the undoped slab, so the tabulated $D_0$ and $Q$ (measured on nominally pure crystals) apply unchanged in the slab.

[V_{Mg}]-[V_O]=d,\qquad [V_{Mg}][V_O]=K_S=e^{-H_S/k_BT},\qquad d=V_{bg}+\sum_i C_i\frac{Z_i-2}{2} $$ $$ V_{factor}=\frac{[V_{Mg}]\left(1-f_{assoc}\tanh(\max(0,d-V_{bg})/V_{sat})\right)}{[V_{Mg}]_{d=V_{bg}}}

$V_{bg}=10^{-4}$ is the extrinsic background of cation vacancies compensating trivalent impurities in nominally pure MgO (commercial crystals: roughly 15–1000 ppm). It is not an intrinsic limit.
$H_S = 6.7$ eV (≈650 kJ/mol Schottky-pair enthalpy) gives a temperature-dependent intrinsic level $\sqrt{K_S}$, ≈3×10⁻¹² at 1473 K. It matters only where acceptors outweigh the background.
Donors ($Z>2$): one cation vacancy per two M³⁺ or per M⁴⁺. A pure trivalent film therefore reaches $V_{factor}\approx 2500$ (after association), a pure Ti film ≈5000.
Acceptors ($Z<2$) are assumed compensated by oxygen vacancies only; hole compensation is not modelled, so cation-vacancy suppression (and the slowdown) in acceptor-rich regions is an upper bound.
$f_{assoc}=0.5$, $V_{sat}=0.02$: phenomenological vacancy–donor association, applied equally to all cations. Real M³⁺–vacancy pairs speed up the trivalent cation and slow the others (Van Orman et al. 2009), which this model does not separate.

Default Q matrix. $Q_{ij}=Q_i$ for every environment $j$, applied the same way to Mg and to the dopants. This reproduces the measured dilute tracer value in MgO and assumes no composition dependence of the barrier. Composition dependence then comes only from $V_{factor}$ and $\Omega$. Edit the Q matrix in the Interaction dialog when environment-dependent data exist.

2.3 Thermodynamics: The Ω-Matrix

All symmetric pair interactions, including Mg–dopant pairs, enter the regular-solution free energy. Index 0 denotes Mg; all cation fractions sum to one. The chemical-potential gradients drive the flux:

g = k_B T\sum_i C_i\ln C_i + h, \qquad h=\sum_{i<j}\Omega_{ij}C_iC_j $$ $$ F_i = \nabla C_i + \frac{C_i}{k_BT}\left(\sum_j\Omega_{ij}\nabla C_j-\nabla h\right) $$ $$ J_i = -D_i^*F_i + C_i\sum_k D_k^*F_k

For one dopant, the thermodynamic factor reduces to $1-2\Omega C(1-C)/(k_BT)$. Unstable or marginal compositions stop with a message; phase separation requires a gradient-energy model. The vacancy correction and most of the parameter database remain phenomenological.

Stability is checked with a Cholesky factorization of the free-energy Hessian on the composition simplex. The cation eliminated from the reduced coordinates is the most abundant one at each face (fraction ≥ 1/(n+1)), so the $1/C_{ref}$ term stays small even inside the Mg-free film and no precision is lost. Positive definiteness does not depend on this choice.

Fluxes are stored on N+1 faces: the sealed boundaries at $x=0$ and $x=L$ (always zero) and the N−1 interior faces at $(k-\tfrac12)\Delta x$. Both end cells have half width.

The finite-volume solver uses upwind drift for every cation, including Mg, with a common lattice velocity that enforces zero total cation flux. The time step follows current outgoing transport rates; failed trial steps are retried without clipping or renormalizing concentrations. There is no step-count cutoff. The main Run calculation button becomes Stop while the worker runs and immediately cancels it when clicked. Upwind transport adds numerical smoothing, so compare grid resolutions for sharp profiles. Reference: NIST FiPy manual, numerical schemes.

Reference: MIT lecture on spinodal decomposition and gradient energy.

3. Reference Values (MgO Host Matrix)

The defaults below are the exact values used in the simulation code. These are supplied model parameters; diffusivity also depends on temperature, composition and the vacancy factor.

3.1 Complete Cation Database (D₀ & Q)

Full list of cations available in the "Dopants" dropdown.

Ion	D₀ (m²/s)	Q (eV)	Source

3.2 Interaction Energies ($\Omega$)

All values are model estimates, generated from the code. Pairs not listed use the Default (0.05 eV). Negative Ω is attractive (favors mixing).

Defined Pair	$\Omega$ (eV)	Basis

4. Bibliography

Blank, S. L., & Pask, J. A. (1969). "Diffusion of Iron and Nickel in Magnesium Oxide Single Crystals." J. Am. Ceram. Soc. 52. doi
Crispin, K. L., & Van Orman, J. A. (2010). "Influence of the crystal field effect on chemical transport in Earth's mantle: Cr³⁺ and Ga³⁺ diffusion in periclase." Phys. Earth Planet. Inter. record
Doman, R. C., Barr, J. B., McNally, R. N., & Alper, A. M. (1963). "Phase Equilibria in the System CaO–MgO." J. Am. Ceram. Soc. 46.
Gourdin, W. H., & Kingery, W. D. (1979). "The defect structure of MgO containing trivalent cation solutes: shell model calculations." J. Mater. Sci. 14.
Harding, B. C. (1972). "Diffusion of Ni and Ba in MgO up to 2500 °C." Phys. Status Solidi B 50, 135.
Harding, B. C., & Price, D. M. (1972). "Cation self-diffusion in MgO up to 2350 °C." Phil. Mag. 26, 253–260. doi (²⁸Mg radiotracer, 1265–2350 °C)
Osenbach, J. W., & Stubican, V. S. (1983). "Grain-boundary diffusion of ⁵¹Cr in MgO and Cr-doped MgO." J. Am. Ceram. Soc. 66, 191–195. doi
Rungis, J., & Mortlock, A. J. (1966). "The diffusion of calcium in magnesium oxide." Phil. Mag. 14, 821–827. doi
Van Orman, J. A., & Crispin, K. L. (2010). "Diffusion in Oxides." Rev. Mineral. Geochem. 72, 757–825. doi
Van Orman, J. A., Li, C., & Crispin, K. L. (2009). "Aluminum diffusion and Al-vacancy association in periclase." Phys. Earth Planet. Inter. 172, 34–42.
Weber, G. W., Bitler, W. R., & Stubican, V. S. (1977). "Diffusion of ⁵¹Cr in MgO Crystals." J. Am. Ceram. Soc. 60, 61–64.
Wuensch, B. J., Steele, W. C., & Vasilos, T. (1973). "Cation self-diffusion in single-crystal MgO." J. Chem. Phys. 58, 5258–5266. (Stable-isotope mass spectrometry, 1000–2400 °C, argon.)
Wuensch, B. J., & Vasilos, T. (1962). "Diffusion of Transition Metal Ions in Single-Crystal MgO." J. Chem. Phys. 36, 2917–2922. doi

Onsager Diffusion Multicomponent diffusion · MgO

Concentration profile

Theory & Documentation