I asked for σ₁₂ = 1.00 and got 0.78 — where an absent phase leaks into the interface in three-component Allen–Cahn

I asked for σ₁₂ = 1.00 and got 0.78#

I put three phases into a phase-field code and typed 1.00 into the 1-2 surface tension field. Then I integrated the free energy across the relaxed interface. The answer was 0.78.

Phase 3 was the one that should not have mattered. It was nowhere in the domain. I had seeded it at $10^{-3}$ and nothing else. Yet at equilibrium its volume fraction peaked at 0.21 right in the middle of the 1-2 interface.

This post tracks that 21% all the way down. We invert the pairwise $\sigma_{ij}$ into per-component gradient coefficients $\kappa_i$, show that the sign of $\kappa_i$ is the same inequality that closes the Neumann triangle, and find the one term that keeps an absent phase out of the interface. At the end, four closures run on the same grid and the numbers decide.

N volume fractions are not N degrees of freedom#

With $N$ immiscible fluids, the volume fractions $c_i$ carry one constraint.

That leaves $N-1$ degrees of freedom. But code allocates $N$ arrays, and attaching one Allen–Cahn equation to each array is the obvious thing to do. Copy a two-phase solver and bump the component count, and this is exactly what you get.

The sum then breaks. Write each chemical potential $\mu_i$ independently and advance with $\partial_t c_i = -M \mu_i$, and you get

with no reason for the right-hand side to vanish. In the run below $\left|\sum c_i - 1\right|$ opens up to 0.859, the interface dissolves, and the surface tension drops to 0.026.

The textbook fix is a Lagrange multiplier. Subtract $\bar\mu = \frac{1}{N}\sum_j \mu_j$,

and the sum is conserved exactly. Swap between the two closures yourself in the simulation below.

Pick 3 equations, no closure and the $\sum c_i - 1$ trace collapses downward while all three curves smear out. Switch to mean projection and the trace pins to $10^{-13}$, but a red $c_3$ bump rises at the center of the interface and sigma_eff sits below 1.00. Push the Lambda slider up and that bump gets flattened.

Inverting pairwise σᵢⱼ into per-component κᵢ#

The free energy is usually written pairwise.

$\beta_{ij}$ is the depth of the double well for the $i$-$j$ pair and $\kappa_i$ is the gradient energy coefficient of component $i$. The wells are pairwise; the gradient coefficients are per component. That asymmetry causes everything that follows.

Neither $\beta_{ij}$ nor $\kappa_i$ is yours to choose freely. The two-phase limit ties them to physical quantities. Substituting $c_2 = 1 - c_1$ gives

which matches the familiar two-phase form $\beta c^2(1-c)^2 + \frac{\lambda}{2}|\nabla c|^2$ with $\lambda = \kappa_1 + \kappa_2$. Integrating the 1D equilibrium $\lambda c'' = dE_0/dc$ once gives $\frac{\lambda}{2}c'^2 = E_0$, that is $c' = \sqrt{2\beta/\lambda}\, c(1-c)$. The solution is a tanh, and

where $D$ is the interface thickness in lattice cells and $\sigma$ the surface tension. Inverting both relations,

$\beta_{ij}$ falls out pair by pair. $\kappa_i$ needs the three equations solved together. For three phases,

One per-component coefficient is a combination of three pairwise tensions. This is the central identity of the post. In the two-phase conservative Allen–Cahn post on mobility and the anti-diffusion term there was a single $\lambda$. With three phases that $\lambda$ splits into three, and it splits by subtraction.

