Singularities in phase separation models: a spectral element approach for the nonlocal Cahn-Hilliard equation

1. Problem Statement

The paper addresses the numerical approximation of the nonlocal Cahn–Hilliard (CH) equation, a macroscopic model describing phase separation in binary mixtures (e.g., alloys). Unlike the classical local CH equation, the nonlocal variant incorporates long-range interactions via a convolution term with a singular kernel.

The specific system studied is:
$$\frac{\partial \rho}{\partial t} = \nabla \cdot \left( m(\rho) \nabla \mu \right), \quad \mu = F'(\rho) - K * \rho$$
where:

• $\rho$: Phase field variable (concentration).
• $m(\rho)$: Constant mobility (set to $m=1$).
• $F(\rho)$: A logarithmic potential (singular at pure states $\pm 1$), defined as $F(s) = \frac{\theta}{2}[(1+s)\log(1+s) + (1-s)\log(1-s)]$. This ensures the solution stays within the physical interval $(-1, 1)$.
• $K$: A Newtonian interaction kernel (in 2D, $K(x) = \frac{1}{2\pi}\log\|x\|$), which is singular at the origin but integrable.
• Boundary Conditions: No-flux ($\partial_n \mu = 0$) on a bounded domain $\Omega$.

Key Challenges:

1. Singularity: The kernel $K$ has a singularity at the origin, making standard collocation or finite difference methods difficult to apply directly without special treatment.
2. Non-locality: The convolution term $K * \rho$ results in dense operators, leading to high computational costs ($O(N^2)$) for standard discretizations.
3. Singular Potential: The derivative $F'(\rho)$ blows up as $\rho \to \pm 1$, requiring strict numerical control to prevent the solution from leaving the physical domain.

2. Methodology

The authors propose a pseudospectral multishape method combined with a cutoff decomposition technique to handle the singular convolution efficiently.

A. Mathematical Framework & Existence Theory

Before numerical implementation, the authors extend the theoretical analysis of the system:

• Domain Generalization: They prove existence and uniqueness of weak solutions on Lipschitz domains (including polygonal domains with corners), rather than just smooth domains.
• Strict Separation: They utilize the "strict separation property" to prove that solutions stay uniformly away from the pure phases ($\pm 1$) for any $t > 0$, even in 3D.
• Long-time Behavior: They establish convergence rates to a unique equilibrium state using the Łojasiewicz–Simon inequality, showing polynomial decay in the $L^2$ norm.

B. Numerical Approximation of the Singular Convolution

To compute $(K * \rho)(x)$, the integral is split into two parts using a cutoff radius $\varepsilon$:
$$(K * \rho)(x) = \underbrace{\int_{\Omega \setminus B_\varepsilon(x)} K(x-y)\rho(y) dy}_{I_1} + \underbrace{\int_{\Omega \cap B_\varepsilon(x)} K(x-y)\rho(y) dy}_{I_2}$$

1. Handling the Singularity ($I_2$):

   • The inner integral over the small ball $B_\varepsilon(x)$ is approximated using a zero-order Taylor expansion (piecewise constant approximation):
     $$I_2 \approx \rho(x) \int_{x - (\Omega \cap B_\varepsilon(x))} K(u) du$$
   • The authors derive closed-form analytical expressions for the integral of the Newtonian kernel over these geometric regions (using polar coordinates and specific geometric decompositions). This term acts as a correction factor and an error indicator.
   • The error of this approximation is shown to be $O(\varepsilon^2 \log \varepsilon)$.

2. Handling the Regular Part ($I_1$):

   • The outer integral (over the domain excluding the singularity) is computed using a Spectral Element Method (SEM).
   • They utilize the MultiShape software package, which constructs complex domains by patching together quadrilateral elements.
   • For each collocation point $x_i$, the domain $\Omega \setminus B_\varepsilon(x_i)$ is decomposed into a "multishape" of quadrilaterals.
   • The convolution is precomputed as a matrix $M$ acting on the discretized $\rho$. This matrix is computed once (since $K$ is time-independent) and reused, significantly reducing the cost of time-stepping.

C. Time Integration

• The resulting semi-discrete system is solved using a stiff solver (ode15s in MATLAB) with variable step sizes.
• To handle the stiffness introduced by the singular potential $F'$ near pure states, the authors also test a regularized potential ($F_\omega$) that matches the Taylor expansion of $F$ near the boundaries, allowing for larger time steps in long-time simulations.

3. Key Contributions

1. Theoretical Extension: Proved well-posedness and convergence to equilibrium for the nonlocal CH system on Lipschitz (polygonal) domains in 2D and 3D, extending previous results limited to smooth domains.
2. Novel Numerical Scheme: Developed a spectral element approach for singular kernels that combines:
   • Analytical integration of the singular part (via Taylor expansion).
   • High-order spectral integration of the regular part (via multishape decomposition).
   • This is the first high-resolution computation of the nonlocal CH system with a singular kernel and logarithmic potential under no-flux conditions.
3. Error Control: Provided rigorous error estimates showing the method achieves spectral convergence for the regular part and $O(\varepsilon^2 \log \varepsilon)$ for the singular correction, verified numerically.
4. Open Source: Released the code and data (GitHub: DDFT-Modelling/NLCH) to facilitate reproducibility.

4. Numerical Results

The authors performed extensive simulations on a unit square domain ($\Omega = (0,1)^2$) with the Newtonian kernel $K_\eta = \eta K$.

• Parameter $\eta$ (Scaling):
  • $\eta > 0$ (Repulsive): The system behaves diffusively, converging to a homogeneous state ($\rho \equiv 0$).
  • $\eta < 0$ (Attractive): The system exhibits phase separation.
    • For moderate negative $\eta$, slow phase separation occurs.
    • For large negative $\eta$, fast phase separation leads to sharp interfaces and step-function-like profiles.
• Initial Conditions:
  • Tested periodic waves and compactly supported functions.
  • Observed that the system evolves towards stationary states with distinct patterns (e.g., concentration at corners or formation of step functions).
• Regularization: Demonstrated that using the regularized potential $F_\omega$ allows for significantly faster computation of long-time stationary states (reaching machine precision with fewer iterations) without altering the physical outcome, as the solution naturally stays away from the singularities.
• Kernel Mixture: Tested a linear combination of the Newtonian kernel and a mollifier. This created a kernel with mixed positive/negative entries, inducing complex, pattern-rich structures not seen with pure Newtonian kernels.

5. Significance

• Computational Efficiency: The proposed method allows for high-resolution solutions with limited computational resources. By precomputing the convolution matrix and using analytical corrections for the singularity, it avoids the $O(N^2)$ cost of dense matrix-vector multiplications at every time step.
• Physical Accuracy: The use of the logarithmic potential ensures the physical constraint $|\rho| < 1$ is respected, which is crucial for modeling binary alloys and biological systems (e.g., tumor growth models).
• Versatility: The "multishape" approach is flexible and can be adapted to various domains (e.g., disks, cubes) and different singular kernels.
• Future Applications: The authors highlight the potential to apply this framework to coupled systems, such as the Cahn–Hilliard–Navier–Stokes equations or the Hawkins–Daruud model for tumor growth, where nonlocal interactions and singular potentials are critical.

In summary, this paper bridges the gap between rigorous mathematical analysis of singular nonlocal PDEs and practical, high-fidelity numerical simulation, providing a robust tool for studying complex phase separation phenomena.