L3 · Chapter 6

Solute trapping in alloy phase-field models and the anti-trapping current

Sandbox Physics · 2026-10-30

When a dilute alloy freezes, the solid takes in only a fraction $k$ of the solute in the liquid next to it and pushes the rest ahead, where it has to diffuse away. That bottleneck makes the growing front unstable and turns it into dendrites. A phase-field model can follow such shapes only if its interface is thousands of times wider than a real one, and a wide interface traps solute in the solid that a real interface would reject. This chapter states the problem the model must solve, shows the extra current Karma added to fix the trapping, and tests it on a flat front where the right answer is known exactly.

Solute has to get out of the way

Hold a dilute binary alloy at a fixed temperature below its liquidus. The liquid touching the solid has concentration $c_l$ and the solid there takes $c_s=kc_l$ , with a partition coefficient $k<1$ . Each step forward of the front rejects the difference $c_l(1-k)$ per unit volume into the liquid, and the gradient of concentration ahead must carry it away. This is the Stefan condition. Diffusion in the solid is neglected.

The liquid concentration at the front is not fixed. A curved front holds less, through the capillary length $d_0$ , and a moving one less again, through a kinetic coefficient $\beta$ . Far away the liquid has concentration $c_\infty$ , below the equilibrium value $c_l^0$ ; the supersaturation $\Omega$ measures by how much. Karma’s runs use $k=0.15$ and $\Omega=0.55$ , so $c_\infty=(1-0.55\times0.85)\,c_l^0=0.5325\,c_l^0$ .

\partial_tc=D\nabla^2c,\qquad c_l(1-k)V_n=-D\,\partial_nc\big|_{+},\qquad \frac{c_l}{c_l^0}=1-(1-k)d_0\kappa-(1-k)\beta V_n,\qquad c_s=kc_l

The sharp-interface problem for isothermal growth. The first three relations are Echebarria et al., Eqs. 17–19: diffusion in the liquid, the Stefan condition on the liquid side

+

and the Gibbs–Thomson relation with kinetics. The last one uses the partition coefficient

k

to set the composition of the solid.

\Omega=\frac{c_l^0-c_\infty}{c_l^0(1-k)}

The scaled supersaturation of the melt.

Source: B. Echebarria · A. Karma

Why a thick interface traps solute

A phase field $\phi$ runs from $+1$ in the solid to $-1$ in the liquid over a width $W$ . A real interface is about a nanometre thick, but a dendrite is microns to millimetres across, so a simulation must take $W$ far larger than the real thickness. Echebarria, Folch, Karma and Plapp list what goes wrong. When solute diffuses quickly in the liquid but hardly at all in the solid, a wide interface adds three spurious effects, each proportional to $W$ : diffusion along the interface, an extra term from the stretching of a moving curved interface, and a jump of chemical potential across it.

The jump is the one you can see. Solute crossing the interface region has no time to escape before the solid closes over it, so the solid keeps more than $kc_l$ . The size of the effect is set by the Péclet number $p=WV/D$ , the interface width over the diffusion length. Karma notes that in a model without a correction the trapped amount scales as $p\ln p$ , which shrinks only slowly as the interface is made thinner.

Source: B. Echebarria · A. Karma

The anti-trapping current

Karma’s model couples the phase field to a chemical-potential variable $u$ that is zero at equilibrium. The solute obeys a continuity equation with two currents. The first is ordinary diffusion, with a function $q(\phi)$ that switches diffusion off in the solid. The second, the anti-trapping current, points from solid to liquid, acts only where the interface is moving, since it is proportional to $\partial_t\phi$ , and is proportional to $W$ , the same order as the error it cancels.

With $a_t=1/(2\sqrt2)$ the chemical-potential jump vanishes; with $q(\phi)$ from Eq. (20), so does diffusion along the interface. The stretching term vanishes too, apart from small corrections that Karma reports as small in computations. The thin-interface analysis then maps the model onto the sharp-interface problem with $d_0=a_1W/\lambda$ and a kinetic coefficient that vanishes when $\lambda=D\tau/(a_2W^2)$ . The Lab uses these choices, and also runs Karma’s comparison model: no anti-trapping current and $q=(1-\phi)/2$ .

\tau\,\partial_t\phi=W^2\nabla^2\phi+\phi-\phi^3-\frac{\lambda}{1-k}\,(1-\phi^2)^2\,(e^u-1),\qquad e^u=\frac{2c/c_l^0}{1+k-(1-k)\phi}

The phase-field equation without anisotropy, with

\phi=+1

in the solid, and the definition of

u

(Karma, Eqs. 16 and 18, with the usual choice of functions).

\partial_tc=-\nabla\cdot\mathbf j,\qquad \mathbf j=-Dc\,q(\phi)\nabla u-a_tW c_l^0(1-k)\,e^u\,\partial_t\phi\,\frac{\nabla\phi}{|\nabla\phi|},\qquad a_t=\frac{1}{2\sqrt2}

Solute continuity with the diffusive current and the anti-trapping current (Karma, Eqs. 2 and 17).

q(\phi)=\frac{1-\phi}{1+k-(1-k)\phi},\qquad d_0=\frac{a_1W}{\lambda},\qquad \lambda=\frac{D\tau}{a_2W^2}\;\Rightarrow\;\beta=0,\qquad a_1=0.8839,\ a_2=0.6267

The interpolation function that removes surface diffusion, and the relations that fix the capillary length and remove interface kinetics.

