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].$$$$


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.