The copper–nickel ($\text{Cu–Ni}$) system is the canonical metallurgical model of a binary isomorphous system. Both elements crystallise in the face-centred cubic ($\text{fcc}$) structure and possess complete mutual solubility across the entire phase diagram in the solid state.
Consider a semi-infinite diffusion couple formed by welding a semi-infinite block of pure copper to a semi-infinite block of pure nickel along a planar contact surface at $t = 0$. Let the laboratory coordinate axis $x$ be normal to the welding interface, oriented such that:
When the couple is held at an elevated isothermal annealing temperature $T$ below the melting point of copper ($T_{\text{m, Cu}} = 1085\ ^\circ\text{C}$), atoms exchange sites with vacancies. Because copper has a significantly lower activation energy for self-diffusion than nickel, copper atoms jump into vacancies with a higher frequency than nickel atoms.
This kinetic asymmetry leads directly to:
To model diffusion in a crystal where lattice planes may move, two distinct coordinate frames must be defined:
In the lattice frame (denoted by the prime superscript $'$), the intrinsic fluxes of copper ($J_{\text{Cu}}'$) and nickel ($J_{\text{Ni}}'$) are proportional to their respective concentration gradients:
\begin{align} J_{\text{Cu}}' &= - D_{\text{Cu}} \frac{\partial C_{\text{Cu}}}{\partial x'} \label{eq:jcu_lat} \\ J_{\text{Ni}}' &= - D_{\text{Ni}} \frac{\partial C_{\text{Ni}}}{\partial x'} = + D_{\text{Ni}} \frac{\partial C_{\text{Cu}}}{\partial x'} \label{eq:jni_lat} \end{align}where $C_{\text{Cu}}$ and $C_{\text{Ni}}$ are molar concentrations per unit volume ($C_{\text{Cu}} + C_{\text{Ni}} = C_{\text{total}} = 1/\Omega_{\text{a}}$, with $\Omega_{\text{a}}$ being the average molar volume, assumed constant), and $D_{\text{Cu}}$ and $D_{\text{Ni}}$ are the intrinsic diffusion coefficients of copper and nickel in the alloy.
The intrinsic mobilities obey Arrhenius relationships:
\begin{equation} D_{\text{Cu}}(T) = D_{0, \text{Cu}} \exp\left(-\frac{Q_{\text{Cu}}}{R T}\right), \quad D_{\text{Ni}}(T) = D_{0, \text{Ni}} \exp\left(-\frac{Q_{\text{Ni}}}{R T}\right) \end{equation}Because copper has a lower melting point ($1085\ ^\circ\text{C} = 1358\text{ K}$) than nickel ($1455\ ^\circ\text{C} = 1728\text{ K}$), its cohesive energy and vacancy activation barrier are substantially lower:
\begin{equation} Q_{\text{Cu}} \approx 211\text{ kJ/mol} \quad < \quad Q_{\text{Ni}} \approx 280\text{ kJ/mol} \end{equation}At all metallurgical temperatures relevant to this diffusion couple (e.g. $700\text{--}1050\ ^\circ\text{C}$), $D_{\text{Cu}}$ exceeds $D_{\text{Ni}}$ by a factor of approximately $4$ to $6$.
Because $D_{\text{Cu}} > D_{\text{Ni}}$ and $\frac{\partial C_{\text{Cu}}}{\partial x} < 0$, the number of copper atoms leaving the copper side per second exceeds the number of nickel atoms arriving to take their place. To prevent spontaneous microscopic voids from forming instantly throughout the bulk lattice, the crystal planes must shift to preserve constant atomic site density.
Let $v$ be the drift velocity of the local crystal lattice planes relative to the laboratory frame. The total flux of species $i$ in the laboratory frame ($J_i$) is the sum of its intrinsic diffusive flux and the convective drift of the lattice:
\begin{align} J_{\text{Cu}} &= J_{\text{Cu}}' + C_{\text{Cu}} v = - D_{\text{Cu}} \frac{\partial C_{\text{Cu}}}{\partial x} + C_{\text{Cu}} v \label{eq:jcu_lab} \\ J_{\text{Ni}} &= J_{\text{Ni}}' + C_{\text{Ni}} v = - D_{\text{Ni}} \frac{\partial C_{\text{Ni}}}{\partial x} + C_{\text{Ni}} v \label{eq:jni_lab} \end{align}Under the assumption of constant total molar density ($C_{\text{Cu}} + C_{\text{Ni}} = C$), the net flux of matter across any plane fixed in the laboratory frame must balance ($J_{\text{Cu}} + J_{\text{Ni}} = 0$):
\begin{equation} \left( - D_{\text{Cu}} \frac{\partial C_{\text{Cu}}}{\partial x} + C_{\text{Cu}} v \right) + \left( - D_{\text{Ni}} \frac{\partial C_{\text{Ni}}}{\partial x} + C_{\text{Ni}} v \right) = 0 \end{equation}Recognising that $\frac{\partial C_{\text{Ni}}}{\partial x} = -\frac{\partial C_{\text{Cu}}}{\partial x}$ and $C_{\text{Cu}} + C_{\text{Ni}} = 1$ (expressed in mole fractions):
\begin{equation} - (D_{\text{Cu}} - D_{\text{Ni}}) \frac{\partial C_{\text{Cu}}}{\partial x} + v = 0 \end{equation}Rearranging yields Darken's first equation for the lattice velocity:
\begin{equation} \boxed{v = (D_{\text{Cu}} - D_{\text{Ni}}) \frac{\partial C_{\text{Cu}}}{\partial x}} \label{eq:darken_first} \end{equation}Substituting the velocity expression \eqref{eq:darken_first} back into the laboratory flux equation \eqref{eq:jcu_lab} for copper:
\begin{align} J_{\text{Cu}} &= - D_{\text{Cu}} \frac{\partial C_{\text{Cu}}}{\partial x} + C_{\text{Cu}} (D_{\text{Cu}} - D_{\text{Ni}}) \frac{\partial C_{\text{Cu}}}{\partial x} \nonumber \\ &= - \left[ (1 - C_{\text{Cu}}) D_{\text{Cu}} + C_{\text{Cu}} D_{\text{Ni}} \right] \frac{\partial C_{\text{Cu}}}{\partial x} \nonumber \\ &= - \left( C_{\text{Ni}} D_{\text{Cu}} + C_{\text{Cu}} D_{\text{Ni}} \right) \frac{\partial C_{\text{Cu}}}{\partial x} \end{align}Comparing this directly to Fick's first law in the laboratory frame, $J_{\text{Cu}} = - \tilde{D} \frac{\partial C_{\text{Cu}}}{\partial x}$, yields Darken's second equation:
\begin{equation} \boxed{\tilde{D} = C_{\text{Ni}} D_{\text{Cu}} + C_{\text{Cu}} D_{\text{Ni}}} \label{eq:darken_second} \end{equation}where $\tilde{D}$ is the interdiffusion coefficient. Fick's second law in the laboratory frame then takes the standard form:
\begin{equation} \frac{\partial C_{\text{Cu}}}{\partial t} = \frac{\partial}{\partial x}\left( \tilde{D} \frac{\partial C_{\text{Cu}}}{\partial x} \right) \label{eq:fick_second} \end{equation}In an experimental diffusion couple, the laboratory frame of reference is established by the Matano interface ($x = 0$). This plane represents the original contact position before diffusion took place, defined rigorously by conservation of mass:
\begin{equation} \int_{-\infty}^{x_{\text{M}}} \left[ C_{\text{Cu}}(-\infty) - C_{\text{Cu}}(x, t) \right] dx = \int_{x_{\text{M}}}^{+\infty} \left[ C_{\text{Cu}}(x, t) - C_{\text{Cu}}(+\infty) \right] dx \end{equation}For our boundary conditions ($C_{\text{Cu}}(-\infty) = 1$ and $C_{\text{Cu}}(+\infty) = 0$), setting $x_{\text{M}} = 0$ requires:
\begin{equation} \boxed{\int_{-\infty}^{0} (1 - C_{\text{Cu}}) \, dx = \int_{0}^{+\infty} C_{\text{Cu}} \, dx} \label{eq:matano_condition} \end{equation}This condition guarantees that the mass of copper that has diffused across $x = 0$ into the right-hand bar equals the mass of copper depleted from the left-hand bar.
When evaluating the system under the average interdiffusivity approximation $\tilde{D} \approx \text{constant}$ evaluated at the interface ($C = 0.50$), equation \eqref{eq:fick_second} reduces to the linear diffusion equation:
\begin{equation} \frac{\partial C_{\text{Cu}}}{\partial t} = \tilde{D} \frac{\partial^2 C_{\text{Cu}}}{\partial x^2} \label{eq:linear_diff} \end{equation}subject to the step initial conditions:
\begin{equation} C_{\text{Cu}}(x, 0) = \begin{cases} 1, & x < 0 \\ 0, & x > 0 \end{cases} \end{equation}and infinite boundary conditions $C_{\text{Cu}}(-\infty, t) = 1$, $C_{\text{Cu}}(+\infty, t) = 0$ for all $t > 0$.
Introducing the dimensionless similarity variable $\lambda \equiv \frac{x}{2\sqrt{\tilde{D} t}}$, partial derivatives transform as:
\begin{equation} \frac{\partial C_{\text{Cu}}}{\partial t} = -\frac{\lambda}{2 t} \frac{d C_{\text{Cu}}}{d \lambda}, \quad \frac{\partial^2 C_{\text{Cu}}}{\partial x^2} = \frac{1}{4 \tilde{D} t} \frac{d^2 C_{\text{Cu}}}{d \lambda^2} \end{equation}Substituting these into \eqref{eq:linear_diff} converts the partial differential equation into an ordinary differential equation:
\begin{equation} \frac{d^2 C_{\text{Cu}}}{d \lambda^2} + 2 \lambda \frac{d C_{\text{Cu}}}{d \lambda} = 0 \end{equation}Integrating once:
\begin{equation} \frac{d C_{\text{Cu}}}{d \lambda} = A \exp(-\lambda^2) \end{equation}Integrating a second time from $\lambda = 0$ to $\lambda$:
\begin{equation} C_{\text{Cu}}(\lambda) = C_{\text{Cu}}(0) + A \int_0^{\lambda} \exp(-\xi^2) \, d\xi = C_{\text{Cu}}(0) + A \frac{\sqrt{\pi}}{2} \operatorname{erf}(\lambda) \end{equation}Applying the boundary conditions:
Adding both equations immediately yields:
\begin{equation} 2 C_{\text{Cu}}(0) = 1 \implies \boxed{C_{\text{Cu}}(0) = 0.50} \end{equation}Subtracting gives $A = -\frac{1}{\sqrt{\pi}}$. The analytical concentration profiles for copper and nickel are therefore:
\begin{align} \boxed{C_{\text{Cu}}(x, t) = \frac{1}{2} \operatorname{erfc}\left( \frac{x}{2\sqrt{\tilde{D} t}} \right) = \frac{1}{2}\left[ 1 - \operatorname{erf}\left( \frac{x}{2\sqrt{\tilde{D} t}} \right) \right]} \label{eq:erfc_cu} \\ \boxed{C_{\text{Ni}}(x, t) = 1 - C_{\text{Cu}}(x, t) = \frac{1}{2}\left[ 1 + \operatorname{erf}\left( \frac{x}{2\sqrt{\tilde{D} t}} \right) \right]} \label{eq:erf_ni} \end{align}In 1947, Ernest Kirkendall placed inert molybdenum wire markers at the interface of a copper–brass diffusion couple and observed that the markers moved relative to the outer ends of the bar during annealing.
Because the markers are inert, they do not dissolve in the lattice; they remain anchored to the local crystal lattice planes. Consequently, the velocity of the markers ($v_{\text{K}}$) in the laboratory frame is identical to the local lattice drift velocity at the marker plane ($x = x_{\text{K}}$):
\begin{equation} v_{\text{K}} = \left. (D_{\text{Cu}} - D_{\text{Ni}}) \frac{\partial C_{\text{Cu}}}{\partial x} \right|_{x = x_{\text{K}}} \label{eq:marker_vel} \end{equation}Differentiating equation \eqref{eq:erfc_cu} with respect to $x$:
\begin{equation} \frac{\partial C_{\text{Cu}}}{\partial x} = -\frac{1}{2\sqrt{\pi \tilde{D} t}} \exp\left( -\frac{x^2}{4 \tilde{D} t} \right) \label{eq:dc_dx} \end{equation}At the interface ($x \approx 0$):
\begin{equation} \left. \frac{\partial C_{\text{Cu}}}{\partial x} \right|_{x=0} = -\frac{1}{2\sqrt{\pi \tilde{D} t}} \end{equation}Substituting this into \eqref{eq:marker_vel}:
\begin{equation} v_{\text{K}}(t) = \frac{dx_{\text{K}}}{dt} = - \frac{D_{\text{Cu}} - D_{\text{Ni}}}{2\sqrt{\pi \tilde{D} t}} \label{eq:dxk_dt} \end{equation}Integrating equation \eqref{eq:dxk_dt} from $t = 0$ (where $x_{\text{K}}(0) = 0$) to time $t$:
\begin{equation} x_{\text{K}}(t) = \int_0^t v_{\text{K}}(t') \, dt' = - \frac{D_{\text{Cu}} - D_{\text{Ni}}}{2\sqrt{\pi \tilde{D}}} \int_0^t (t')^{-1/2} \, dt' \end{equation} \begin{equation} \boxed{x_{\text{K}}(t) = - \frac{D_{\text{Cu}} - D_{\text{Ni}}}{\sqrt{\pi \tilde{D}}} \sqrt{t} = 2 \beta \sqrt{t}} \label{eq:marker_shift} \end{equation}where $\beta = -\frac{D_{\text{Cu}} - D_{\text{Ni}}}{2\sqrt{\pi \tilde{D}}} < 0$ is a constant of parabolic growth.
In classical Darken analysis, it is assumed that vacancies are maintained in perfect thermal equilibrium everywhere ($C_{\text{v}} = C_{\text{v}}^{\text{eq}}$) by rapid creation and destruction at edge dislocation jogs. However, in real crystals with finite dislocation density, this equilibrium cannot be sustained instantaneously during interdiffusion.
In an interstitial-free substitutional lattice, every atomic jump is accompanied by an equal and opposite vacancy jump. Therefore, the intrinsic flux of vacancies in the lattice frame ($J_{\text{v}}'$) must balance the total intrinsic atomic flux:
\begin{equation} J_{\text{v}}' + J_{\text{Cu}}' + J_{\text{Ni}}' = 0 \implies J_{\text{v}}' = - (J_{\text{Cu}}' + J_{\text{Ni}}') \end{equation}Substituting equations \eqref{eq:jcu_lat} and \eqref{eq:jni_lat}:
\begin{equation} \boxed{J_{\text{v}}' = (D_{\text{Cu}} - D_{\text{Ni}}) \frac{\partial C_{\text{Cu}}}{\partial x}} \label{eq:j_vac} \end{equation}Because $\frac{\partial C_{\text{Cu}}}{\partial x} < 0$ and $D_{\text{Cu}} > D_{\text{Ni}}$, $J_{\text{v}}' < 0$. Vacancies stream continuously from the nickel side into the copper side.
The time rate of change of the vacancy concentration $C_{\text{v}}(x, t)$ is governed by the divergence of the vacancy flux and the rate of vacancy creation/annihilation at dislocation climb jogs:
\begin{equation} \frac{\partial C_{\text{v}}}{\partial t} = - \frac{\partial J_{\text{v}}'}{\partial x} - \dot{\rho}_{\text{sink}} \end{equation}Following the standard linear kinetic model, the sink/source term is proportional to the deviation of the vacancy concentration from its local thermal equilibrium value $C_{\text{v}}^{\text{eq}}(T)$:
\begin{equation} \dot{\rho}_{\text{sink}} = \frac{C_{\text{v}} - C_{\text{v}}^{\text{eq}}}{\tau_{\text{v}}} \end{equation}where $\tau_{\text{v}}$ is the vacancy relaxation time, related to the mean dislocation density $\rho_{\text{d}}$ by:
\begin{equation} \tau_{\text{v}} \approx \frac{1}{D_{\text{v}} \rho_{\text{d}}} \end{equation}($D_{\text{v}}$ being the vacancy diffusivity). Under quasi-steady-state conditions ($\frac{\partial C_{\text{v}}}{\partial t} \approx 0$ relative to rapid vacancy diffusion):
\begin{equation} - \frac{\partial J_{\text{v}}'}{\partial x} - \frac{C_{\text{v}} - C_{\text{v}}^{\text{eq}}}{\tau_{\text{v}}} = 0 \end{equation}Rearranging gives the excess vacancy concentration $\Delta C_{\text{v}} = C_{\text{v}} - C_{\text{v}}^{\text{eq}}$:
\begin{equation} \Delta C_{\text{v}}(x, t) = - \tau_{\text{v}} \frac{\partial J_{\text{v}}'}{\partial x} = - \tau_{\text{v}} (D_{\text{Cu}} - D_{\text{Ni}}) \frac{\partial^2 C_{\text{Cu}}}{\partial x^2} \label{eq:delta_cv} \end{equation}Differentiating the concentration gradient \eqref{eq:dc_dx} once more with respect to $x$ yields the second spatial derivative:
\begin{equation} \frac{\partial^2 C_{\text{Cu}}}{\partial x^2} = \frac{\partial}{\partial x}\left[ -\frac{1}{2\sqrt{\pi \tilde{D} t}} \exp\left( -\frac{x^2}{4 \tilde{D} t} \right) \right] = \frac{x}{4 \sqrt{\pi} (\tilde{D} t)^{3/2}} \exp\left( -\frac{x^2}{4 \tilde{D} t} \right) \label{eq:d2c_dx2} \end{equation}Substituting \eqref{eq:d2c_dx2} into the excess vacancy formula \eqref{eq:delta_cv}:
\begin{equation} \boxed{\frac{\Delta C_{\text{v}}(x, t)}{C_{\text{v}}^{\text{eq}}} = - \frac{\tau_{\text{v}} (D_{\text{Cu}} - D_{\text{Ni}})}{C_{\text{v}}^{\text{eq}}} \cdot \frac{x}{4 \sqrt{\pi} (\tilde{D} t)^{3/2}} \exp\left( -\frac{x^2}{4 \tilde{D} t} \right)} \label{eq:vac_ratio_final} \end{equation}Letting $z \equiv \frac{x}{2\sqrt{\tilde{D} t}}$, the spatial shape is dictated by the dimensionless function $-z \exp(-z^2)$:
Differentiating $-z \exp(-z^2)$ with respect to $z$ and equating to zero:
\begin{equation} \frac{d}{dz}\left( -z e^{-z^2} \right) = -(1 - 2z^2) e^{-z^2} = 0 \implies z = \pm \frac{1}{\sqrt{2}} \end{equation}The physical locations of peak vacancy supersaturation ($x_{\text{max}}$) and peak vacancy subsaturation ($x_{\text{min}}$) are:
\begin{equation} \boxed{x_{\text{max}} = - \sqrt{2 \tilde{D} t}, \quad x_{\text{min}} = + \sqrt{2 \tilde{D} t}} \label{eq:peak_positions} \end{equation}The numerical values used in the simulation and their physical provenance are summarised in the table below:
| Symbol | Physical quantity | Numerical value / expression | Source / remarks |
|---|---|---|---|
| $D_{0, \text{Cu}}$ | Pre-exponential factor for $\text{Cu}$ | $7.8 \times 10^{-5}\text{ m}^2/\text{s}$ | Tracer self-diffusion in $\text{Cu–Ni}$ alloys |
| $Q_{\text{Cu}}$ | Activation enthalpy for $\text{Cu}$ | $211.0\text{ kJ/mol}$ ($2.19\text{ eV}$) | Lower activation barrier (lower melting point) |
| $D_{0, \text{Ni}}$ | Pre-exponential factor for $\text{Ni}$ | $1.9 \times 10^{-4}\text{ m}^2/\text{s}$ | Tracer self-diffusion in $\text{Cu–Ni}$ alloys |
| $Q_{\text{Ni}}$ | Activation enthalpy for $\text{Ni}$ | $280.0\text{ kJ/mol}$ ($2.90\text{ eV}$) | Higher activation barrier (higher melting point) |
| $R$ | Universal gas constant | $8.31446\text{ J/(mol K)}$ | Standard physical constant |
| $T_{\text{m, Cu}}$ | Melting point of pure copper | $1085\ ^\circ\text{C}$ ($1358.15\text{ K}$) | Defines upper thermal stability limit of couple |
| $T_{\text{m, Ni}}$ | Melting point of pure nickel | $1455\ ^\circ\text{C}$ ($1728.15\text{ K}$) | Demonstrates higher cohesive energy of nickel |
| $\tilde{D}$ | Darken interdiffusivity | $0.5(D_{\text{Cu}} + D_{\text{Ni}})$ | Evaluated at equiatomic interface ($C = 0.50$) |
| $x_{\text{K}}(t)$ | Kirkendall marker shift | $- \frac{D_{\text{Cu}} - D_{\text{Ni}}}{\sqrt{\pi \tilde{D}}} \sqrt{t}$ | Displacement into copper bar ($x < 0$) |
| $x_{\text{peak}}$ | Position of maximum voiding | $- \sqrt{2 \tilde{D} t}$ | Peak of vacancy supersaturation profile |
In summary, the full mathematical model implemented in the interactive simulator comprises:
\begin{align} C_{\text{Cu}}(x, t) &= \begin{cases} \Theta(-x), & t = 0 \\ \frac{1}{2}\operatorname{erfc}\left(\frac{x}{2\sqrt{\tilde{D} t}}\right), & t > 0 \end{cases} \\[6pt] C_{\text{Ni}}(x, t) &= 1 - C_{\text{Cu}}(x, t) \\[6pt] v_{\text{K}}(t) &= (D_{\text{Cu}} - D_{\text{Ni}}) \left.\frac{\partial C_{\text{Cu}}}{\partial x}\right|_{x = x_{\text{K}}} \\[6pt] x_{\text{K}}(t) &= - \frac{D_{\text{Cu}} - D_{\text{Ni}}}{\sqrt{\pi \tilde{D}}} \sqrt{t} \\[6pt] \frac{C_{\text{v}}(x, t)}{C_{\text{v}}^{\text{eq}}} &= 1 - \tau_{\text{v}} (D_{\text{Cu}} - D_{\text{Ni}}) \frac{\partial^2 C_{\text{Cu}}}{\partial x^2} \end{align}