A Convergent BDF2-SAV Scheme for Fractional Cahn-Hilliard Equations with High-Dimensional Integral Fractional Laplacian

1. Introduction

Phase-field models provide a standard baseline for simulating various free boundary problems [1]. Back in 1958, Cahn and Hilliard built the CH equation by minimizing a free energy functional to track the thermodynamic and kinetic evolution of binary alloys [1]. Evaluating solid-liquid interfacial tensions represents a core application, since this physical quantity largely controls solidification and phase separation inside porous media [2]. Modern imaging tasks, specifically binary image inpainting, also borrow this framework [3]. A continuous scalar phase field tracks the interface implicitly. That variable stays almost flat inside the bulk region but changes sharply across the diffuse layer [1]. To set up our current study, we define the coupled system on a high-dimensional bounded domain $\Omega \subset {\mathbb{R}}^d$ ($d\ge 2$) with homogeneous Dirichlet boundary conditions.

Physically, particle interactions at long ranges naturally unfold as spatial convolutions or non-local integrals. Squeezing these global dynamics into purely localized derivatives exposes a fundamental flaw inherent in traditional integer-order CH formulations. That arbitrary truncation exacts a severe physical toll [1]. By their very nature, fractional-order frameworks step in to faithfully map out these exact non-local transport mechanisms, recovering the global interactions that standard models filter out [1]. The Riesz fractional operator ${\mathcal{L}}_{\alpha }$ directly handles our spatial discretization.