With $D=4$ the prefactor $\frac{3D}{4}$ is $3$. On the second row of the table, $\sigma=(1.00,\,0.80,\,0.60)$ gives $\kappa_1=3\times1.20=3.60$, $\kappa_2=3\times0.80=2.40$, and $\kappa_3=3\times0.40=1.20$. The three coefficients sum to $7.20$. The $1$-$2$ pair alone sums to $6.00$. The two-phase relation $\kappa_i+\kappa_j=\frac{3}{2}\sigma_{ij}D$ with $1.00$ and $4$ also returns $6.00$. Put $\lambda=6$ into $\sigma=\frac{2\lambda}{3D}$ and the quotient is $\frac{12}{12}=1.00$. The continuum formula matches the typed input. Case A on the grid returns $0.9665$, and that $3.3\%$ gap is what four cells across the interface cost. The same factor $3$ on $\sigma=(1.00,\,0.40,\,0.40)$ gives $\kappa_3=-0.60$. That is the row where the triangle refuses to close. At $(1.00,\,0.55,\,0.55)$, $\kappa_3=0.30$ is still positive and the wedges are $155.4^\circ$, $155.4^\circ$, $49.2^\circ$. The margin before the sign dies is $0.30$, and the phase-3 angle is already down to $49^\circ$. $\kappa_3=1.20$ never enters the $\sigma_{12}$ budget. It shows up as gradient energy only after phase 3 occupies the interface. Equal tensions $\sigma=(1.00,\,1.00,\,1.00)$ put every $\kappa$ at $3.00$ and every wedge at $120^\circ$. The sum is $360^\circ$. The $49.2^\circ$ wedge on the $0.55$ row is $70.8^\circ$ narrower than that $120^\circ$. The angle pinches before the sign flips.

Where κ₃ crosses into negative#

Subtraction means the sign can be lost. If $\sigma_{13} + \sigma_{23} < \sigma_{12}$ then $\kappa_3 < 0$ and the gradient energy opens downward. Sharpening the interface now lowers the energy, so the solution collapses to the grid spacing. The model does not merely degrade: the functional is unbounded below from the start.

The same condition is known under another name. At a triple junction the three tensions balance and fix the wedge angle $\theta_i$ of each phase.

$\theta_1$ is the angle occupied by phase 1; rotate the indices for the other two. The condition for this cosine to leave $[-1, 1]$ is exactly the condition $\kappa_i < 0$. In other words, whether the Neumann triangle closes and whether the model is well posed are the same inequality. Hold $\sigma_{12}=1.00$ and let $\sigma_{13}+\sigma_{23}$ fall below $\sigma_{12}$: the triangle breaks and the phase-3 wedge opens to 180°. That is phase 3 wetting the 1-2 interface completely. At $\sigma=(1.00, 0.80, 0.60)$ the wedges come out as 143.1° / 126.9° / 90.0°. A 3-4-5 right triangle, so the phase-3 angle is exactly a right angle.

Four closures on the same grid#

A 1D grid of 200 cells, $D=4$, $\sigma=(1.00, 0.80, 0.60)$. Phase 1 on the left, phase 2 on the right, phase 3 seeded at $10^{-3}$ in the interface only. After 8000 relaxation steps I measure $\left|\sum c_i - 1\right|$, $\max c_3$, and the free-energy integral across the interface, which is the realized surface tension.

import math
import numpy as np
 
D, MOB, DT, NX = 4.0, 1.0 / 3.0, 0.2, 200
 
 
def pair_to_component(s12, s13, s23):
    """Pairwise surface tensions -> per-component gradient coefficients kappa_i and wells beta_ij."""
    f = 0.75 * D
    kap = np.array([f * (s12 + s13 - s23),
                    f * (s12 + s23 - s13),
                    f * (s13 + s23 - s12)])
    bet = 12.0 / D * np.array([[0, s12, s13], [s12, 0, s23], [s13, s23, 0]])
    return kap, bet
 
 
def neumann_wedges(s12, s13, s23):
    """Wedge angles of the three phases in degrees if the Neumann triangle closes, else None."""
    args = [(s23 ** 2 - s12 ** 2 - s13 ** 2) / (2 * s12 * s13),
            (s13 ** 2 - s12 ** 2 - s23 ** 2) / (2 * s12 * s23),
            (s12 ** 2 - s13 ** 2 - s23 ** 2) / (2 * s13 * s23)]
    if any(abs(a) > 1 for a in args):
        return None
    return [math.degrees(math.acos(a)) for a in args]
 
 