Source: A. Karma · B. Echebarria

Worked example: does a flat front leave the right solute behind?

Problem. Take $W=\tau=1$ , $k=0.15$ , $\Omega=0.55$ and $d_0/W=0.277$ , the coarser of Karma’s two ratios. Find $\lambda$ , $D$ and the unit of time $d_0^2/D$ , then predict what a flat front leaves in the solid.

Parameters. $\lambda=a_1/(d_0/W)=0.8839/0.277=3.191$ and $D=a_2\lambda=0.6267\times3.191=2.000$ . The unit of time is $0.277^2/2.000=0.0384\,\tau$ .

Prediction. A flat front has no curvature, and the kinetic coefficient is zero, so local equilibrium puts $c_l=c_l^0$ at the front and the solid must hold $kc_l^0=0.150\,c_l^0$ however fast the front moves.

Check. The route benchmark grows a flat front into the melt on a grid of $\Delta x=0.4\,W$ . As it slows, the Péclet number falls from $0.066$ to $0.0087$ . With the anti-trapping current the solid holds $0.1500$ to $0.1506\,c_l^0$ throughout. Without it the solid holds $0.267$ at $p=0.11$ , falling to $0.166$ at $p=0.0098$ . For $p\le0.043$ the excess over $k$ , divided by $kp|\ln p|$ , stays between $2.32$ and $2.40$ , consistent with Karma’s $p\ln p$ scaling.

Answer. $\lambda=3.19$ , $D=2.00$ , a time unit of $0.0384\,\tau$ , and $0.150\,c_l^0$ in the solid, which the anti-trapping model delivers to within $0.4\%$ and the standard model misses by $11\%$ to $78\%$ .

Source: A. Karma

Try it in the Lab

Compare the solute left by the two interface models

Choose A four-armed dendrite. Predict the solute readout once the arms are long: a little below $k=0.15$ because the tip is curved.
Choose Solute gets trapped, which switches off the anti-trapping current. Predict whether the solid holds more or less solute and whether the tip grows faster or slower, then watch.
Explain both changes from trapping: what happens to the solute that the front no longer rejects, and why that lets the tip run ahead.

Exercises with explained answers

What is the far-field concentration of the melt for $k=0.15$ and $\Omega=0.45$ ?
Show the answer
$c_\infty=[1-\Omega(1-k)]\,c_l^0=(1-0.45\times0.85)\,c_l^0=0.6175\,c_l^0$ . The melt is closer to equilibrium, so the crystal grows more slowly.
A dendrite tip grows at $Vd_0/D=0.018$ with $d_0/W=0.277$ . What is its Péclet number $WV/D$ , and what does the flat-front test say about trapping in the standard model at that value?
Show the answer
$WV/D=(W/d_0)(Vd_0/D)=0.018/0.277=0.065$ . On a flat front at $p=0.066$ the anti-trapping model holds $1.004\,k$ , while the standard model at $p=0.043$ already holds $1.31\,k$ . At the tip’s Péclet number the standard model traps a lot.
Repeat the parameter calculation for Karma’s finer ratio $d_0/W=0.544$ . By how much does the unit of time change?
Show the answer
$\lambda=0.8839/0.544=1.625$ , $D=0.6267\times1.625=1.018$ , and $d_0^2/D=0.544^2/1.018=0.291\,\tau$ , about $7.6$ times longer than at $0.277$ . A run to the same scaled time therefore needs many more steps.
Local equilibrium at a curved tip predicts that the solid holds $c_s/c_l^0=k[1-(1-k)d_0/\rho]$ . Evaluate it for a tip radius $\rho=21.3\,d_0$ .
Show the answer
$0.15\times(1-0.85/21.3)=0.15\times0.960=0.144$ . The route benchmark measures $0.1447$ along the axis of the anti-trapping dendrite with that tip radius.

Sources and what we recomputed

A. Karma, Phase-field formulation for quantitative modeling of alloy solidification, arXiv:cond-mat/0103289v1 (2001)
Version 1, full text. Eqs. (16)–(18) and (22): the isothermal dilute-alloy model with the anti-trapping current. Eqs. (19)–(20): the equilibrium profile and the choice of $q(\phi)$ that removes surface diffusion. The parameter paragraph: $\Delta x=0.4$ , $\Delta t=0.008$ , $k=0.15$ , $\Omega=0.55$ , $\epsilon_4=0.02$ , $a_1=0.8839$ , $a_2=0.6267$ , a seed of radius $22\,d_0$ . Fig. 1: tip velocity against time. Fig. 2: solute in the solid along the dendrite axis.
B. Echebarria, R. Folch, A. Karma and M. Plapp, Quantitative phase-field model of alloy solidification, arXiv:cond-mat/0404164v2 (2004)
Version 2. Section I B, Eq. (2): with a mesoscopic interface thickness, solute trapping, surface diffusion and interface stretching all grow in proportion to $W$ . Section I D: the anti-trapping current as a non-variational fix. Section I E, Eqs. (7)–(13): the full model, the current of Eq. (10) with $a=1/(2\sqrt2)$ , and the relations $d_0=a_1W/\lambda$ and $\beta$ of Eq. (13). Section II A, Eqs. (17)–(21): the sharp-interface problem for isothermal growth, with the Stefan condition and the Gibbs–Thomson relation. Section IV gives the thin-interface analysis; it was skimmed, not rederived.