$$\left\{\begin{array}{cc}{\delta }_{t}u={\mathcal{L}}_{\alpha }v,\hfill & \mathbf{x}\in \Omega ,t\in (0,T],\hfill \\ v=-{\epsilon }^2{\mathcal{L}}_{\alpha }u+f\left(u\right),\hfill & \mathbf{x}\in \Omega ,t\in (0,T],\hfill \\ u=v=0,\hfill & \mathbf{x}\in \delta \Omega ,t\in (0,T],\hfill \\ u(\mathbf{x},0)=u_0\left(\mathbf{x}\right),\hfill & \mathbf{x}\in \Omega ,\hfill \end{array}\right.$$

(1)

In the above system, the positive parameter $\epsilon$ controls the thickness of the diffuse interface [4], and the nonlinear function $f\left(u\right)=F^{\prime }\left(u\right)$ denotes the derivative of the macroscopic double-well potential $F\left(u\right)={\displaystyle \frac{1}{4}}{(u^2-1)}^2$ [5].

From a mathematical perspective, this coupled system represents the gradient flow of the total free energy $E\left(u\right)={\int }_{\Omega }\left({\displaystyle \frac{{\epsilon }^2}{2}}{\left|{\mathcal{L}}_{\alpha }^{1/2}u\right|}^2+F\left(u\right)\right)d\mathbf{x}$ evaluated in the dual Sobolev space [6]. Given the strictly dissipative nature of the CH framework, its continuous solution inherently obeys an energy dissipation law, meaning the global energy of the system decreases monotonically as time evolves [7]. Rigorously inheriting this dissipative behavior at the fully discrete level—which corresponds to developing unconditionally energy-stable algorithms—is a fundamental prerequisite for reliably simulating long-time macroscopic phase evolutions [8].

In the pursuit of maintaining energy stability while discretizing the stiff nonlinearities, various strategies have been proposed, including the classical convex splitting technique [9] and the invariant energy quadratization (IEQ) methodology [8]. Widespread attention currently surrounds the scalar auxiliary variable (SAV) formulation and its generalized extensions (gPAV) as highly efficient alternatives [10]. As a standard numerical tool for these gradient flows, the method introduces a single scalar variable to replace the square root of the nonlinear term. That substitution effectively linearizes the entire system, meaning one only needs to solve linear systems with constant coefficients at any arbitrary time step. First-order temporal accuracy, unfortunately, acts as a strict bottleneck for the majority of these current SAV-based algorithms. Extended simulation runs easily trigger unphysical high-frequency oscillations under such limited precision. Coupling the SAV approach with the second-order backward differentiation formula (BDF2) directly resolves this specific problem. It successfully pushes the temporal accuracy to second order while rigorously guaranteeing unconditional energy dissipation at the discrete level [7].

Scaling these high-order numerical schemes to multidimensional scenarios ($d\ge 2$) inevitably triggers critical bottlenecks in terms of computational speed and memory consumption [1] Discretization on the spatial grid restructures the fractional problem into an enormous linear algebraic system. A notoriously ill-conditioned, $2\times 2$ symmetric indefinite block matrix governs this setup [5]. Because fractional derivatives operate non-locally, the resulting coefficient matrices become completely dense. They naturally exhibit a multi-level block Toeplitz layout [12]. Inverting these blocks directly is computationally prohibitive. To bypass this obstacle, preconditioned Krylov subspace solvers step in as the mainstream workaround. Built on the discrete sine transform, the $\tau$-preconditioner compresses the spectral distribution. This action safely anchors the eigenvalues in a region away from zero [13]. Strict theoretical bounds cap the per-step complexity at $\mathcal{O}(NlogN)$ [5]. Practically, non-contiguous memory access patterns deep inside these high-dimensional tensor computations inevitably trigger an avalanche of severe cache misses. Under real hardware execution, those low-level penalties heavily dilute the actual algorithmic efficiency.

Against the dual demands of temporal accuracy and high-dimensional computational efficiency, we build a fully discrete algorithmic framework. Breakdowns of the core contributions follow:

• Second-order temporal discretization and unconditional stability. Hooking up the BDF2 method handles time stepping. That specific choice pushes the temporal precision straight to $\mathcal{O}\left({\tau }^2\right)$, natively supplying excellent stiff decay to damp out high-frequency wiggles. Formal proofs secure unconditional energy stability strictly at the discrete grid level.
• VTP-accelerated preconditioned fast solver. Spatially, a newly constructed sixth-order Riesz fractional centered difference scheme pairs with a sine-transform-based $\tau$-preconditioner. To break through the memory access bottlenecks inherent in high-dimensional tensor math, a vectorized tensor processing (VTP) layout steps in. Customized tensor shuffling allows this VTP approach to repackage preconditioned matrix-vector multiplications into strictly memory-contiguous batched transforms. Such low-level hardware alignment effectively cuts down on severe cache misses. While firmly preserving the $\mathcal{O}(NlogN)$ theoretical complexity bounds, actual computational throughput under high spatial resolutions experiences a massive boost.

The mathematical and computational framework of this study unfolds across four subsequent sections. Section 2 (Materials and Methods) constructs the discrete foundation by deriving the sixth-order Riesz fractional spatial scheme via double generating functions, explicitly verifying the structural properties of the weighting coefficients, and detailing the synergistic integration of the $\tau$-preconditioner with the continuous memory-access VTP architecture. Section 3 (Results) subjects this proposed framework to strict theoretical and empirical scrutiny, providing rigorous proofs for unconditional energy stability alongside multidimensional numerical validations of the $\mathcal{O}\left(h^6\right)$ spatial convergence and 3D phase separation morphologies under extreme fractional anisotropy. Moving to Section 4 (Discussion), we critically evaluate the physical and computational implications of these outcomes, isolating the memory access bottlenecks of traditional DST algorithms and demonstrating how tensor matricization within the VTP framework mitigates cache starvation to sustain optimal high-dimensional execution speedups. Finally, Section 5 (Conclusions) closes the manuscript by summarizing our central algorithmic contributions.

2. Materials and Methods

2.1. A Sixth-Order Approximation to the Riesz Fractional Derivative

To aggressively cut down spatial truncation errors, we introduce a specific generating function $G_{{\alpha }_i}\left(z\right)$ for the continuous Riesz fractional operator ${\mathcal{L}}_{x_i}^{{\alpha }_i}$ ($1<{\alpha }_i<2$).which is specifically constructed to eliminate the leading truncation errors in the Taylor expansion and achieve the desired sixth-order convergence rate. That mathematical setup explicitly wipes out the leading truncation errors during Taylor expansion to secure the sixth-order convergence rate.

That mathematical setup explicitly wipes out the leading truncation errors during Taylor expansion to secure the sixth-order convergence rate. Coefficient symmetry natively guarantees that the discrete fractional operator produces a symmetric Toeplitz-like matrix across the interior grid. Algebraically, this exact structure opens the door for spectral diagonalization using the discrete sine transform. Furthermore, its strict zero-sum property enforces the negative semi-definiteness of the discrete operator. Strict adherence to the discrete energy dissipation law naturally follows.

Let $\Omega$ be uniformly partitioned into square elements of side length h, and let M denote the total number of grid points. To obtain a high-order approximation of the Riesz derivative, we construct a numerical differentiation formula based on an appropriately designed generating function. Specifically, we introduce a new generating function $G\left(z\right)$ of the form

$$\begin{array}{c}\hfill g\left(z\right)={\displaystyle \frac{3}{2}}{(z^{-1}-2+z)}^{\alpha /2}-{\displaystyle \frac{3\left(2^{-\alpha }\right)}{5}}{(z^{-2}-2+z^2)}^{\alpha /2}+{\displaystyle \frac{3^{-\alpha }}{10}}{(z^{-3}-2+z^3)}^{\alpha /2},\end{array}$$

(2)

In summary, the discrete scheme obtained via the sixth-order generating function is presented as follows:

$$\begin{array}{c}\hfill {\delta }_h^{{\alpha }_i}u\left(x_i\right):=-{\displaystyle \frac{1}{h^{{\alpha }_i}}}\sum _{k=-\infty }^{\infty }{g}_k^{{\alpha }_i}u(x_i-k{h}_i),\left|\omega \right|\le 1\end{array}$$

(3)

Setting $z=e^{-i\theta }$ in (2) and applying the inverse Fourier transform, the explicit analytical expression for the corresponding weighting coefficients $g_k^{\left({\alpha }_i\right)}$ can be deduced as:

$$g_k^{{\alpha }_i}={\displaystyle \frac{3}{2}}{a}_k^{{\alpha }_i}-{\mathbf{1}}_{\left\{2\right|k\}}·{\displaystyle \frac{3\left(2^{-{\alpha }_i}\right)}{5}}{a}_{k/2}^{{\alpha }_i}+{\mathbf{1}}_{\left\{3\right|k\}}·{\displaystyle \frac{3^{-{\alpha }_i}}{10}}{a}_{k/3}^{{\alpha }_i},$$

(4)

where $\mathsf{\Gamma }(·)$ denotes the Euler gamma function and $a_k^{{\alpha }_i}={\displaystyle \frac{{(-1)}^k\mathsf{\Gamma }({\alpha }_i+1)}{\mathsf{\Gamma }({\alpha }_i/2-k+1)\mathsf{\Gamma }({\alpha }_i/2+k+1)}}$. Furthermore, we define a class of special functions as follows:

$${\mathbf{1}}_{\left\{k\right|n\}}=\left\{\begin{array}{cc}1,\hfill & n\equiv 0(modk);\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(5)

To ensure the energy stability of the subsequent fully discrete scheme and the applicability of the fast transform-based preconditioner, it is crucial to establish the algebraic properties of the coefficients $g_k^{\left({\alpha }_i\right)}$.

Proposition 1.  

For any fractional order ${\alpha }_i\in (1,2)$, the sixth-order weighting coefficients $g_k^{\left({\alpha }_i\right)}$ satisfy the following properties:

• ${\displaystyle 1<g_0^{{\alpha }_i}<3}$,
• ${\displaystyle g_k^{{\alpha }_i}=g_{-k}^{{\alpha }_i}<0,k=1,2,\dots }$,
• ${\displaystyle |g_{k+1}^{{\alpha }_i}|\le |g_k^{{\alpha }_i}|,k=1,2,\dots }$,
• ${\displaystyle \sum _{k=-\infty }^{\infty }{g}_k^{{\alpha }_i}=0,\sum _{k\ne 0}\left|g_k^{{\alpha }_i}\right|\le g_0^{{\alpha }_i}}$.

where $g_k^{\alpha }={\displaystyle \frac{{(-1)}^k\mathsf{\Gamma }(\alpha +1)}{\mathsf{\Gamma }(\alpha /2-k+1)\mathsf{\Gamma }(\alpha /2+k+1)}}$.

The symmetry property (i) guarantees that the discrete fractional operator yields a symmetric Toeplitz-like matrix in the interior domain, which is essential for the spectral diagonalization using the discrete sine transform. The zero-sum property (v) structurally ensures the semi-negative definiteness of the discrete operator, preserving the energy dissipation law at the discrete level.

Next, we take the sixth-order approximation to demonstrate that the format can achieve six-order accuracy. The remaining discrete formulas can be proved in a similar manner.

Proposition 2.  

Let

$$u\in {\mathcal{L}}^{6+\alpha }\left(\mathbb{R}\right)=\left\{u|u\in L^1\left(\mathbb{R}\right)and{\int }_{\mathbb{R}}{(1+|\omega \left|\right)}^{6+\alpha }\left|\widehat{u}\left(\omega \right)\right|d\omega <\infty \right\},1<\alpha <2,$$

then we have that

$${\delta }_{x,\mathbb{R}}^{\alpha }u\left(x\right)={\delta }_{h,\mathbb{R}}^{\alpha }u\left(x\right)+O\left(h^6\right),$$

(6)

uniformly in $x\in \mathbb{R}$.

Proof. 

Let we denote $g\left(i\varpi h\right)={\lambda }_1{\left(4sin\left({\displaystyle \frac{\varpi h}{2}}\right)\right)}^{\alpha }+{\lambda }_2{\left(4sin\left({\displaystyle \frac{2\varpi h}{2}}\right)\right)}^{\alpha }+{\lambda }_3{\left(4sin\left({\displaystyle \frac{3\varpi h}{2}}\right)\right)}^{\alpha }.$ and $u(\varpi ,h)=\mathcal{F}\left\{{\delta }_{x,\mathbb{R}}^{\alpha }u\left(x\right)-{\delta }_{h,\mathbb{R}}^{\alpha }u\left(x\right)\right\}.$ According to the properties of Fourier transform, there is

$$\mathcal{F}\left\{{\delta }_{x,\mathbb{R}}^{\alpha }u\right\}=-{\left|\varpi \right|}^{\alpha }\widehat{u}\left(\varpi \right),$$

and

$$\mathcal{F}\left\{{\delta }_{h,\mathbb{R}}^{\alpha }u\right\}=-{\left|\varpi \right|}^{\alpha }{\displaystyle \frac{g\left(i\varpi h\right)}{{(\varpi h/2)}^{\alpha }}}\widehat{u}\left(\varpi \right).$$

Together with

$$\mathcal{F}\left\{{\left(4sin\left({\displaystyle \frac{\varpi h}{2}}\right)\right)}^{\alpha }\right\}=-{\left|\varpi \right|}^{\alpha }\left(1-{\displaystyle \frac{\alpha {\left(\varpi h\right)}^2}{24}}+\left(-{\displaystyle \frac{\alpha }{2880}}+{\displaystyle \frac{{\alpha }^2}{1152}}\right){\left(\varpi h\right)}^4+O{\left(\varpi h\right)}^6\right),$$

we have

$$\begin{array}{cc}\hfill u(\varpi ,h)=& {\left|\varpi \right|}^{\alpha }\left\{-1+\left[({\lambda }_1+2^{\alpha }{\lambda }_2+3^{\alpha }{\lambda }_3)-{\displaystyle \frac{1}{24}}\left(\alpha ({\lambda }_1+2^{2+\alpha }{\lambda }_2+3^{2+\alpha }{\lambda }_3)\right){\left(\varpi h\right)}^2\right.\right.\hfill \\ & \left.\left.+{\displaystyle \frac{1}{5760}}\left(\alpha (-2+5\alpha )({\lambda }_1+2^{4+\alpha }{\lambda }_2+3^{4+\alpha }{\lambda }_3){\left(\varpi h\right)}^4\right)+O{\left(\varpi h\right)}^6\right]\right\}\widehat{u}\left(\varpi \right).\hfill \end{array}$$

Solving the equation

$$\left\{\begin{array}{c}{\lambda }_1+2^{\alpha }{\lambda }_2+3^{\alpha }{\lambda }_3=1,\hfill \\ \alpha ({\lambda }_1+2^{2+\alpha }{\lambda }_2+3^{2+\alpha }{\lambda }_3)=0,\hfill \\ \alpha (-2+5\alpha )({\lambda }_1+2^{4+\alpha }{\lambda }_2+3^{4+\alpha }{\lambda }_3=0,\hfill \end{array}\right.$$

yields

$${\lambda }_1={\displaystyle \frac{3}{2}},{\lambda }_2=-{\displaystyle \frac{3\left(2^{-\alpha }\right)}{5}},{\lambda }_3={\displaystyle \frac{3^{-\alpha }}{10}},$$

which in turn enables us to obtain

$$u(\varpi ,h)={\left|\varpi \right|}^{\alpha }\left[-{\displaystyle \frac{-16\alpha +42{\alpha }^2-35{\alpha }^3}{80640}}{\left(\varpi h\right)}^6+\mathcal{O}{\left(\varpi h\right)}^8\right]\widehat{u}\left(\varpi \right).$$

Thus, there exists a constant $C_{h1}>0$ such that

$$\left|u(\varpi ,h)\right|\le C_{h1}{\left|\varpi \right|}^{\alpha +6}\left|\widehat{u}\left(\varpi \right)\right|h^6.$$

From the inverse Fourier transform,

$$\begin{array}{cc}\hfill \left|{\delta }_{x,\mathbb{R}}^{\alpha }u\left(x\right)-{\delta }_{h,\mathbb{R}}^{\alpha }u\left(x\right)\right|& =\left|{\displaystyle \frac{1}{2\pi }}{\int }_{\mathbb{R}}u(\varpi ,h)e^{i\varpi x}d\varpi \right|\hfill \\ & \le {\displaystyle \frac{1}{2\pi }}{\int }_{\mathbb{R}}\left|u(\varpi ,h)\right|d\varpi \hfill \\ & \le {\displaystyle \frac{C_{h1}}{2\pi }}\left({\int }_{\mathbb{R}}{\left|\varpi \right|}^{\alpha +6}\left|\widehat{u}\left(\varpi \right)\right|d\varpi \right)h^6.\hfill \end{array}$$

Since $u\in C^{6+\alpha }\left(\mathbb{R}\right)$, we have${\int }_{\mathbb{R}}{(1+|\varpi |}^{\alpha +6}\left)\right|\widehat{u}\left(\varpi \right)|d\varpi <\infty$

therefore $\left|{\delta }_{x,\mathbb{R}}^{\alpha }u\left(x\right)-{\delta }_{h,\mathbb{R}}^{\alpha }u\left(x\right)\right|\le C_{h2}{h}^6,$ Here, $C_{h2}={\displaystyle \frac{C_{h1}}{2\pi }}{\int }_{\mathbb{R}}{(1+|\varpi |}^{\alpha +6}\left)\right|\widehat{u}\left(\varpi \right)|d\varpi$. □

We now extend this approximation to construct the multi-dimensional grid system. Consider a d-dimensional rectangular domain $\Omega ={\prod }_{i=1}^d(a_i,b_i)$. For each spatial direction $i\in \{1,2,\dots ,d\}$, let $M_i$ be a positive integer representing the number of interior grid points. We define the uniform mesh size $h_i={\displaystyle \frac{b_i-a_i}{M_i+1}}$ and generate a set of equidistant nodes $x_{i,j}=a_i+j{h}_i$ for $j=0,1,\dots ,M_i+1$.

Let $\mathbb{Z}$ denote the set of integers. For each direction i, we define the one-dimensional interior and global index sets as follows:

$$\begin{array}{cc}\hfill I_i& =\{j\in \mathbb{Z}\mid 1\le j\le M_i\},\hfill \\ \hfill {\widehat{I}}_i& =\{j\in \mathbb{Z}\mid 0\le j\le M_i+1\}.\hfill \end{array}$$

Subsequently, we construct the d-dimensional multi-index sets via the Cartesian product $K={\prod }_{i=1}^d{I}_i$ and $\widehat{K}={\prod }_{i=1}^d{\widehat{I}}_i$. A multi-index $J=(j_1,j_2,\dots ,j_d)\in \widehat{K}$ uniquely identifies a spatial grid point $X_J=(x_{1,j_1},x_{2,j_2},\dots ,x_{d,j_d})$. The collection of all such nodes forms the d-dimensional spatial mesh ${\Omega }_h=\{X_J\mid J\in \widehat{K}\}$.

2.2. Fully Discrete Scheme

In this subsection, we reformulate the multidimensional spatial fractional Cahn-Hilliard (FCH) equation based on the Scalar Auxiliary Variable (SAV) approach. First, we restate the model equations and rewrite the two-dimensional FCH model into the following equivalent coupled system:

$$\left\{\begin{array}{c}{\displaystyle \frac{\delta u(x,t)}{\delta t}}={\mathcal{L}}_x^{\alpha }v(x,t),\hfill \\ v(x,t)=-{\epsilon }^2{\mathcal{L}}_x^{\alpha }u(x,t)+{\displaystyle \frac{r\left(t\right)}{\sqrt{{\int }_{\Omega }F\left(u\right)dx+C_0}}}f\left(u(x,t)\right),\hfill \\ {\displaystyle \frac{dr\left(t\right)}{dt}}={\displaystyle \frac{1}{2\sqrt{{\int }_{\Omega }F\left(u\right)dx+C_0}}}{\int }_{\Omega }f\left(u\right)u_{t}dx.\hfill \end{array}\right.$$

(7)

2.2.1. Function Spaces and Continuous Operators

Consider the continuous Riesz fractional derivative operator ${\mathcal{L}}_{x_i}^{{\alpha }_i}$ ($1<{\alpha }_i\le 2$) in the $x_i$-direction defined on the bounded domain $\Omega$.Let ${\mathcal{L}}_x^{\alpha }={\sum }_{i=1}^d{\mathcal{L}}_{x_i}^{{\alpha }_i}$. Under homogeneous Dirichlet boundary conditions, the operator $-{\mathcal{L}}_{x_i}^{{\alpha }_i}$ is symmetric positive definite. Accordingly, there exists a continuous root operator denoted by ${\tilde{\mathcal{L}}}_{x_i}^{{\alpha }_i}$, which satisfies the following inner product identity:

$$(-{\mathcal{L}}_{x_i}^{{\alpha }_i}v,v)=({\tilde{\mathcal{L}}}_{x_i}^{{\alpha }_i}v,{\tilde{\mathcal{L}}}_{x_i}^{{\alpha }_i}v)={\parallel {\tilde{\mathcal{L}}}_{x_i}^{{\alpha }_i}v\parallel }^2.$$

(8)

Throughout this paper, $\parallel ·\parallel$ and $(·,·)$ represent the standard $L^2$ norm and inner product, respectively.

2.2.2. Multi-Dimensional Fully Discrete Scheme

Based on the d-dimensional spatial mesh, we proceed with the temporal discretization of the fractional Cahn-Hilliard equation. Let $\tau =T/N$ be the time step size and $t_k=k\tau$. We denote $U_J^k$ and $V_J^k$ as the numerical approximations of the phase-field variable $u(X,t)$ and the chemical potential $v(X,t)$ at the spatial node $X_J$ and time level $t_k$, respectively. Furthermore, we introduce the auxiliary scalar approximation $r^k\approx r\left(t_k\right)$.

To spatially discretize the continuous Riesz fractional derivative, we adopt the fractional centered difference formula. Let $U_J$ and $V_J$ be the numerical approximations of $u(X,t)$ and $v(X,t)$ at the spatial grid point $X_J$, respectively. The discrete fractional operator ${\delta }_{x_i}^{{\alpha }_i}$ with respect to the i-th spatial variable $x_i$ at point $X_J$ is expressed as:

$${\delta }_{x_i}^{{\alpha }_i}{U}_J:=-{\displaystyle \frac{1}{h_i^{{\alpha }_i}}}\sum _{L\in {\mathcal{N}}_i\left(J\right)}{g}_{\left|L\right(i)-J(i\left)\right|}^{\left({\alpha }_i\right)}{U}_L+O\left(h_i^6\right),J\in \mathbb{K},i=1,2,\dots ,d,$$

(9)

where ${\mathcal{N}}_i\left(J\right)$ denotes the set of multi-indices varying exclusively along the i-th dimension, and $L\left(i\right),J\left(i\right)$ represent their corresponding index components in that specific direction.

Let $M={\prod }_{i=1}^d{M}_i$ be the total number of interior spatial degrees of freedom. Specifically, for the one-dimensional case (where $M=M_1$), the discrete operator (9) simplifies to the following scalar summation:

$${\delta }_x^{\alpha }{U}_J=\sum _{i=1}^d{\delta }_{x_i}^{{\alpha }_i}{U}_J+O\left(h_i^6\right)(h_i\to 0)$$

(10)

where the centered difference weights $g_k^{\left({\alpha }_i\right)}$ are explicitly defined as follows:

$$g_k^{{\alpha }_i}={\displaystyle \frac{3}{2}}{a}_k^{{\alpha }_i}-{\mathbf{1}}_{\left\{2\right|n\}}·{\displaystyle \frac{3\left(2^{-{\alpha }_i}\right)}{5}}{a}_{k/2}^{{\alpha }_i}+{\mathbf{1}}_{\left\{3\right|n\}}·{\displaystyle \frac{3^{-{\alpha }_i}}{10}}{a}_{k/3}^{{\alpha }_i},$$

(11)

It is mathematically evident from (11) that for any fractional order $\alpha \in (1,2)$, the sequence $g_k^{\left({\alpha }_i\right)}$ does not truncate to zero for $k\ne 0$. Consequently, applying this numerical scheme yields a fully dense and symmetric Toeplitz coefficient matrix, which rigorously captures the non-local physical nature of the fractional phase-field model.

To handle the extreme stiffness of the nonlinear term explicitly while preserving unconditional energy stability, we define the second-order extrapolation ${\overline{U}}_J^{k+1}=2U_J^k-U_J^{k-1}$ and the corresponding nonlinear functional approximation:

$$W\left({\overline{U}}_J^{k+1}\right)={\displaystyle \frac{f\left({\overline{U}}_J^{k+1}\right)}{\sqrt{\left({\prod }_{i=1}^d{h}_i\right){\sum }_{J\in K}F\left({\overline{U}}_J^{k+1}\right)+C_0}}}.$$

(12)

Applying the SAV approach combined with the second-order backward differentiation formula (BDF2), we obtain the following d-dimensional fully discrete scheme:

$$\left\{\begin{array}{c}{\displaystyle \frac{3U_J^{k+1}-4U_J^k+U_J^{k-1}}{2\tau }=\sum _{i=1}^d{\delta }_{x_i}^{{\alpha }_i}{V}_J^{k+1},}\hfill \\ {\displaystyle V_J^{k+1}=-{\epsilon }^2\sum _{i=1}^d{\delta }_{x_i}^{{\alpha }_i}{U}_J^{k+1}+r^{k+1}W\left({\overline{U}}_J^{k+1}\right),}\hfill \\ {\displaystyle \frac{3r^{k+1}-4r^k+r^{k-1}}{2\tau }=\frac{1}{2}\left(\prod _{i=1}^d{h}_i\right)\sum _{J\in K}W\left({\overline{U}}_J^{k+1}\right)\frac{3U_J^{k+1}-4U_J^k+U_J^{k-1}}{2\tau },}\hfill \end{array}\right.$$

(13)

which is subject to the following homogeneous boundary and initial conditions:

$$\left\{\begin{array}{cc}{U}_J^k=V_J^k=0,\hfill & J\in \delta K,0\le k\le N,\hfill \\ U_J^0=U_0\left(X_J\right),U_J^1=U_1\left(X_J\right),\hfill & J\in K,\hfill \\ r^0=\sqrt{\left({\prod }_{i=1}^d{h}_i\right){\sum }_{J\in K}F\left(U_J^0\right)+C_0}\hfill \\ r^1=\sqrt{\left({\prod }_{i=1}^d{h}_i\right){\sum }_{J\in K}F\left(U_J^1\right)+C_0}.\hfill \end{array}\right.$$

(14)

2.3. Fast Preconditioner and Its Tensor Batching Optimization

In this subsection, we introduce the spectral-equivalent $\tau$-preconditioner to decouple the severely ill-conditioned linear system derived from the SAV-BDF2 scheme. To address the curse of dimensionality encountered by this preconditioner in high-dimensional spatial problems, we further propose a Vectorized Tensor Processing (VTP) architecture for low-level computational acceleration.

2.3.1. Preconditioned Symmetrized System Under SAV-BDF2

Following the spatial discretization of the multidimensional Riesz fractional operator, the intrinsic mechanism of the SAV-BDF2 scheme yields an asymmetric coupled algebraic system. To legitimately invoke Krylov subspace solvers such as MINRES, this structural asymmetry must be algebraically resolved.

By introducing the specific scaling parameters $c_1={\epsilon }^2$ and $c_2={\displaystyle \frac{2\tau }{3}}$, alongside the auxiliary multiplier $\theta =\sqrt{c_1/c_2}$ and the equivalent aggregated coefficient $c=\sqrt{c_1{c}_2}$, we strictly recast the original coupled partial differential system into a highly symmetric $2\times 2$ block linear topology:

$$\mathcal{A}\left[\begin{array}{c}{\tilde{\mathbf{U}}}^{k+1}\\ {\mathbf{V}}^{k+1}\end{array}\right]=\left[\begin{array}{cc}-cG& I\\ I& cG\end{array}\right]\left[\begin{array}{c}\theta {\mathbf{U}}^{k+1}\\ {\mathbf{V}}^{k+1}\end{array}\right]=\left[\begin{array}{c}{\mathbf{f}}_2\\ \theta {\displaystyle \frac{2\tau }{3}}{\mathbf{f}}_1\end{array}\right].$$

(15)

As the spatial grid step $h\to 0$ and the physical interface width $\epsilon \to 0$, the dense coefficient matrix G inevitably drives the global block matrix $\mathcal{A}$ into severe ill-conditioning. To fundamentally circumvent this computational hurdle in multidimensional scenarios, we adopt the robust block diagonal preconditioner $\mathcal{P}$, which is inspired by the sine-transform-based technique proposed in [6]:

$$\mathcal{P}=\left[\begin{array}{cc}c\tau \left(G\right)+I& 0\\ 0& c\tau \left(G\right)+I\end{array}\right],$$

(16)

where $\tau \left(G\right)={\sum }_{i=1}^d\left(I_{M_i^{-}}\otimes \tau \left(G_{{\alpha }_i}\right)\otimes I_{M_i^{+}}\right)$ denotes the d-dimensional global approximation matrix.

For clarity and practical implementation, based on our spatial domain discretization, each sub-block $\tau \left(G_{{\alpha }_i}\right)$ can be rigorously diagonalized by the 1D discrete sine transform (DST) matrix $S_{M_i}$, satisfying $\tau \left(G_{{\alpha }_i}\right)=S_{M_i}{\Lambda }_{{\alpha }_i}{S}_{M_i}$. The internal elements of the orthogonal transformation matrix $S_{M_i}$ and the diagonal eigenvalue matrix ${\Lambda }_{{\alpha }_i}$ are explicitly defined as follows [6]:

$$\begin{array}{cc}\hfill {\left[S_{M_i}\right]}_{j,l}& =\sqrt{{\displaystyle \frac{2}{M_i+1}}}sin\left({\displaystyle \frac{jl\pi }{M_i+1}}\right),\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\Lambda }_{{\alpha }_i}{]}_j& =g_0^{\left({\alpha }_i\right)}+2\sum _{k=1}^{M_i-1}{g}_k^{\left({\alpha }_i\right)}cos\left({\displaystyle \frac{kj\pi }{M_i+1}}\right),\hfill \end{array}$$

(18)

for $j,l\in \{1,2,\dots ,M_i\}$, where $g_k^{\left({\alpha }_i\right)}$ represent the generation coefficients derived from the fractional difference scheme. This flawless tensor diagonalization mechanism guarantees that the inverse action ${\mathcal{P}}^{-1}$ can be executed with quasi-optimal complexity via the fast Fourier transform.

Remark. 

It has been strictly proven in [5] that the eigenvalues of the preconditioned matrix ${\mathcal{P}}^{-1}\mathcal{A}$ are uniformly bounded within the disjoint intervals $[-{\displaystyle \frac{3}{2}},-{\displaystyle \frac{1}{2\sqrt{2}}}]\cup [{\displaystyle \frac{1}{2\sqrt{2}}},{\displaystyle \frac{3}{2}}]$. From an algebraic perspective, this spectral distribution is intrinsically immune to the specific value of the temporal scaling factor $c>0$. Consequently, the introduction of our new coefficient $c_2={\displaystyle \frac{2\tau }{3}}$, which is exclusively dictated by the BDF2 discretization, preserves these theoretical bounds intact. This mathematical property ensures that the block diagonal preconditioner remains highly robust, ultimately yielding a mesh-independent linear convergence rate for the asymmetric SAV-BDF2 scheme.

2.3.2. Vectorized Tensor Processing (VTP) Architecture

From a purely algebraic perspective, the theoretical time complexity of applying the d-dimensional preconditioner ${\mathcal{P}}^{-1}$ is mathematically bounded by $\mathcal{O}(N{\sum }_{i=1}^{d}log{M}_i)$.

Nevertheless, under standard routines, executing a d-dimensional joint sine transform usually means slicing the global tensor into isolated one-dimensional vectors. The system then launches exactly $N/M_k$ independent scalar DST calls along the k-th dimension. This mechanism tears memory continuity apart. Such granular data extraction inevitably causes steep instruction overhead and massive cache misses. Hampered by these limits, computational efficiency in high-dimensional scenarios simply collapses.

To fundamentally dismantle this hardware-level barrier, we propose the Vectorized Tensor Processing (VTP) architecture. Let $\mathcal{U}\in {\mathbb{R}}^{M_1\times M_2\times \dots \times M_d}$ denote the d-dimensional physical state tensor isomorphic to the discrete algebraic system. When the sine transform is required along the k-th dimension, an orthogonal permutation operator ${\mathcal{T}}_k$ is first introduced to advance the target computational dimension to the leading position in memory. Subsequently, this permuted tensor is compulsorily reshaped into a highly dense two-dimensional matrix manifold:

$${\widehat{\mathcal{U}}}^{\left(k\right)}=\mathrm{Reshape}\left({\mathcal{T}}_k\left(\mathcal{U}\right),\left[M_k,{\displaystyle \frac{N}{M_k}}\right]\right).$$

(19)

Governed by this absolutely contiguous memory mapping mechanism, the initially fragmented $N/M_k$ scalar transforms are instantaneously aggregated into a single, high-density batched execution block. To strictly avoid the complexity regression that would be caused by explicit dense matrix multiplication, the transformation is executed holistically via a Batched Fast Discrete Sine Transform (Batched FDST) operator ${\mathcal{F}}_{\mathrm{DST}}^{\left(k\right)}$ over the continuous column space of the manifold:

$${\mathbf{Z}}^{\left(k\right)}={\mathcal{F}}_{\mathrm{DST}}^{\left(k\right)}\left({\widehat{\mathcal{U}}}^{\left(k\right)}\right).$$

(20)

Upon the completion of this batched fast operation, the inverse permutation operator ${\mathcal{T}}_k^{-1}$ is immediately applied, accurately mapping the data back to its original physical coordinate system.

By reconstructing the underlying dataflow, the VTP architecture elegantly collapses the staggering $\mathcal{O}(N/M_k)$ fragmented function calls within a single iteration step into strictly $\mathcal{O}\left(d\right)$ macroscopic batched operations. Pinning the spatial footprint at $\mathcal{O}\left(N\right)$ while maintaining the optimal algebraic time complexity, this tensorization strategy decisively breaks through the execution bottlenecks inherent in high-dimensional discrete schemes right down at the bare-metal hardware level.

While strictly maintaining a minimal spatial complexity of $\mathcal{O}\left(N\right)$ and the optimal algebraic time complexity, this tensorization strategy decisively breaks through the execution bottlenecks inherent in high-dimensional discrete schemes at the bare-metal hardware level.

3. Results

3.1. Energy Stability Analysis

In this subsection, we initially formulate the necessary definitions ensuring the regularity of the d-dimensional discrete operations, and subsequently provide a rigorous proof establishing the unconditional energy stability of the proposed scheme.

3.1.1. Necessary Definitions

Let U be the space of grid functions defined on all interior mesh points, arranged in lexicographical order. For any grid functions $u,v\in U$, we define the discrete inner product and the corresponding $L^2$ norm in the d-dimensional space as:

$${(u,v)}_h=\left(\prod _{i=1}^d{h}_i\right)\sum _{J\in K}{u}_J{v}_J,{\parallel v\parallel }_h=\sqrt{(v,v)}.$$

(21)

For each spatial direction $i\in \{1,\dots ,d\}$, the discrete fractional operator $-{\delta }_{x_i}^{{\alpha }_i}$ is symmetric positive definite on the space U. Correspondingly, there exists a discrete root operator, denoted by ${\tilde{\delta }}_{x_i}^{{\alpha }_i}$, such that the inner product identity $(-{\delta }_{x_i}^{{\alpha }_i}u,u)=({\tilde{\delta }}_{x_i}^{{\alpha }_i}u,{\tilde{\delta }}_{x_i}^{{\alpha }_i}u)={\parallel {\tilde{\delta }}_{x_i}^{{\alpha }_i}u\parallel }_h^2$ holds natively under homogeneous boundary conditions.

3.1.2. Proof of Energy Decay

To establish the energy dissipation property of the fully discrete SAV-BDF2 scheme in high-dimensional settings, we first introduce the crucial algebraic identity for the BDF2 formulation.

Lemma 3.1 For any sequence of vectors $\left\{a^k\right\}$ in a real inner product space, the following identity holds unconditionally:

$$\begin{array}{cc}\hfill 2(a^{k+1},3a^{k+1}-4a^k+a^{k-1})=& \parallel a^{k+1}{\parallel }_h^2-\parallel a^k{\parallel }_h^2+\parallel 2a^{k+1}-a^k{\parallel }_h^2-{\parallel 2a^k-a^{k-1}\parallel }_h^2\hfill \\ \hfill & +\parallel a^{k+1}-2a^k+a^{k-1}{\parallel }_h^2.\hfill \end{array}$$

(22)

Based on the SAV-BDF2 time discretization, we define the modified discrete total energy ${\tilde{E}}^{k+1}$ at the $(k+1)$-th time level as:

$${\tilde{E}}^{k+1}={\displaystyle \frac{{\epsilon }^2}{2}}\sum _{i=1}^d\left(\parallel {\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1}{\parallel }_h^2+{\parallel 2{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1}-{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^k\parallel }_h^2\right)+\left({\left(r^{k+1}\right)}^2+{(2r^{k+1}-r^k)}^2\right).$$

(23)

Theorem 3.1: The d-dimensional fully discrete SAV-BDF2 scheme guarantees that the modified total energy decreases monotonically, meaning ${\tilde{E}}^{k+1}\le {\tilde{E}}^k$ holds unconditionally for any time step $\tau >0$.

Proof: Taking the discrete inner product of the first equation in the scheme with $V^{k+1}$, and utilizing the discrete fractional summation-by-parts equivalent, we obtain:

$$\left({\displaystyle \frac{3U^{k+1}-4U^k+U^{k-1}}{2\tau }},V^{k+1}\right)=-\sum _{i=1}^d{∥{\tilde{\delta }}_{x_i}^{{\alpha }_i}{V}^{k+1}∥}_h^2.$$

(24)

Taking the discrete inner product of the second equation with the time derivative term ${\displaystyle \frac{3U^{k+1}-4U^k+U^{k-1}}{2\tau }}$ yields:

$$\begin{array}{cc}\hfill \left(V^{k+1},{\displaystyle \frac{3U^{k+1}-4U^k+U^{k-1}}{2\tau }}\right)& ={\epsilon }^2\sum _{i=1}^d\left({\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1},{\tilde{\delta }}_{x_i}^{{\alpha }_i}{\displaystyle \frac{3U^{k+1}-4U^k+U^{k-1}}{2\tau }}\right)\hfill \\ \hfill & +r^{k+1}\left(W\left({\overline{U}}^{k+1}\right),{\displaystyle \frac{3U^{k+1}-4U^k+U^{k-1}}{2\tau }}\right).\hfill \end{array}$$

(25)

Equating (24) and (25), and multiplying both sides by $2\tau$, we deduce:

$$\begin{array}{cc}\hfill -2\tau \sum _{i=1}^d{∥{\tilde{\delta }}_{x_i}^{{\alpha }_i}{V}^{k+1}∥}^2& ={\epsilon }^2\sum _{i=1}^d\left({\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1},3{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1}-4{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^k+{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k-1}\right)\hfill \\ \hfill & +r^{k+1}\left(W\left({\overline{U}}^{k+1}\right),3U^{k+1}-4U^k+U^{k-1}\right).\hfill \end{array}$$

(26)

From the third equation of the SAV-BDF2 scheme, the nonlinear inner product term can be exactly substituted by the auxiliary scalar variables:

$$\left(W\left({\overline{U}}^{k+1}\right),3U^{k+1}-4U^k+U^{k-1}\right)=2(3r^{k+1}-4r^k+r^{k-1}).$$

(27)

Substituting (27) into (26), the equation transforms into a strictly quadratic form:

$$-2\tau \sum _{i=1}^d{∥{\tilde{\delta }}_{x_i}^{{\alpha }_i}{V}^{k+1}∥}_h^2={\epsilon }^2\sum _{i=1}^d\left({\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1},3{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1}-4{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^k+{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k-1}\right)+2r^{k+1}(3r^{k+1}-4r^k+r^{k-1}).$$

(28)

Using the BDF2 polarization identity (22) to both the fractional gradient term and the scalar term, (28) rigorously expands to:

$$\begin{array}{cc}\hfill -2\tau \sum _{i=1}^d{∥{\tilde{\delta }}_{x_i}^{{\alpha }_i}{V}^{k+1}∥}_h^2& ={\displaystyle \frac{{\epsilon }^2}{2}}\sum _{i=1}^d[\parallel {\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1}{\parallel }_h^2-\parallel {\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^k{\parallel }_h^2+{\parallel 2{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1}-{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^k\parallel }_h^2\hfill \\ \hfill & -\parallel 2{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^k-{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k-1}{\parallel }_h^2+{\parallel {\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1}-2{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^k+{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k-1}\parallel }_h^2]\hfill \\ \hfill & +[{\left(r^{k+1}\right)}^2-{\left(r^k\right)}^2+{(2r^{k+1}-r^k)}^2\hfill \\ \hfill & -{(2r^k-r^{k-1})}^2+{(r^{k+1}-2r^k+r^{k-1})}^2].\hfill \end{array}$$

(29)

Rearranging the terms according to the time levels $k+1$ and k, and utilizing the definition of the modified discrete total energy (23), we arrive at the discrete dissipation equation:

$${\tilde{E}}^{k+1}-{\tilde{E}}^k+{\displaystyle \frac{{\epsilon }^2}{2}}\sum _{i=1}^d{\parallel {\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k+1}-2{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^k+{\tilde{\delta }}_{x_i}^{{\alpha }_i}{U}^{k-1}\parallel }_h^2+{(r^{k+1}-2r^k+r^{k-1})}^2=-2\tau \sum _{i=1}^d{∥{\tilde{\delta }}_{x_i}^{{\alpha }_i}{V}^{k+1}∥}_h^2.$$

(30)

Since the second-order numerical dissipation terms on the left side of (30) are strictly non-negative, and the right side satisfies $-2\tau {\sum }_{i=1}^d{\parallel {\tilde{\delta }}_{x_i}^{{\alpha }_i}{V}^{k+1}\parallel }^2\le 0$, we logically deduce ${\tilde{E}}^{k+1}\le {\tilde{E}}^k$. This rigorously completes the proof.

3.2. Global Convergence Analysis

Rigorous global convergence proofs for the fully discrete SAV-BDF2 scheme unfold throughout this subsection. By confining the analysis to a zero-mean subspace, the mathematical setup natively permits direct energy estimation inside the discrete negative fractional Sobolev space, essentially laying the prerequisite groundwork for all subsequent derivations. Strategically, computational pairing tricks like the polarization identity extract inner products into measurable norms for comparative scaling. To guarantee an exact algebraic cancellation of those notoriously stiff nonlinear components, the cross-terms remain deliberately untouched as inner products rather than suffering premature expansion. These residual truncation errors get cordoned off into independent norms prior to expansion. A closed-loop error inequality system immediately follows. Working through explicit coefficient matching and reliable scaling bounds ultimately secures the theoretical justification for the numerical convergence.

3.2.1. Discrete Matrix, Equivalent Norms, and Core Lemmas

Definition 3  

(Discrete Fractional Matrix and Equivalent Norms). Let ${\Omega }_h$ denote the d-dimensional uniform spatial mesh. Let $\mathcal{L}={\sum }_{i=1}^d{\delta }_{x_i}^{{\alpha }_i}$ be the original symmetric negative definite fractional operator corresponding to the continuous Laplacian $\Delta$. To establish a valid discrete Sobolev topology, we define the corresponding symmetric positive definite (SPD) matrix as $\mathcal{A}=-\mathcal{L}=-{\sum }_{i=1}^d{\delta }_{x_i}^{{\alpha }_i}$.

Within the discrete zero-mean subspace ${\mathring{U}}_h=\{U\in {\mathbb{R}}^M\mid {\sum }_{J\in K}{U}_J=0\}$, all eigenvalues of $\mathcal{A}$ are strictly positive, ensuring the unique existence of its inverse matrix ${\mathcal{A}}^{-1}$. In addition to the standard discrete $L^2$-norm ${\parallel U\parallel }_h=\sqrt{{\langle U,U\rangle }_h}$, the positive definiteness induces the equivalent fractional norms:

$${\parallel U\parallel }_H=\sqrt{{\langle U,\mathcal{A}U\rangle }_h},{\parallel U\parallel }_{H^{-1}}=\sqrt{{\langle U,{\mathcal{A}}^{-1}U\rangle }_h}.$$

(31)

Given that ${\mathring{U}}_h$ constitutes a real Hilbert space and $\mathcal{A}$ is self-adjoint, the generalized Cauchy-Schwarz inequality ${\langle U,V\rangle }_h\le {\parallel U\parallel }_H{\parallel V\parallel }_{H^{-1}}$ holds unconditionally.

Lemma 4  

(Mesh-Independent Discrete Interpolation and Poincaré Inequalities). For any given penalty parameter $\gamma >0$ independent of the temporal step size $\tau$ and the spatial step size h, the following discrete interpolation inequality unconditionally holds in the subspace ${\mathring{U}}_h$:

$${\parallel U\parallel }_h^2\le {\displaystyle \frac{\gamma }{2}}{\parallel U\parallel }_H^2+{\displaystyle \frac{1}{2\gamma }}{\parallel U\parallel }_{H^{-1}}^2,\forall U\in {\mathring{U}}_h.$$

(32)

Furthermore, let ${\lambda }_1=\mathcal{O}\left(1\right)>0$ be the minimum non-zero eigenvalue of $\mathcal{A}$ on ${\mathring{U}}_h$. There exists a global constant $C_1={\lambda }_1^{-1/2}$, strictly independent of the mesh size h, such that ${\parallel U\parallel }_{H^{-1}}\le C_1{\parallel U\parallel }_h$ and ${\parallel U\parallel }_h\le C_1{\parallel U\parallel }_H$. This implies the discrete Poincaré inverse control inequality:

$${\parallel U\parallel }_{H^{-1}}\le C_1^2{\parallel U\parallel }_H,\forall U\in {\mathring{U}}_h.$$

(33)

Lemma 5  

(Discrete Lipschitz Continuity of the SAV Functional). The modified nonlinear functional takes the composite form $W\left(u\right)={\displaystyle \frac{f\left(u\right)}{\sqrt{E\left(u\right)+C_0}}}$, where $f\left(u\right)=u^3-u$, the truncation constant is $C_0>0$, and $E\left(u\right)={\langle F\left(u\right),1\rangle }_h$ acts as a discrete spatial scalar.

Assume that up to the current time step $t_m$, the numerical solutions are confined within a compact set ${\parallel u\parallel }_{L^{\infty }},{\parallel U^k\parallel }_{L^{\infty }}\le C_2$. Under this constraint, $W(·)$ satisfies a global Lipschitz continuity in the discrete space. Introducing the second-order temporal extrapolation $u\left(t_{k+1}\right)-{\overline{U}}^{k+1}={\xi }_{\tau }^{k+1}+2E_u^k-E_u^{k-1}$, where $\parallel {\xi }_{\tau }^{k+1}{\parallel }_H\le C{\tau }^2$ by Taylor’s theorem. By the discrete Poincaré inequality, we have $\parallel {\xi }_{\tau }^{k+1}{\parallel }_h\le C_1{\parallel {\xi }_{\tau }^{k+1}\parallel }_H\le \mathcal{O}\left({\tau }^2\right)$. Since the denominator is an integral scalar independent of spatial coordinates, the fractional matrix directly penetrates it: ${\mathcal{A}}^{1/2}W\left(u\right)=\left({\mathcal{A}}^{1/2}f\left(u\right)\right)/\sqrt{E\left(u\right)+C_0}$. Leveraging the discrete fractional product rule, there exist independent constants $C_3,C_4>0$ such that the nonlinear bounds are governed by:

$$\begin{array}{cc}\hfill \parallel W\left(u\left(t_{k+1}\right)\right)-W\left({\overline{U}}^{k+1}\right){\parallel }_h^2& \le C_3\left(\parallel {\xi }_{\tau }^{k+1}{\parallel }_h^2+\parallel E_u^k{\parallel }_h^2+{\parallel E_u^{k-1}\parallel }_h^2\right),\hfill \\ \hfill \parallel W\left(u\left(t_{k+1}\right)\right)-W\left({\overline{U}}^{k+1}\right){\parallel }_H^2& \le C_4\left(\parallel {\xi }_{\tau }^{k+1}{\parallel }_H^2+\parallel E_u^k{\parallel }_H^2+{\parallel E_u^{k-1}\parallel }_H^2\right).\hfill \end{array}$$

(34)

3.2.2. Discrete Error Evolution System

By evaluating the exact continuous solutions at the grid points $(X_J,t_{k+1})$ and substituting them into the fully discrete SAV-BDF2 scheme, we define the numerical errors as $E_u^{k+1}=u\left(t_{k+1}\right)-U^{k+1}$, $E_v^{k+1}=v\left(t_{k+1}\right)-V^{k+1}$, and $e_r^{k+1}=r\left(t_{k+1}\right)-R^{k+1}$. Applying the backward difference operator ${\delta }_t{\varphi }^{k+1}={\displaystyle \frac{3{\varphi }^{k+1}-4{\varphi }^k+{\varphi }^{k-1}}{2\tau }}$, the consistency equations yield the coupled error evolution system:

$$\left\{\begin{array}{cc}\hfill {\delta }_t{E}_u^{k+1}& =-\mathcal{A}{E}_v^{k+1}+{\rho }_u^{k+1},\hfill \\ \hfill E_v^{k+1}& ={\epsilon }^2\mathcal{A}{E}_u^{k+1}+Q_N^{k+1}+{\rho }_v^{k+1},\hfill \\ \hfill {\delta }_t{e}_r^{k+1}& =Q_r^{k+1}+{\rho }_r^{k+1}.\hfill \end{array}\right.$$

(35)

The nonlinear phase-field residual and the auxiliary scalar residual are rigorously defined as:

$$\begin{array}{cc}\hfill Q_N^{k+1}& =r\left(t_{k+1}\right)W\left(u\left(t_{k+1}\right)\right)-R^{k+1}W\left({\overline{U}}^{k+1}\right),\hfill \\ \hfill Q_r^{k+1}& ={\displaystyle \frac{1}{2}}\left[{\langle W\left(u\left(t_{k+1}\right)\right),{\delta }_{t}u\left(t_{k+1}\right)\rangle }_h-{\langle W\left({\overline{U}}^{k+1}\right),{\delta }_t{U}^{k+1}\rangle }_h\right].\hfill \end{array}$$

(36)

Assumption 6  

(Truncation Errors and Initialization). According to Taylor expansion and difference approximation theories, the local truncation errors maintain consistent spatiotemporal accuracy (second-order in time, sixth-order in space), bounded by a global constant $C_5$:

$$\parallel {\rho }_u^{k+1}{\parallel }_{H^{-1}}^2+\parallel {\rho }_v^{k+1}{\parallel }_H^2+{\left|{\rho }_r^{k+1}\right|}^2\le C_5({\tau }^4+h^{12}).$$

(37)

As BDF2 is a two-step method, the initial physical conditions are prescribed exactly, implying $E_u^0=0$ and $e_r^0=0$. The numerical solutions at the first time step, $U^1$ and $R^1$, are assumed to be computed by a sufficiently accurate single-step method. There exists a constant $C_{10}$ such that the initialization errors strictly conform to the global convergence rate:

$$\parallel E_u^1{\parallel }_H\le C_{10}({\tau }^2+h^6),\left|e_r^1\right|\le C_{10}({\tau }^2+h^6).$$

(38)

3.2.3. Proof of the Global Convergence

Provided that Assumption 6 holds, the multidimensional fully discrete SAV-BDF2 scheme converges unconditionally. There exist global constants $C_{12}$ and $C_{13}$, independent of $\tau$ and h, such that for any valid time step m:

$$\parallel u\left(t_m\right)-U^m{\parallel }_{H^{-1}}\le C_{13}exp\left(C_{12}T\right)\left({\tau }^2+h^6\right).$$

(39)

Proof. 

Substituting the chemical potential equation into the phase-field evolution equation, we recast the system into a weak variational form acting on the inverse matrix ${\mathcal{A}}^{-1}$:

$${\mathcal{A}}^{-1}{\delta }_t{E}_u^{k+1}=-{\epsilon }^2\mathcal{A}{E}_u^{k+1}-Q_N^{k+1}-{\rho }_v^{k+1}+{\mathcal{A}}^{-1}{\rho }_u^{k+1}.$$

(40)

Taking the discrete $L^2$ inner product of both sides with the temporal increment test function $4\tau {\delta }_t{E}_u^{k+1}$ within ${\mathring{U}}_h$, we obtain:

$$\begin{array}{cc}\hfill 4\tau {\langle {\mathcal{A}}^{-1}{\delta }_t{E}_u^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_h& =-4\tau {\epsilon }^2{\langle \mathcal{A}{E}_u^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_h\hfill \\ \hfill & -4\tau {\langle Q_N^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_h\hfill \\ \hfill & -4\tau {\langle {\rho }_v^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_h\hfill \\ \hfill & +4\tau {\langle {\mathcal{A}}^{-1}{\rho }_u^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_h.\hfill \end{array}$$

(41)

The left-hand side yields the positive definite term $4\tau \parallel {\delta }_t{E}_u^{k+1}{\parallel }_{H^{-1}}^2$. For the first term on the right-hand side, the symmetric property translates it to $-4\tau {\epsilon }^2{\langle E_u^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_H$. Applying the algebraic BDF2 polarization identity $2{\langle 3A-4B+C,A\rangle }_H={\parallel A\parallel }_H^2-{\parallel B\parallel }_H^2{+\parallel 2A-B\parallel }_H^2{-\parallel 2B-C\parallel }_H^2+{\parallel A-2B+C\parallel }_H^2$, and dropping the non-negative second-order dissipation term $-{\epsilon }^2{\parallel E_u^{k+1}-2E_u^k+E_u^{k-1}\parallel }_H^2\le 0$, we define the baseline phase-field functional ${\mathbb{E}}_u^k={\epsilon }^2(\parallel E_u^k{\parallel }_H^2+\parallel 2E_u^k-E_u^{k-1}{\parallel }_H^2)$. Rearranging the terms yields:

$$\begin{array}{cc}\hfill {\mathbb{E}}_u^{k+1}-{\mathbb{E}}_u^k+4\tau {\parallel {\delta }_t{E}_u^{k+1}\parallel }_{H^{-1}}^2& \le -4\tau {\langle Q_N^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_h\hfill \\ \hfill & -4\tau {\langle {\rho }_v^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_h\hfill \\ \hfill & +4\tau {\langle {\rho }_u^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_{H^{-1}}.\hfill \end{array}$$

(42)

Concurrently, multiplying the auxiliary scalar evolution equation by $8\tau e_r^{k+1}$ ensures precise alignment with the $1/2$ factor in $Q_r^{k+1}$. Following a parallel BDF2 polarization and defining the scalar functional ${\mathbb{E}}_r^k=2|e_r^k{|}^2+2{|2e_r^k-e_r^{k-1}|}^2$, dropping the positive residual gives:

$${\mathbb{E}}_r^{k+1}-{\mathbb{E}}_r^k\le 8\tau e_r^{k+1}{Q}_r^{k+1}+8\tau e_r^{k+1}{\rho }_r^{k+1}.$$

(43)

To handle the nonlinear matrix cancellation, we define the total energy functional ${\mathcal{H}}^k={\mathbb{E}}_u^k+{\mathbb{E}}_r^k$. Summing (42) and (43) isolates the implicit nonlinear array $I_{NL}=-4\tau {\langle Q_N^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_h+8\tau e_r^{k+1}{Q}_r^{k+1}$. Substituting the definition of $Q_N^{k+1}$ and applying ${\delta }_t{U}^{k+1}={\delta }_{t}u\left(t_{k+1}\right)-{\delta }_t{E}_u^{k+1}$ within $Q_r^{k+1}$, the array expands as:

$$\begin{array}{cc}\hfill I_{NL}& =-4\tau e_r^{k+1}{\langle W\left(u\left(t_{k+1}\right)\right),{\delta }_t{E}_u^{k+1}\rangle }_h\hfill \\ \hfill & -4\tau R^{k+1}{\langle W\left(u\left(t_{k+1}\right)\right)-W\left({\overline{U}}^{k+1}\right),{\delta }_t{E}_u^{k+1}\rangle }_h\hfill \\ \hfill & +4\tau e_r^{k+1}{\langle W\left(u\left(t_{k+1}\right)\right)-W\left({\overline{U}}^{k+1}\right),{\delta }_{t}u\left(t_{k+1}\right)\rangle }_h\hfill \\ \hfill & +4\tau e_r^{k+1}{\langle W\left({\overline{U}}^{k+1}\right),{\delta }_t{E}_u^{k+1}\rangle }_h.\hfill \end{array}$$

(44)

The first and fourth terms exactly cancel. Extracting the common inner product factor with the $R^{k+1}$ term and invoking the mapping $e_r^{k+1}+R^{k+1}\equiv r\left(t_{k+1}\right)$, the implicit nonlinear terms algebraically collapse to:

$$\begin{array}{cc}\hfill I_{NL}& =-4\tau r\left(t_{k+1}\right){\langle W\left(u\left(t_{k+1}\right)\right)-W\left({\overline{U}}^{k+1}\right),{\delta }_t{E}_u^{k+1}\rangle }_h\hfill \\ \hfill & +4\tau e_r^{k+1}{\langle W\left(u\left(t_{k+1}\right)\right)-W\left({\overline{U}}^{k+1}\right),{\delta }_{t}u\left(t_{k+1}\right)\rangle }_h.\hfill \end{array}$$

(45)

To bound the remaining terms, we apply the weighted Young’s inequality ($ab\le a^2+b^2/4$). Utilizing the Cauchy-Schwarz inequality, the three terms coupled with ${\delta }_t{E}_u^{k+1}$ are decomposed as follows:

$$\begin{array}{cc}\hfill 4\tau \left|r\left(t_{k+1}\right){\langle W\left(u\right)-W\left(\overline{U}\right),{\delta }_t{E}_u^{k+1}\rangle }_h\right|& \le 4\tau |r\left(t_{k+1}\right){|}^2\parallel W\left(u\right)-W\left(\overline{U}\right){\parallel }_H^2+\tau {\parallel {\delta }_t{E}_u^{k+1}\parallel }_{H^{-1}}^2,\hfill \\ \hfill 4\tau \left|{\langle {\rho }_v^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_h\right|& \le 4\tau \parallel {\rho }_v^{k+1}{\parallel }_H^2+\tau {\parallel {\delta }_t{E}_u^{k+1}\parallel }_{H^{-1}}^2,\hfill \\ \hfill 4\tau \left|{\langle {\rho }_u^{k+1},{\delta }_t{E}_u^{k+1}\rangle }_{H^{-1}}\right|& \le 4\tau \parallel {\rho }_u^{k+1}{\parallel }_{H^{-1}}^2+\tau {\parallel {\delta }_t{E}_u^{k+1}\parallel }_{H^{-1}}^2.\hfill \end{array}$$

(46)

The summation of the extracted norms equals precisely $3\tau \parallel {\delta }_t{E}_u^{k+1}{\parallel }_{H^{-1}}^2$. This quantity is fully absorbed by the positive dissipation term $4\tau \parallel {\delta }_t{E}_u^{k+1}{\parallel }_{H^{-1}}^2$ on the left-hand side, leaving a non-negative remainder $\tau \parallel {\delta }_t{E}_u^{k+1}{\parallel }_{H^{-1}}^2$. For the terms independent of ${\delta }_t{E}_u^{k+1}$, subject to the physical bounds $\parallel {\delta }_t{u\parallel }_{L^{\infty }}\le C_6$ and $\left|r\left(t\right)\right|\le C_7$:

$$\begin{array}{cc}\hfill 4\tau e_r^{k+1}{\langle W\left(u\right)-W\left(\overline{U}\right),{\delta }_{t}u\rangle }_h& \le \tau |e_r^{k+1}{|}^2+4\tau C_6^2{\parallel W\left(u\right)-W\left(\overline{U}\right)\parallel }_h^2,\hfill \\ \hfill 8\tau e_r^{k+1}{\rho }_r^{k+1}& \le \tau |e_r^{k+1}{|}^2+16\tau {\left|{\rho }_r^{k+1}\right|}^2.\hfill \end{array}$$

(47)

Substituting the global Lipschitz bounds from Lemma 5, let $C_8$ denote the maximum combined constant arising from the Lipschitz bounds and the Young inequality coefficients. Aggregating these terms exposes the phase-field error history for both the k-th and $(k-1)$-th steps in strong and weak topologies: $\tau C_8(\parallel E_u^k{\parallel }_H^2+\parallel E_u^{k-1}{\parallel }_H^2)+\tau C_8(\parallel E_u^k{\parallel }_h^2+\parallel E_u^{k-1}{\parallel }_h^2)$. To decouple the weak norm’s dependence on the mesh layout, we invoke Lemma 4 strictly on the $L^2$ components, assigning a penalty parameter $\gamma =1/{\epsilon }^2$:

$$\tau C_8{\parallel E_u^j\parallel }_h^2\le \tau C_8\left({\displaystyle \frac{1}{2{\epsilon }^2}}\parallel E_u^j{\parallel }_H^2+{\displaystyle \frac{{\epsilon }^2}{2}}{\parallel E_u^j\parallel }_{H^{-1}}^2\right),\mathrm{for}j\in \{k,k-1\}.$$

(48)

We cascade the resulting dual norms into the discrete Poincaré inequality $\parallel E_u^j{\parallel }_{H^{-1}}^2\le C_1^4{\parallel E_u^j\parallel }_H^2$ to elevate their topologies to the H-norm space. Combining these elevated terms with the pre-existing high-order terms, the total phase-field historical bound is consolidated to

$$\tau C_8\left[1+{\displaystyle \frac{1}{2{\epsilon }^2}}+{\displaystyle \frac{C_1^4{\epsilon }^2}{2}}\right](\parallel E_u^k{\parallel }_H^2+\parallel E_u^{k-1}{\parallel }_H^2).$$

(49)

Concurrently, the scalar residual previously bounded by $2\tau |e_r^{k+1}{|}^2$ natively satisfies $2\tau |e_r^{k+1}{|}^2\le \tau {\mathcal{H}}^{k+1}$. Factoring out the physical parameter ${\epsilon }^2$ to match the energy functional structure (${\mathcal{H}}^j\ge {\epsilon }^2{\parallel E_u^j\parallel }_H^2$), we instantiate the global system constant $C_9=max\left(1,{\displaystyle \frac{C_8}{{\epsilon }^2}}\left[1+{\displaystyle \frac{1}{2{\epsilon }^2}}+{\displaystyle \frac{C_1^4{\epsilon }^2}{2}}\right]\right)$. Grouping all exterior truncation errors and Taylor remainders into a unified base $b=C_5({\tau }^4+h^{12})$, the system flawlessly transforms into an explicit recurrence governed by $\mathcal{H}$:

$${\mathcal{H}}^{k+1}-{\mathcal{H}}^k+\tau {\parallel {\delta }_t{E}_u^{k+1}\parallel }_{H^{-1}}^2\le \tau C_9\left({\mathcal{H}}^{k+1}+{\mathcal{H}}^k+{\mathcal{H}}^{k-1}\right)+\tau b.$$

(50)

Since $\tau \parallel {\delta }_t{E}_u^{k+1}{\parallel }_{H^{-1}}^2\ge 0$, we obtain the weaker inequality without this term. Rearranging the terms configures the implicit relation $(1-\tau C_9){\mathcal{H}}^{k+1}\le (1+\tau C_9){\mathcal{H}}^k+\tau C_9{\mathcal{H}}^{k-1}+\tau b$. Imposing a temporal constraint $\tau \le {\displaystyle \frac{1}{2C_9}}$ guarantees $1-\tau C_9\ge 1/2>0$. Dividing both sides by this factor and subtracting ${\mathcal{H}}^k$, the identity ${\displaystyle \frac{1+\tau C_9}{1-\tau C_9}}-1={\displaystyle \frac{2\tau C_9}{1-\tau C_9}}$ yields:

$${\mathcal{H}}^{k+1}-{\mathcal{H}}^k\le {\displaystyle \frac{2\tau C_9}{1-\tau C_9}}{\mathcal{H}}^k+{\displaystyle \frac{\tau C_9}{1-\tau C_9}}{\mathcal{H}}^{k-1}+{\displaystyle \frac{\tau b}{1-\tau C_9}}.$$

(51)

Applying the summation operator ${\sum }_{k=1}^m$ triggers a telescoping sum ${\sum }_{k=1}^m({\mathcal{H}}^{k+1}-{\mathcal{H}}^k)={\mathcal{H}}^{m+1}-{\mathcal{H}}^1$. For the lagged sum, we define the initial state for the fictitious step $E_u^{-1}=0,e_r^{-1}=0$, which, coupled with the exact physical initialization $E_u^0=0,e_r^0=0$, algebraically enforces ${\mathcal{H}}^0\equiv 0$. Consequently, ${\sum }_{k=1}^m{\mathcal{H}}^{k-1}\le {\sum }_{j=0}^m{\mathcal{H}}^j$ holds strictly. Bounding the multipliers via ${\displaystyle \frac{3\tau C_9}{1-\tau C_9}}\le 6\tau C_9$ and the constant sum via ${\displaystyle \frac{1}{1-\tau C_9}}\le 2$ condenses the expression to:

$${\mathcal{H}}^{m+1}\le \left({\mathcal{H}}^1+2Tb\right)+\sum _{j=0}^{m}6C_9\tau {\mathcal{H}}^j.$$

(52)

In accordance with Assumption 6, squaring the initial constraints provides the bound ${\mathcal{H}}^1\le C_{11}{({\tau }^2+h^6)}^2$, where $C_{11}$ is derived from $C_{10}$. Given the condition ${\tau }^4+h^{12}\le {({\tau }^2+h^6)}^2$, the constant terms fuse perfectly into $(C_{11}+2T{C}_5){({\tau }^2+h^6)}^2$. By the discrete Gronwall lemma,we have :

$${\mathcal{H}}^{m+1}\le exp\left(6C_9{t}_{m+1}\right)(C_{11}+2T{C}_5){({\tau }^2+h^6)}^2.$$

(53)

Recalling the coercivity of the energy functional ${\mathcal{H}}^{m+1}\ge {\epsilon }^2\parallel E_u^{m+1}{\parallel }_H^2\ge {\epsilon }^2{C}_1^{-2}{\parallel E_u^{m+1}\parallel }_{H^{-1}}^2$. Executing the square root operation symmetrically reduces the exponent $6C_9$ to $3C_9$. Defining $C_{12}=3C_9$ and $C_{13}={\epsilon }^{-1}{C}_1\sqrt{C_{11}+2T{C}_5}$, we explicitly recover the optimal convergence precision:

$$\parallel u\left(t_{m+1}\right)-U^{m+1}{\parallel }_{H^{-1}}\le C_{13}exp\left(C_{12}T\right)\left({\tau }^2+h^6\right).$$

(54)

This completes the proof. □

Remark. 

To verify the a priori uniform boundedness assumption ($\parallel E_u^{m+1}{\parallel }_{L^{\infty }}\le 1$), we apply the discrete inverse inequality:

$$\parallel E_u^{m+1}{\parallel }_{L^{\infty }}\le C_{14}{h}^{-(d+\alpha )/2}{\parallel E_u^{m+1}\parallel }_H.$$

(55)

Enforcing the strict asymptotic mesh coupling constraint ${\tau }^2+h^6=o\left(h^{(d+\alpha )/2}\right)$ assures that in the infinitesimal limit as $\tau ,h\to 0$, the mesh singularity amplification is completely dominated by the high-order convergence rate. Consequently, the maximum norm error is compelled to zero, ensuring $\parallel E_u^{m+1}{\parallel }_{L^{\infty }}\to 0\le 1$. This validates authenticates the assumption and resolves the stiff nonlinearities.

3.3. Numerical Experiments

To see how the algorithm actually works, this subsection runs numerical tests on the fractional Cahn-Hilliard equation from two main angles.

Table 1. Spatial convergence orders of the proposed sixth-order scheme under highly anisotropic configurations in 3D, 4D, and 5D spaces.

Dim	Config ($\mathit{\alpha }$)	h Refinement	Max Error (${\mathit{L}}_{\mathit{\infty }}$)	Order
3D	$[1.2,1.5,1.8]$	$1/32\to 1/64$	3.41485e-03	—
		$1/64\to 1/128$	5.62345e-05	5.9242
		$1/128\to 1/256$	8.90324e-07	5.9810
4D	$[1.2,1.5,1.8,1.3]$	$1/32\to 1/64$	3.57993e-03	—
		$1/64\to 1/128$	5.88997e-05	5.9255
		$1/128\to 1/256$	9.32308e-07	5.9813
5D	$[1.2,1.5,1.8,1.3,1.7]$	$1/32\to 1/64$	5.22353e-03	—
		$1/64\to 1/128$	8.59100e-05	5.9261
		$1/128\to 1/256$	1.35972e-06	5.9814

Table 1. Spatial convergence orders of the proposed sixth-order scheme under highly anisotropic configurations in 3D, 4D, and 5D spaces.

Dim	Config ($\mathit{\alpha }$)	h Refinement	Max Error (${\mathit{L}}_{\mathit{\infty }}$)	Order
3D	$[1.2,1.5,1.8]$	$1/32\to 1/64$	3.41485e-03	—
		$1/64\to 1/128$	5.62345e-05	5.9242
		$1/128\to 1/256$	8.90324e-07	5.9810
4D	$[1.2,1.5,1.8,1.3]$	$1/32\to 1/64$	3.57993e-03	—
		$1/64\to 1/128$	5.88997e-05	5.9255
		$1/128\to 1/256$	9.32308e-07	5.9813
5D	$[1.2,1.5,1.8,1.3,1.7]$	$1/32\to 1/64$	5.22353e-03	—
		$1/64\to 1/128$	8.59100e-05	5.9261
		$1/128\to 1/256$	1.35972e-06	5.9814

4. Discussion

4.1. Accuracy Verification, Energy Stability, and Phase Morphologies

This subsection analysis the accuracy and physical reliability of the proposed scheme. First, the spatial convergence of the sixth-order central difference scheme under highly anisotropic fractional operators is verified. Subsequently, the unconditional energy stability is tested. The computational domain is set to a 3D space $\Omega ={(0,1)}^3$ with a grid resolution of $M=32$. The system is initialized with a multi-dimensional sinusoidal perturbation $u_0(x,y,z)=0.05sin\left(2\pi x\right)sin\left(2\pi y\right)sin\left(2\pi z\right)$. The interface width parameter is set to ${ϵ}^2=0.04$, and the time step size is fixed at $\tau =0.01$. The simulation is advanced to a terminal time $T=3.0$ to capture the delayed spinodal decomposition process induced by extreme fractional anisotropy.

Checking if the physics hold up involves running six specific fractional setups in a 3D grid. These chosen cases cover symmetric setups that act the same in all directions (e.g., ${\alpha }_1={\alpha }_2={\alpha }_3=1.5$) alongside setups that change depending on the direction (e.g., ${\alpha }_1=1.2,{\alpha }_2=1.5,{\alpha }_3=1.8$). Our experiments keep track of the modified discrete energy. Across every single test, the energy levels fall sharply when the phase separation first starts, and then they level out into a slow slide during the coarsening stage. Importantly, having different fractional parameters does not break this steady downward trend. That behavior closely matches the theoretical energy dissipation rule, as verified in Figure 1. Figure 2 points out the vastly different 3D phase separation patterns driven by these varying parameters.

4.2. Underlying Execution Efficiency Analysis: VTP vs. DST Under the $\tau$-Preconditioner

This subsection directly compares the underlying single-step execution efficiency of the VTP architecture against the conventional Standard DST within the identical mathematical framework of the $\tau$-preconditioner, isolating the core computational operation of the Krylov solver ($z=P^{-1}x$). The traditional Standard DST relies on multi-layer nested loops to extract data fibers. As the system dimension d increases, the memory strides generated by this non-contiguous access rapidly exceed the physical capacity of CPU caches. In contrast, the VTP architecture utilizes continuous tensor permutation to reshape the computational dimensions into contiguous memory blocks, facilitating the SIMD vectorization pipelines of modern processors.

As presented in Table 2, the speedup ratio reflects the dynamic shifting of hardware bottlenecks. Because the traditional DST extracts short data fibers sequentially, severe SIMD pipeline starvation and deep nested loop overheads dominate its computation time. VTP mitigates this loop overhead via matricization, reaching a speedup of $35.55\times$ in the 5D space with a minimal grid ($M=4$). However, as the grid scaling induces exponential data expansion, the dense tensor permutation operations in VTP begin to encounter the theoretical upper limits of memory bandwidth. Consequently, when the grid increases to $M=8$, the single-step speedup stabilizes in the range of $14.27\times \sim 25.20\times$. Despite this hardware-imposed constraint on memory throughput, VTP strictly maintains contiguous memory access and continues to exhibit a substantial performance advantage.

Furthermore, Table 2 indicates the decoupling characteristic of the anisotropic parameter $\alpha$ on hardware execution efficiency. Under the same grid (e.g., 4D, $M=12$), whether $\alpha$ is symmetric $[1.5,\dots ,1.5]$ or highly asymmetric $[1.2,1.4,1.6,1.8]$, the single-step execution time of VTP fluctuates within a narrow interval of $0.0026\sim 0.0039$ seconds. The fundamental reason is that changing the directionality of $\alpha$ merely alters the values in the diagonal preconditioner matrix; the underlying instruction graph remains physically unchanged.

As illustrated in Figure 3, when approaching the extreme spaces from 6D to 8D, the execution time of the traditional DST architecture exhibits an exponential deterioration trend due to cache degradation. In contrast, the VTP architecture maintains a steady polynomial-level scaling in the log-log coordinate system, validating its hardware robustness against high-dimensional scaling.

4.3. Total Solver Time Comparison: $\tau$-Preconditioned VTP vs. Baseline Algorithm

The total solver time is synergistically determined by the single-step execution time and the algebraic Krylov iteration count.

Pitting the $\tau$-preconditioned VTP setup against a baseline skew-circulant preconditioner forms the core of this subsubsubsection.Boundary disturbances near the edges drive up the condition number. Such boundary mismatches force the iteration count to climb visibly as the grid gets finer, a trend clearly captured in Table 3.Mathematical alignment with the homogeneous Dirichlet boundaries comes naturally to the $\tau$-preconditioner.

Algebraically, it safely fences in the eigenvalues of the preconditioned matrix. The Krylov steps stay pinned at a rock-bottom level of roughly 15. This low count holds steady no matter what you do with the dimension d, the grid size M, or the fractional parameters $\alpha$.

Table 3 puts hard numbers on the VTP acceleration.In the 4D tests with $M=12$, the baseline iteration count gets worse and falls back to 29 steps. Combined speedup ratios for the VTP framework land anywhere from $55.0\times$ up to $73.2\times$. To our disappointment, pushing M from 4 to 8 in the 5D runs causes the peak speedup to slide from $61.1\times$ down to $53.2\times$. Data grows exponentially alongside grid refinement. Heavy tensor shuffling starts hitting the physical wall of memory bandwidth. Those bare-metal limits inevitably pull the acceleration back a bit. Giving up a little ground to the hardware only eats up a tiny fraction of the time saved by taking fewer steps. Our VTP setup easily holds onto a clear performance lead over the older baseline. Interestingly, picking drastically different $\alpha$ values for each spatial axis exposes another weakness in the baseline. Raw execution times for VTP stay entirely flat. That older baseline reacts badly to the messy spectral gaps caused when these $\alpha$ numbers do not match.

Figure 4 maps out the scaling behavior from 6D all the way up to 8D. Under these brutal setups, total execution times for the baseline spiral out of control rather than simply crashing from RAM capacity limits. Choking on severe cache thrashing and mounting iteration penalties, that older solver eventually forces a reliance on theoretical time projections just to finish the evaluation. Strict log-linear scaling naturally emerges from our VTP design. It delivers a highly workable engine to push through massive fractional PDEs without hitting those same hardware walls.

5. Conclusion

This paper set out to build an optimized, highly efficient fully discrete algorithmic framework for the multi-dimensional space-fractional Cahn-Hilliard equation. To capture the nonlocal phase separation physics with enough accuracy, we put together a high-order numerical scheme. The temporal side uses a second-order SAV-BDF2 discretization; the spatial side employs a sixth order Riesz fractional centered difference method that we derived from a pair of generating functions. We supplied rigorous theoretical proofs to back up two claims: the scheme is uncon ditionally energy stable, and it hits optimal spatio-temporal convergence rates inside the fully discrete setting. Beyond the discretization itself, we also had to deal with the severe ill-conditioning and the memory access bottlenecks that come with high-dimensional fractional operators. For that, we implemented a preconditioned MINRES solver and introduced a hardware-level optimization: the Vectorized Tensor Processing (VTP) architecture. The VTP idea is to reshape the tensors on the fly, turning multi-dimensional preconditioned matrix-vector products into batched trans form operations that access memory contiguously. This pretty much wipes out cache thrashing. Extensive numerical experiments checked the theoretical convergence orders, confirmed that the discrete energy really does dissipate, and showed the substantial speedup the VTP-accelerated solver gets over conventional implementations. Taken together, these results back up the effi ciency and the robustness of the framework we have laid out.