def _ext(a):
    return np.concatenate(([a[0]], a, [a[-1]]))
 
 
def lap1d(a):
    e = _ext(a)
    return e[:-2] - 2 * a + e[2:]
 
 
def grad1d(a):
    e = _ext(a)
    return (e[2:] - e[:-2]) / 2.0
 
 
def chem_potential(c, kap, bet, lam3):
    c1, c2, c3 = c
    dbulk = np.array([
        2 * bet[0, 1] * c1 * c2 ** 2 + 2 * bet[0, 2] * c1 * c3 ** 2 + 2 * lam3 * c1 * (c2 * c3) ** 2,
        2 * bet[0, 1] * c2 * c1 ** 2 + 2 * bet[1, 2] * c2 * c3 ** 2 + 2 * lam3 * c2 * (c1 * c3) ** 2,
        2 * bet[0, 2] * c3 * c1 ** 2 + 2 * bet[1, 2] * c3 * c2 ** 2 + 2 * lam3 * c3 * (c1 * c2) ** 2,
    ])
    return np.array([dbulk[i] - kap[i] * lap1d(c[i]) for i in range(3)])
 
 
def line_tension(c, kap, bet, lam3):
    """Free-energy integral across the interface = the realized surface tension."""
    c1, c2, c3 = c
    e = (bet[0, 1] * (c1 * c2) ** 2 + bet[0, 2] * (c1 * c3) ** 2
         + bet[1, 2] * (c2 * c3) ** 2 + lam3 * (c1 * c2 * c3) ** 2)
    for i in range(3):
        e = e + 0.5 * kap[i] * grad1d(c[i]) ** 2
    return float(e.sum())
 
 
def initial_profile(seed):
    x = np.arange(NX) - NX / 2
    c3 = seed * np.exp(-(x / D) ** 2)
    c1 = 0.5 * (1 - np.tanh(2 * x / D)) * (1 - c3)
    c2 = 0.5 * (1 + np.tanh(2 * x / D)) * (1 - c3)
    return np.array([c1, c2, c3])
 
 
def relax(kap, bet, closure, lam3=0.0, steps=8000, seed=1e-3):
    c = initial_profile(0.0 if closure == 'two-phase' else seed)
    for _ in range(steps):
        mu = chem_potential(c, kap, bet, lam3)
        if closure == 'two-phase':
            m = 0.5 * (mu[0] + mu[1])
            rhs = np.array([mu[0] - m, mu[1] - m, np.zeros(NX)])
        elif closure == 'independent':
            rhs = mu
        else:
            rhs = mu - mu.mean(axis=0)
        c = c - MOB * DT * rhs
        if not np.isfinite(c).all() or np.abs(c).max() > 1e3:
            return None, None, None
    return c, float(np.abs(c.sum(axis=0) - 1.0).max()), line_tension(c, kap, bet, lam3)
 
 
print("== pairwise sigma_ij -> per-component kappa_i (D = %.1f) ==" % D)
print("%-20s %-26s %-8s %s" % ("sigma_12,13,23", "kappa_1,2,3", "closes", "wedge angles [deg]"))
for s in [(1.0, 1.0, 1.0), (1.0, 0.8, 0.6), (1.0, 0.55, 0.55), (1.0, 0.4, 0.4)]:
    kap, _ = pair_to_component(*s)
    th = neumann_wedges(*s)
    print("%-20s %-26s %-8s %s" % (
        "%.2f %.2f %.2f" % s,
        "%+.2f %+.2f %+.2f" % tuple(kap),
        "yes" if kap.min() > 0 else "NO",
        "-" if th is None else "%.1f %.1f %.1f (sum %.1f)" % (th[0], th[1], th[2], sum(th))))
 
