The ACMultiInterface kernel implements a multiphase interface gradient term with cross terms between all order parameters from Nestler, B., and A. A. Wheeler. “Anisotropic Multi-Phase-Field Model: Interfaces and Junctions.” Physical Review E 57, no. 3 (1998): 2602. Eq. (7) and (8) (also see footnote 1)

f_{int} = \sum_{\substack{a,b \\ b\neq a}} \frac12 \kappa_{ab} \left| \eta_a\nabla\eta_b - \eta_b\nabla\eta_a\right|^2,

Where the sum is taken over unique tuples a,b (i.e. without the permutations b,a). We take the functional derivative taken using the lemma

\frac{\delta f}{\delta\eta} = \frac{\partial f}{\partial\eta} - \nabla\frac{\partial f}{\partial\nabla\eta}.

We obtain a one dimensional sum for each of the $$\eta$$-derivatives. \begin{aligned} \frac{\delta f_{int}}{\delta\eta_a} & = \sum_b \kappa_{ab} \left[ (\eta_a\nabla\eta_b - \eta_b\nabla\eta_a)\nabla\eta_b + \nabla\left((\eta_a\nabla\eta_b - \eta_b\nabla\eta_a)\eta_b\right) \right] \\ &= \sum_{\substack{b\\b\neq a}} \kappa_{ab} \left[ \underbrace{2(\eta_a\nabla\eta_b - \eta_b\nabla\eta_a)\nabla\eta_b}_{\text{order 1}} + \underbrace{\eta_b(\eta_a\nabla^2\eta_b - \eta_b\nabla^2\eta_a)}_{\text{order 2}} \right] \end{aligned} We transform this expression into the weak form and see that the derivative order on the order 2 term has to be reduced by shifting a gradient onto the test function by applying the product rule v \nabla\cdot\mathbf{w} = -\nabla v \cdot \mathbf{w} + \nabla\cdot (v\mathbf{w}), after multiplying with the test function $$\psi$$ and integrating over the volume $$\Omega$$$. We identify $$v$$$ and $$w$$$as follows \int_\Omega\psi\eta_b(\eta_a\nabla^2\eta_b - \eta_b\nabla^2\eta_a) = \int_\Omega\underbrace{\psi\eta_a\eta_b}_{v_1}\nabla\cdot\underbrace{\nabla\eta_b}_{\mathbf{w}_1} - \int_\Omega\underbrace{\psi\eta_b^2}_{v_2}\nabla\cdot\underbrace{\nabla\eta_a}_{\mathbf{w}_2} \int_\Omega\left[ -\nabla(\psi\eta_a\eta_b)\cdot\nabla\eta_b\right] - \int_\omega\left[-\nabla(\psi\eta_b^2)\cdot\nabla\eta_a\right] + \int_\Omega \nabla\cdot(\psi\eta_a\eta_b\nabla\eta_b) - \int_\Omega \nabla\cdot(\psi\eta_b^2\nabla\eta_a). We get rid of the last two terms by applying the divergence theorem and obtain \underbrace{ \int_\Omega\left[ -\nabla(\psi\eta_a\eta_b)\cdot\nabla\eta_b\right] - \int_\Omega\left[-\nabla(\psi\eta_b^2)\cdot\nabla\eta_a\right] }_{\text{volume terms}} + \underbrace{ \int_{\partial\Omega} \mathbf{n}\cdot(\psi\eta_a\eta_b\nabla\eta_b) - \int_{\partial\Omega} \mathbf{n}\cdot(\psi\eta_b^2\nabla\eta_a). }_{\text{boundary terms}} to convert them from volume to surface/boundary integrals. We again apply the product rule to expand the gradient of the product in the volume terms and obtain \int_\Omega\left[ -\left( \eta_a\eta_b\nabla\psi + \psi\eta_b\nabla\eta_a + \psi\eta_a\nabla\eta_b \right) \cdot\nabla\eta_b\right] - \int_\Omega\left[ -\left( \eta_b^2\nabla\psi + 2\psi\eta_b\nabla\eta_b \right)\cdot\nabla\eta_a\right]. ## Residual The total residual $$R_a$$$ is then

