## 2 Numerical Method  
In this section, the Fourier spectral method is introduced for the proposed model. This method offers an efficient and highly accurate computational method for solving partial differential equations, particularly when the equations are defined on periodic domains. This method expands the solution in terms of trigonometric basis functions, converting spatial derivatives into simple algebraic operations in Fourier space. By using the Fast Fourier Transform (FFT), it significantly reduces computational cost while preserving spectral accuracy. This method is especially effective for problems with smooth solutions and offers exponential convergence. However, its use is limited in some specific domains unless modified with appropriate extensions or combined with other techniques. The Fourier spectral method has seen broad applications in simulating fluid dynamics, phase-field models, and nonlinear pattern formation [31]. Now, we describe a computational scheme for Equation (4) in the domain \(\Omega=(L_{x},\ R_{x})\times(L_{y},\ R_{y})\), which is discretized as \(\Omega_{d}=\{(x_{m},y_{n})|x_{m}=L_{x}+(m-0.5)h_{x},\)\(y_{n}=L_{y}+(n-0.5)h_{y},\)\(1\leq m\leq N_{x},\)\(1\leq n\leq N_{y}\}.\) Let \(\phi_{mn}^{k}=\phi(x_{m},y_{n},t_{k})\), where \(t_{k}=k\Delta t\), and \(\Delta t\) is the temporal step size. The discrete Fourier transform is defined as  
\[\hat{\phi}_{pq}^{k}=\sum_{m=1}^{N_{x}}\sum_{n=1}^{N_{y}}\phi_{mn}^{k}e^{-i( \xi_{p}x_{m}+\eta_{q}y_{n})},\ \ -\frac{N_{x}}{2}+1\leq p\leq\frac{N_{x}}{2},\ \ -\frac{N_{y}}{2}+1\leq q\leq \frac{N_{y}}{2}, \tag{6}\]  
where \(\xi_{p}=2\pi p/L_{x}\), \(\eta_{q}=2\pi q/L_{y}\) and \(\phi_{mn}^{k}\) denotes the numerical approximation of \(\phi(x_{m},y_{n},t_{k})\). The inverse discrete Fourier transform can be expressed as follows:  
\[\phi_{mn}^{k}=\frac{1}{N_{x}N_{y}}\sum_{p=-\frac{N_{x}}{2}+1}^{\frac{N_{y}}{ 2}}\sum_{q=-\frac{N_{y}}{2}+1}^{\frac{N_{y}}{2}}\phi_{pq}^{k}e^{i(\xi_{p}x_{m} +\eta_{q}y_{n})}. \tag{7}\]Let \(\phi(x,y,t)\) be a continuous version of \(\phi^{k}_{mn}\). Equation (8) represents the Fourier series expansion of the function \(\phi(x,y,t)\):  
\[\phi(x,y,t)=\frac{1}{N_{x}N_{y}}\sum_{p=-\frac{N_{x}}{2}+1}^{\frac{N_{y}}{2}}\sum _{q=-\frac{N_{y}}{2}+1}^{\frac{N_{y}}{2}}\hat{\phi}(\xi_{p},\eta_{q},t)e^{i(\xi_ {p}x+\eta_{q}y)}. \tag{8}\]  
Based on this spectral representation, both the Laplacian \(\Delta\phi(x,y,t)\) and the bi-Laplacian \(\Delta^{2}\phi(x,y,t)\) can be conveniently computed in Fourier space through simple multiplicative operations. Accordingly, Equations (9) and (10) express the inverse transforms that yield the Laplacian and bi-Laplacian in the physical space:  
\[\Delta\phi(x,y,t) = \frac{1}{N_{x}N_{y}}\sum_{p=-\frac{N_{x}}{2}+1}^{\frac{N_{x}}{2}} \sum_{q=-\frac{N_{y}}{2}+1}^{\frac{N_{y}}{2}}-(\xi_{p}^{2}+\eta_{q}^{2})\hat{ \phi}(\xi_{p},\eta_{q},t)e^{i(\xi_{p}x+\eta_{q}y)}, \tag{9}\] \[\Delta^{2}\phi(x,y,t) = \frac{1}{N_{x}N_{y}}\sum_{p=-\frac{N_{x}}{2}+1}^{\frac{N_{x}}{2}} \sum_{q=-\frac{N_{y}}{2}+1}^{\frac{N_{y}}{2}}(\xi_{p}^{2}+\eta_{q}^{2})^{2} \hat{\phi}(\xi_{p},\eta_{q},t)e^{i(\xi_{p}x+\eta_{q}y)}. \tag{10}\]  
To efficiently solve Equation (11), we use the linearly stabilized splitting method [32], which is known to improve numerical stability and allow for larger time steps without compromising the stability of the solution. The linearly stabilized splitting method is a semi-implicit time integration scheme designed for stiff or nonlinear PDEs, especially gradient-flow-type models. It adds and subtracts a linear stabilization term to the nonlinear component, which allows the linear part to be treated implicitly and the remaining nonlinear part explicitly. This approach improves stability and permits larger time steps without sacrificing accuracy. The method is particularly effective for phase-field equations such as the Allen-Cahn and Cahn-Hilliard models, and offers unconditional energy stability and computational efficiency by avoiding nonlinear solvers.