kap, bet = pair_to_component(1.0, 0.8, 0.6)
print()
print("== 1-2 interface at equilibrium, phase 3 seeded at 1e-3 (input sigma_12 = 1.00) ==")
print("%-30s %-13s %-9s %s" % ("closure", "max|sum c-1|", "max c_3", "sigma_eff"))
cases = [("A two-phase (c_3 frozen)", 'two-phase', 0.0),
         ("B three equations, no closure", 'independent', 0.0),
         ("C mean projection", 'projected', 0.0),
         ("D mean projection + Lambda=400", 'projected', 400.0)]
for tag, closure, lam3 in cases:
    c, drift, sig = relax(kap, bet, closure, lam3)
    if c is None:
        print("%-30s diverged" % tag)
        continue
    print("%-30s %-13.2e %-9.4f %.4f" % (tag, drift, c[2].max(), sig))
 
print()
print("== triple-penalty Lambda scan (mean projection) ==")
print("%-10s %-10s %-10s %s" % ("Lambda", "max c_3", "sigma_eff", "vs input"))
for lam3 in (0.0, 10.0, 30.0, 60.0, 150.0, 400.0):
    c, drift, sig = relax(kap, bet, 'projected', lam3)
    print("%-10.0f %-10.4f %-10.4f %+.1f%%" % (lam3, c[2].max(), sig, (sig - 1.0) * 100))

== pairwise sigma_ij -> per-component kappa_i (D = 4.0) ==
sigma_12,13,23       kappa_1,2,3                closes   wedge angles [deg]
1.00 1.00 1.00       +3.00 +3.00 +3.00          yes      120.0 120.0 120.0 (sum 360.0)
1.00 0.80 0.60       +3.60 +2.40 +1.20          yes      143.1 126.9 90.0 (sum 360.0)
1.00 0.55 0.55       +3.00 +3.00 +0.30          yes      155.4 155.4 49.2 (sum 360.0)
1.00 0.40 0.40       +3.00 +3.00 -0.60          NO       -
 
== 1-2 interface at equilibrium, phase 3 seeded at 1e-3 (input sigma_12 = 1.00) ==
closure                        max|sum c-1|  max c_3   sigma_eff
A two-phase (c_3 frozen)       1.78e-15      0.0000    0.9665
B three equations, no closure  8.59e-01      0.0000    0.0257
C mean projection              2.59e-13      0.2114    0.7809
D mean projection + Lambda=400 8.77e-15      0.0512    0.9715
 
== triple-penalty Lambda scan (mean projection) ==
Lambda     max c_3    sigma_eff  vs input
0          0.2114     0.7809     -21.9%
10         0.1887     0.8031     -19.7%
30         0.1448     0.8377     -16.2%
60         0.0963     0.8704     -13.0%
150        0.0515     0.9111     -8.9%
400        0.0512     0.9715     -2.9%

Moving the three $\sigma$ sliders flips the sign of the $\kappa$ bars and the wedge angles together. Drop them to $(1.00,\,0.40,\,0.40)$ and the closes column in the table turns over.

Hold $\sigma_{12}$ at $1.00$ and lower $\sigma_{13}$ and $\sigma_{23}$. The triangle breaks the moment the red $\kappa_3$ crosses $0$. That is the $\kappa_3=-0.60$ row at $(1.00,\,0.40,\,0.40)$. At $(1.00,\,0.80,\,0.60)$ the angles stay $143.1^\circ$ / $126.9^\circ$ / $90.0^\circ$.

Σc holds to 1e-13 while the interface loses 21.9%#

A freezes phase 3 outright, so it is a two-phase run. $\sigma_{\rm eff} = 0.9665$, which is 3.3% below the input 1.00. That is not a model defect but discretization error at four cells across the interface. It confirms the relation $\kappa_1 + \kappa_2 = \frac{3}{2}\sigma_{12}D$ holds.

B solves the three equations independently. The sum drifts by 0.859 and the surface tension vanishes to 0.0257. The interface melted.

C is the interesting line. The mean projection pins the sum at $2.6\times10^{-13}$. But $\max c_3 = 0.2114$ and $\sigma_{\rm eff} = 0.7809$, which is 21.9% below the input. The constraint is satisfied and the physics is wrong.