\begin{aligned} R_a = L_a\sum_{\substack{b\\b\neq a}}\kappa_{ab} \Bigg[ &\,& \int_\Omega\left[ 2\psi(\eta_a\nabla\eta_b - \eta_b\nabla\eta_a)\nabla\eta_b \right]& \\ &+& \int_\Omega\left[ -\left( \eta_a\eta_b\nabla\psi + \psi\eta_b\nabla\eta_a + \psi\eta_a\nabla\eta_b \right) \cdot\nabla\eta_b\right] & \\ &-& \int_\Omega\left[ -\left( \eta_b^2\nabla\psi + 2\psi\eta_b\nabla\eta_b \right)\cdot\nabla\eta_a\right]& \Bigg]. \end{aligned}

## On-diagonal Jacobian

The on-diagonal jacobian $$J_a$$$is obtained by taking the derivative with respect to $$\eta_{aj}$$$, where $$\frac{\partial \eta_a}{\partial \eta_{aj}} = \phi_j$$$and $$\frac{\partial \nabla\eta_{aj}}{\partial \eta_{aj}} = \nabla\phi_j$$$

\begin{aligned} J_a = L_a\sum_{\substack{b\\b\neq a}}\kappa_{ab} \Bigg[ &\,&\int_\Omega2\psi\left[ (\phi_j\nabla\eta_b - \eta_b\nabla\phi_j)\nabla\eta_b \right]& \\ &+& \int_\Omega\left[ -\left( \phi_j\eta_b\nabla\psi + \psi\eta_b\nabla\phi_j + \psi\phi_j\nabla\eta_b \right)\cdot\nabla\eta_b \right]& \\ &-& \int_\Omega\left[ -\left( \eta_b^2\nabla\psi + 2\psi\eta_b\nabla\eta_b \right)\cdot\nabla\phi_j \right]& \Bigg]. \end{aligned}

## Off-diagonal jacobian

For the off diagonal Jacobian entry $$J_{ab}$$$we take the derivative $$\frac{\partial R_a}{\partial \eta_{bj}}$$$ and obtain

\begin{aligned} J_{ab} = &\,& L_a\kappa_{ab}\int_\Omega2\psi\left[ (\eta_a\nabla\phi_j - \phi_j\nabla\eta_a)\nabla\eta_b + (\eta_a\nabla\eta_b - \eta_b\nabla\eta_a)\nabla\phi_j \right] \\ &+& \int_\Omega\left[ -\left( \eta_a\phi_j\nabla\psi + \psi\phi_j\nabla\eta_a + \psi\eta_a\nabla\phi_j \right) \cdot\nabla\eta_b -\left( \eta_a\eta_b\nabla\psi + \psi\eta_b\nabla\eta_a + \psi\eta_a\nabla\eta_b \right) \cdot\nabla\phi_j \right] \\ &-& \int_\Omega\left[ %-\left( \eta_b^2\nabla\psi + 2\psi\eta_b\nabla\eta_b \right)\cdot\nabla\eta_a -\left( 2\eta_b\phi_j\nabla\psi + 2\psi(\phi_j\nabla\eta_b + \eta_b\nabla\phi_j) \right)\cdot\nabla\eta_a \right] \end{aligned}

1) Note, that in the two-phase case with $$\eta_b=1-\eta_a$$\$ this reduces to

\begin{aligned} f_{int} & = \frac12 \kappa_{ab} \left| \eta_a\nabla(1-\eta_a) - (1-\eta_a)\nabla\eta_a\right|^2 \\ & = \frac12 \kappa_{ab} \left| -\eta_a\nabla\eta_a - \nabla\eta_a + \eta_a\nabla\eta_a\right|^2 \\ & = \frac12 \kappa_{ab} \left| \nabla\eta_a \right|^2, \end{aligned}

which is the familiar form implemented by ACInterface.