A constraint violation announces itself. A 21.9% error does not. That makes it the expensive one.

The arithmetic of an absent phase leaking in#

Strip away the gradient term and the leak is obvious. Assume symmetry at the center of the 1-2 interface and set $c_3 = p$, $c_1 = c_2 = (1-p)/2$.

Differentiate in $p$ and evaluate at $p=0$. The second term starts at $p^2$, so it contributes nothing to the slope. Only the first term survives.

Always negative, no matter what $\sigma_{13}$ and $\sigma_{23}$ are. In a pairwise quartic well, adding a third phase to a binary interface is always downhill. With $\beta_{12} = 3.0$ the slope is $-0.75$.

The numbers agree. At $p=0$, $E_0 = \beta_{12}/16 = 0.1875$. At the measured $p = 0.2114$, $E_0 = 0.1018$, a 45.7% drop. The bulk energy welcomes phase 3. The only thing holding it back is $\kappa_3 |\nabla c_3|^2$, and with $\kappa_3 = 1.20$ the balance only arrives at 0.21.

There is a family resemblance to parasitic currents, where a drop that should sit still starts to flow. When the tension discretization is off, the sum can be right and the interface quietly wrong. Here $\sum c_i$ sits at $2.59\times10^{-13}$ while $\sigma_{\rm eff}$ is $0.7809$.

Where the triple-well penalty Λ belongs#

The remedy is a term that penalizes three phases overlapping at one point.

In the scan, raising $\Lambda$ pushes $\max c_3$ from 0.2114 down to 0.0512 and brings $\sigma_{\rm eff}$ back from 0.7809 to 0.9715. At $\Lambda = 400$ the deficit is −2.9%, essentially the same place as A's discretization error of −3.3%. It is indistinguishable from freezing phase 3.

The term does not plug the hole, though. Differentiating $\Lambda c_1^2c_2^2c_3^2$ in $p$ also gives zero at $p=0$. The downhill slope $-\beta_{12}/4$ is still there. What $\Lambda$ does is build a wall after the bump has already grown, which is why $c_3$ stops at 0.05 instead of going to zero however high you push $\Lambda$.

To close the hole itself, reduce to $N-1$ degrees of freedom with $c_N = 1 - \sum_{i<N} c_i$, or simply do not solve the equation of a locally absent phase. A is that case, and that is why its $\max c_3$ is exactly zero.

Whether to trust the typed 1.00 or the integrated 0.7809#

Typing the pairwise table freezes the sign of $\kappa_i = \frac{3D}{4}(\sigma_{ij}+\sigma_{ik}-\sigma_{jk})$ and the wedge angles together. At $\sigma=(1.00,\,0.80,\,0.60)$ the coefficients are $+3.60$, $+2.40$, $+1.20$ and the angles are $143.1^\circ$ / $126.9^\circ$ / $90.0^\circ$, summing to $360^\circ$. At $\sigma=(1.00,\,0.40,\,0.40)$, $\kappa_3=-0.60$ and the triangle does not close. Five phases make ten such triples. One negative entry and the energy opens downward in the cell where that triple meets. Asking whether the three $\cos\theta_i$ sum to $360^\circ$ is the same question as asking whether every $\kappa_i$ is positive. An angle can be laid on the picture.

The number to verify against is not the typed $\sigma_{12}=1.00$. It is C's integral, $0.7809$. The sum holds at $2.59\times10^{-13}$ while the tension is $21.9\%$ low. As with a capillary time step that remains after surface tension is treated implicitly, that $21.9\%$ rides straight into the Laplace pressure and the capillary-wave speed. A check that only watches $\sum c_i=1$ passes. Even at $\Lambda=400$ the deficit is $-2.9\%$, or $0.9715$, which differs from A's frozen-phase value $0.9665$ by $0.0050$. Leave $1.00$ in the input box as the reference and that gap never reaches the log.