Structure-Preserving Adaptive Spectral–Finite Element Method for Nonlocal Quantum–Fractional PDEs with Variable Coefficients

Full text

Text extracted from the article PDF for reading and search. Tables, figures and equations are in the PDF, which remains the authoritative version.

1. Introduction

Mathematical modeling through nonlocal quantum fractional partial differential equations (PDEs) that employ variable coefficients has become a rapidly growing methodology for the development of highly accurate representations of complex systems characterized by memory effects, spatial nonlocality, heterogeneous materials, and interactions at the scale of quantum mechanics. Whereas traditional integer-order PDEs can only model local behavior through differential equations that relate to the present state of a system, nonlocal and fractional PDEs allow for an extension of the concept of locality to include both long-range and historical effects that can impact the present state of a system. Such features make these types of PDEs particularly applicable to modeling the following: anomalous diffusion, quantum transport, fractional Schrödinger-type dynamics, propagation of waves in non-locality, complex materials and dissipative systems where the physical characteristics of matter vary as a function of space and time (Atman & Şirin, 2020; Godinho & Vancea, 2025; Mostafanejad, 2021). The use of variable coefficients and variable-order operators in these models also adds greater physical realism, since many quantum and fractional systems develop in (chemical) environments that exhibit variable degrees of diffusion, present different types of potential energy fields, respond to external forces nonlinearly, or have fractional orders that vary with both location and time (Kheybari, 2021; Xiang et al., 2021; Zhao et al., 2019). Ultra-quick convergence of smooth input types allows for efficient applications of spectral techniques, particularly when dealing with overall operators, multiple branches, and high-regulation output types (ZRJ & ZHO & TGR, 2021; AGG & DWS, 2018; ZR & GES, 2020). Unless you only need to process using a limited set of geometrically identifiable objects, spectral approaches provide you fewer limitations regarding form, i.e. allow you to be at a generic position (DZ & APP, 2017; DA & ET & EL, 2017; DA & SUE, 2018). However, you will have to adapt to changes in your application domain when using these hybrid implementations. Thus, regardless of whether you are implementing a hybrid add-on or a more "traditional" version, the end result should be that you will not only obtain significant advances in solving your application, but you will also find that your method will perform extremely different in regards to how the various structures actually function. Many problems demonstrate irregular behavior at different rates because of local anomalies, local nonlinear behavior, as well as nonuniformity in coefficient and/or quantum nature; thus, refining the resulting numerical solutions will generally not yield a meaningful-cost solution until effective h and/or p adaptive methods are successfully employed (Bulle etal 2023; Dorfler 1996; Sun etal 2007). For spectral/finite element discretizations, the use of adaptive refinement techniques relies on adaptations based upon a posteriori error estimates, providing a systematic approach for determining which area is in need of refinement while keeping the cost of computation under control and still achieving accurate results (Bulle etal 2023; Dorfler 1996; Sun etal 2007). For example, when utilizing spectral/finite element discretizations, h refinement, p refinement, hp refinement, basis shifting, size scaling, domain partitioning, and/or local enhancement of the approximation space can all potentially be incorporated into adaptive refinements based upon previous error estimates (Chou etal 2023; Feng & Mavriplis 2002; Mavriplis 1994; Valenciano & Owens 2000).

Certification is a method of achieving numerical error control using rigorous high standards in error estimation, convergence theory and stop criteria. Many adaptive finite element methods in the literatures provide evidence of convergence, reliability, efficiency, and optimality with respect to several classes of elliptic and associated PDEs (Binev et al. 2004; Dörfler 1996; Stevenson 2007). More recent work has utilized certification for the development of inexact solvers, nonlinear PDEs, fractional passive models (Laplace models), eigenvalue problems, and adaptive algorithms that maintain quasi-optimal computational complexities (Bringmann et al. 2024; Bringmann et al. 2025; Bulle et al. 2023; Giani & Graham 2009).

Nonlocal operators introduce global coupling, singular integral terms, dense matrix structures, and boundary interactions that are not captured by standard local residual estimators for fractional Laplacian and time-fractional parabolic problems, specialized a posteriori estimators and reconstruction-based techniques are required to obtain reliable error control (Ainsworth & Glusa, 2017; Bulle et al., 2023; Cao et al., 2025). In addition, nonlinear variable-order fractional PDEs introduce further complications because the differential order itself change over the computational domain, requiring adaptive approximation of both the solution and the operator structure (Azghay et al., 2025; Xiang et al., 2021; Zhao et al., 2019).

Quantum and fractional systems possess conservation laws, stability constraints, orthogonality conditions, or energy structures that must be respected by the numerical scheme If these properties are not preserved, a highly accurate approximation in the algebraic sense still produce physically misleading results over long-time simulations. Structure-preserving discretizations have therefore become important in computational quantum mechanics, nonlocal diffusion, Kohn–Sham density functional theory, and fractional Schrödinger-type models (Du et al., 2013; Hendy & Zaky, 2020; Motamarri et al., 2013; Zhan et al., 2025) which studies shwed that energy-stable and orthonormality-preserving schemes in quantum models demonstrate the importance of aligning numerical discretization with the mathematical structure of the governing equation (Zhan et al., 2025). Similarly, penalty methods and spectral element formulations provide mechanisms for enforcing boundary and interface conditions while maintaining stability in high-order discretizations (Hesthaven, 2000; Song & Jung, 2020). Operational matrix approaches have been used to transform fractional PDEs and nonlocal constraints into algebraic systems, offering computational efficiency and simplified implementation for certain classes of problems (Khalil et al., 2019; Khalil et al., 2021). Multi-domain spectral collocation and domain decomposition techniques provide additional flexibility by dividing the computational domain into subregions and enforcing coupling through interface conditions, penalty terms, or nonlocal constraints (Feng & Mavriplis, 2002; Zhao et al., 2019).

Recent advances in variational quantum algorithms have introduced an additional dimension to the numerical solution of PDEs which the variational quantum approaches have been proposed for nonlinear and multidimensional PDEs and eventually provide advantages for selected high-dimensional or computationally intensive problems (Nguyen, 2025; Sarma et al., 2024). However, their use in nonlinear nonlocal quantum–fractional PDEs with variable coefficients remains at an early stage.

The aim of this study is to develop and analyze a certified, structure-preserving adaptive spectral–finite element method for solving nonlocal quantum–fractional partial differential equations with variable coefficients. In addition, the study aims to preserve the essential mathematical and physical structures of quantum–fractional models, including stability, nonlocal consistency, and boundary-condition integrity, particularly in the presence of variable coefficients, variable-order operators, and nonlinear nonlocal interactions.

2. Mathematical Model and Problem Formulation

Let Ω ⊂ ℝᵈ, d = 1, 2, or 3, be a bounded Lipschitz domain with boundary ∂Ω. We consider a complex-valued wave function ψ: Ω × [0,T] → ℂ governed by a nonlocal quantum– fractional PDE of the form:

4. Proposed Certified Structure-Preserving Adaptive

Spectral–Finite Element Method To resolve the computational challenges induced by the nonlocal fractional operator and to guarantee physical fidelity over long-term simulations, we propose a novel Certified Structure-Preserving Adaptive Spectral–Finite Element Method (CSPA-SFEM). This framework blends the geometric flexibility of finite element methods with the diagonalizing efficiency of spectral calculus, backed by a fully computable, rigorous error estimator.

4.1. Spatial Finite Element Discretization

Let ??????ℎbe a conforming, shape-regular simplicial triangulation of the bounded domain �⊂ℝ??????. We construct the standard continuous piecewise polynomial finite element space ??????ℎ⊂ �01(�)of degree ??????≥1associated with ??????ℎ , defined as:

This construction naturally preserves the self ad jointness, symmetry, and positive definiteness of the continuous fractional operator at the algebraic level, bypassing the dense matrix obstructions encountered in traditional integral fractional Laplacians.

4.2. Spectral Truncation

In practical high-dimensional computations, computing the full spectrum up to �ℎis computationally prohibitive. Therefore, we introduce a spectral truncation parameter ??????≪ �ℎretaining only the lowest $m$ dominant eigenmodes. The truncated discrete fractional operator ℒ??????,ℎ,��:??????ℎ→??????ℎis formulated as:

To certify the reliability of this approximation within our adaptive loop, the spectral truncation error is controlled by the computable high-frequency energy tail estimator �spec defined as:

This explicit component guarantees that our adaptive refinement mechanism monitors both spatial grid partition and spectral mode enrichment simultaneously.

4.3. Structure-Preserving Time Integration

Correspondingly, for a standard Hartree kernel response where �(|??????|2)=|??????|2 the nonlocal interaction discrete gradient ??????ℎDG:??????ℎ×??????ℎ→??????ℎis structured symmetrically as:

4.4. Adaptive Certified Error Estimator

To automatically tune the discretization parameters, we derive a posteriori error estimator. The total computable error estimator $\eta_{\text{total}}$ is decomposed into distinct, orthogonal error contributions:

The explicit formulations for these estimators are defined as follows:

• Spatial Residual Estimator �space: Measures the discretization error from the spatial mesh grid triangulation, defined via element residuals ??????�and boundary jump residual �??????:

• Coefficient Approximation Estimator �coeff: Accounts for the interpolation or projection errors of the variable diffusion tensor ??????(�)on the grid:

• Temporal Residual Estimator (�time): Tracks the truncation error due to the time-stepping progress:

The Adaptive Loop Implementation The practical execution of the multi-adaptive configuration follows an extended adaptive feedback loop:

2. Estimate: Calculate all localized error indicators across each element K.

3. Mark: Apply the collective Dörfler marking criterion with a parameter �∈(0,1). We select a minimal subset of elements ℳℎ⊂??????ℎ such that:

4. Refine / enrich

• Spatial Refinement: Elements marked in ℳℎ are subdivided using standard newest-vertex bisection techniques. • Spectral Enrichment Strategy: If the spectral error dominates the spatial error, i.e �spec2≥ �spec�space2where �spec≈0.1 is a safety threshold), we expand the spectral parameter:

??????⟵??????+????????????, thereby adding higher-frequency eigenfunctions to the active simulation space without needing full spatial grid remeshing.

5. Update: Interpolate solutions onto the modified space, update the matrices, and advance to the next time frame.

5. Theoretical Analysis

In this section, we establish the fundamental theoretical foundations of the proposed Certified Structure-Preserving Adaptive Spectral–Finite Element Method (CSPA-SFEM). We rigorously prove the well-posedness, discrete conservation laws, a priori convergence, and the reliability and efficiency of the adaptive error estimator.

5.1. Well-Posedness of the Semi-Discrete and Fully

Discrete Problems To guarantee that the numerical trajectory is well-defined, we first address the existence and uniqueness of solutions in the finite-dimensional space ??????ℎ.

1. For any initial data ψ0∈Vh, the fully discrete implicit scheme proposed admits a unique solution ψhn+1∈Vhat each time steτn>0p .

2. Under subcritical nonlinear growth constraints 2s > d/2 3. for the cubic term), the local discrete solution can be uniquely extended globally in time up to the final simulation time .

Proof Sketch. The existence of a solution at each discrete time step follows from Brouwer's Fixed-Point Theorem applied to the finite-dimensional mapping on??????ℎ . Letℱ:??????ℎ→??????ℎ be the nonlinear operator representing the implicit step. By taking the inner product with the state vector and invoking the a priori mass bound (established in Section 5.2), we demonstrate that ℱ maps a closed ball into itself, ensuring existence. Uniqueness follows from the monotonicity of the local discrete gradient ??????ℎcombined with a discrete fractional Gronwall inequality under a mild constraint on the time step ??????�.

5.2. Rigorous Proof of Discrete Mass Conservation

A paramount feature of quantum systems is the preservation of the total probability density (mass).

We prove that our fully discrete scheme preserves this constraint exactly without accumulating drift over long-term dynamics Theorem 5.2 (Exact Discrete Mass Preservation). Let??????ℎ�+1 be the solution generated by the structure-preserving midpoint scheme. For any??????=0 and real-valued potentials ??????(�), the discrete mass is strictly invariant: Proof. We multiply the fully discrete governing equation by the complex conjugate of the midpoint solution,??????̅ℎ�+1/2=1 2(??????̅ℎ�+1+??????̅ℎ�) , and integrate over� . Taking the standard inner product yields:

5.3. Discrete Energy Stability and Solver Calibration

When standard schemes solve nonlinear PDEs, the energy often exhibits non-physical dissipation or blowing up due to truncation errors. The proposed discrete-gradient mechanism explicitly counters this.

Theorem 5.3 (Discrete Hamiltonian Conservation). Under identical conditions to Theorem 5.2, the fully discrete solution satisfies an exact energy conservation law at the continuous algebraic level:

Proof. The proof mirrors the mass preservation strategy but evaluates the real part Re{⋅}after testing the equation with the directional difference (??????ℎ�+1−??????ℎ�) . By construction, the local discrete gradient ??????ℎsatisfies the exact identity:

Inexact Solver Energy Defect Bounds In practical execution, the implicit algebraic system is solved using an iterative method (e.g., Newton-Raphson) terminated at a non-zero tolerance . We certify this non-exact energy preservation via the following computable theorem:

5.4. Convergence of the Fully Discrete Scheme

We present a comprehensive a priori error estimate that unifies the five interacting dimensions of discretization error in our framework.

• Spectral mode count (??????−2�/??????): Exhibiting algebraic or exponential convergence rates based on the regular smoothness parameters r.

5.5. Reliability and Efficiency of the Estimator

To establish that the adaptive loop behaves optimally, we prove that our combined total error estimator�total acts as a mathematically valid upper and lower bound for the true approximation error measured in the natural fractional energy norm‖⋅‖???????????? .

Reliability (Upper Bound)

The estimator is certified as reliable if it guarantees that the true computational error does not exceed the estimated tolerance up to a generic constant.

Proof Outline. By inserting the projection of the exact solution into the weak formulation, we establish the error-residual equation. Utilizing the global coercivity of the discrete fractional operator on �(ℒ??????�/2)and applying standard element-by-element interpolation bounds for the spatial residual and jump operators �??????, we bound the state error by the computable terms of �total. The coefficient and kernel deviations are handled via direct application of Cauchy-Schwarz inequalities.

Conversely, the estimator is efficient if it does not overestimate the actual error, preventing over-refinement of the grid or wasteful spectral mode allocation.

6. Numerical Results and Discussion

This section presents a benchmark numerical validation of the proposed Certified Structure-Preserving Adaptive Spectral-Finite Element Method (CSPA-SFEM) for the nonlocal quantum-fractional PDE formulated in the previous sections. The numerical tests were designed to evaluate four central properties of the method: approximation accuracy, certified a posteriori error control, preservation of discrete physical invariants, and computational efficiency under adaptive mesh and spectral enrichment. New developments in variational quantum algorithms have created another approach for numerically solving PDEs, with the potential to help in solving certain high-dimension or computationally intensive PDEs (Nguyen, 2025; Sarma et al., 2024). However, such approaches are only starting to be explored for use with nonlinear nonlocal quantum-fractional PDEs with variable coefficients. Unless stated otherwise, all numerical test cases were generated on the 1-D domain Ω = (0, 1) with final time T = 1. The fractional order was set s = 0.65, variable diffusion as a(x) = 1 + 0.25 sin (2πx), potential as V(x) = 1 + x², and cubic interaction coefficient as β(x) = 0.5 + 0.25x; with the nonlocal Hartree-type kernel defined using K(x, y) = exp[−25(x − y)²]. The initial wave function was a normalized Gaussian-type wave function with localized oscillatory modulation. The combination of these factors leads to a nonuniform spatial solution profile and therefore brings a suitable benchmark for spatially adaptive testing, spectral truncation controls and time integrators that preserve invariants.

6.1. Wave-function resolution and qualitative accuracy

To evaluate the final wave function amplitude | ψ(x,T) | using the proposed numerical method, a first numerical experiment was performed; this was done by comparison with a reference solution for the CSPA-SFEM method. Results were characterized by an excellent correlation between CSPA-SFEM calculations and the reference solution through the entire domain, including for the two-dimensional area containing the large localized concentrated areas of the greatest local residuals resulting from the variable coefficient and nonlocal interaction. The high degree of correlation between the two had been most demonstrated in the main region of concentration, indicating that the adaptive finite element refinement and combined spectral representations produced results that were representative of the global fractional behavior and spatial location of the localized variation.

\ Fig 1: Final wave-function amplitude obtained by the proposed CSPA-SFEM compared with the reference benchmark solution.

6.2. Convergence behavior under adaptive refinement

The use of an adaptive refinement scheme resulted in a consistent decrease in both the L² norm and the fractional energy norm errors. Monotonically decreasing errors were observed as the number of active spatial DOF increased, supporting that the new adaptive spectral-FEM formulation is stable and convergent for variable-coefficient quantum fractional dynamics. The observed converging behavior was in agreement with the theoretical estimate given in Section 5, where total error is dependent on the spatial discretisation size, contributions from the time step, and a truncation tail from the spectral truncation.

Table 1: Adaptive convergence and certified estimator performance for the benchmark problem. Cycle DOF Modes m L² error V_s error η_total Effectivity CPU time (s)

Fig 2: Decay of L² and fractional energy-norm errors as the adaptive degrees of freedom increase.

6.3. Reliability of the certified estimator

The approved estimator η_total stayed greater than the real fractional energy-norm error during all adaptive cycles and therefore maintained an effective, consistent efficiency level. The range of effectiveness was approximately 1.21 to 1.24, indicating the estimator did not overestimate the error by a significant amount, thus not producing excessive amounts of refinement. This behavior supports the proposed certification system as being useful and reliable.

Fig 3: Comparison between the true benchmark error and the total certified estimator across adaptive cycles.

6.4. Contribution of estimator components

The total η decomposition revealed that at the beginning of the adaptive cycles, the largest area for error was related to errors due to spatial residuals. After refinement by local refinement, the spatial residual is reduced, and it began to show the component corresponding to spectral truncation, which activated the spectral enrichment step (m ← m+Δm). At all stages, the contributions from time, coefficient, kernel and quadrature remained controlled, indicating that the estimator has separated the different numerical error mechanisms, rather than viewing global error as a single non-diagnostic entity.

Fig 4: Decomposition of the total adaptive estimator into spatial, spectral, temporal, coefficient, kernel, and quadrature components.

6.5. Preservation of mass and Hamiltonian energy

The presence of structure-preserving Midpoint and Discrete Gradient approaches ensured that the discrete mass and Hamiltonian energy remained approximately equal to an acceptable solver tolerance through the entire time period. Standard non-preserving approaches experienced significant incremental drift of the same invariant; its relative mass drift was around 10⁻¹² and its relative energy drift was about 10⁻¹¹. This further validates the long-term viability of conducting Quantum/Fractional simulations using this method as small algebraically cumulative errors have the potential to result in physically unjust enacts.

Fig 5: Relative mass drift of the proposed method compared with a standard non-preserving scheme.

Fig 6: Relative Hamiltonian energy drift of the proposed method compared with a standard non-preserving scheme.

6.6. Local adaptivity and residual localization

The localized residual indicators identified the regions requiring refinement without uniformly increasing the mesh density over the entire domain. The largest residuals appeared near the localized wave-packet peaks and the regions where the nonlocal kernel generated stronger interaction gradients. After adaptive refinement, the magnitude of the local indicators decreased substantially across the domain, with the strongest reduction occurring near the initially dominant residual peaks. This confirms that the adaptive loop SOLVE– ESTIMATE–MARK–REFINE/ENRICH–UPDATE effectively concentrated computational effort where it was mathematically needed.

Fig 7: Local residual indicator before and after adaptive refinement.

6.7. Computational efficiency

According to the results from our computation scaling trials, the CPU times associated with the adaptive strategy increase sub-quadratically as a function of the active degree of freedoms (DOF). This behaviour is due to using a combination of local mesh refinement and spectral enrichment where spatial (local) refinement will only occur in areas with large residuals, and spectral (global) enrichment will only occur when the estimator for the high-frequency tail indicates insufficient modal resolution. Thus, the method yields increased accuracy without the expense of using complete uniform mesh refinement or complete spectral computation.

Fig 8: Computational scaling of the proposed method across adaptive cycles.

Based on the numerical results, the proposed CSPA-SFEM provides an accurate, stable and computationally efficient modelling framework for variable-coefficient nonlocal quantum-fractional PDEs. Through this method, monotonic reductions in error have been achieved, the behaviour of the estimator is confirmed as reliable, the manner in which spectral truncation is controlled has been demonstrated, and discrete physical invariants are preserved with close to exact accuracy. These results corroborate the previously stated theoretical results, thus illustrating that the proposed methodology is well suited to problems having the characteristics of nonlocality, fractional diffusion, variable coefficients, and long-time quantum stability.

7. Discussion

The numerical findings of the present study demonstrate that the proposed Certified Structure-Preserving Adaptive Spectral–Finite Element Method (CSPA-SFEM) provides a robust and mathematically consistent framework for approximating nonlocal quantum–fractional partial differential equations with variable coefficients. The observed stability, convergence, error reduction, and preservation of discrete invariants are consistent with the theoretical motivation behind combining adaptive finite element discretization, spectral fractional approximation, and structure-preserving time integration. This is particularly important because fractional and nonlocal quantum models are known to involve long-range interactions, memory effects, dense operator structures, and nonstandard boundary behavior, which cannot be treated efficiently by conventional local discretization techniques alone (Atman & Şirin, 2020; Godinho & Vancea, 2025; Mostafanejad, 2021). The qualitative behavior of the computed wave function showed that the proposed scheme was able to capture both global nonlocal effects and localized spatial variations caused by heterogeneous coefficients. This agrees with the broader understanding that fractional and nonlocal operators introduce spatial interactions beyond neighboring mesh elements, making the solution sensitive to both global operator structure and local coefficient variation (Ainsworth & Glusa, 2017; Du et al., 2013). The stable profile of the wave-function amplitude indicates that the spectral representation of the fractional operator contributes to numerical smoothness and reduces artificial oscillations, which is consistent with the advantages of spectral and spectral/hp methods in high-accuracy approximation of smooth or semi-smooth solution structures (Hafeez & Krawczuk, 2024; Xu et al., 2018).

The convergence results confirm that the proposed method benefits from the complementary strengths of finite element and spectral discretization. The finite element component provides geometric flexibility and local mesh refinement, while the spectral component offers an efficient representation of fractional powers of elliptic operators. Similar advantages have been reported in adaptive finite element methods for fractional Laplacian problems, where a posteriori estimators and efficient solvers are essential for controlling the complexity of the nonlocal operator (Ainsworth & Glusa, 2017; Bulle et al., 2023). The decrease in numerical error under mesh refinement also agrees with classical adaptive finite element theory, which establishes that properly designed adaptive schemes can achieve convergence and near-optimal approximation rates (Binev et al., 2004; Dörfler, 1996; Stevenson, 2007). The role of spectral enrichment was also evident in the numerical results. At early refinement levels, the spectral truncation component contributed noticeably to the total error, indicating that the number of active eigenmodes is a decisive factor in the approximation of fractional operators. This agrees with previous work on fractional spectral collocation and rational spectral methods, which shows that spectral representations are highly effective for fractional PDEs but require appropriate mode selection to balance accuracy and computational cost (Tang et al., 2018; Tang et al., 2020; Zayernouri & Karniadakis, 2014). Therefore, the proposed spectral enrichment criterion is important because it prevents under-resolution of the fractional operator while avoiding the unnecessary computation of excessive high-frequency modes. The reliability of the a posteriori estimator represents one of the most important outcomes of the study. The numerical comparison between the true error and the estimated error showed that the estimator followed the same decay trend and provided a practical upper bound for the numerical error. This result supports the certified nature of the proposed framework and is consistent with modern adaptive methods for fractional and nonlocal problems, where standard local residual indicators alone are often insufficient to capture all error sources (Ainsworth & Glusa, 2017; Bulle et al., 2023; Cao et al., 2025). In particular, nonlocal and fractional PDEs require estimators that account not only for spatial residuals but also for spectral truncation, kernel approximation, quadrature effects, and coefficient heterogeneity.

The decomposition of the total error estimator further confirms that nonlocal quantum–fractional PDEs involve multiple interacting sources of numerical error. The dominance of spatial and spectral components at coarse levels suggests that both the physical domain discretization and the fractional operator approximation must be refined together. This supports the use of hybrid adaptive strategies rather than relying exclusively on uniform mesh refinement or fixed spectral truncation. Similar multi-component challenges have been emphasized in variable-order fractional models and nonlinear fractional PDEs, where the operator structure itself vary in space or time and therefore requires flexible numerical treatment (Kheybari, 2021; Xiang et al., 2021; Zhao et al., 2019).

The preservation of discrete mass is a central strength of the proposed method. In quantum-type systems, the (L^2)-norm of the wave function is associated with probability conservation; therefore, numerical drift in mass lead to physically misleading simulations. The present results showed that the midpoint structure-preserving formulation maintained the mass almost exactly throughout the simulation. This behavior is consistent with the principle that numerical schemes for Schrödinger-type and quantum models should preserve the fundamental invariants of the continuous system whenever possible (Hendy & Zaky, 2020; Zhan et al., 2025). The result also demonstrates that the proposed discretization does not merely approximate the algebraic form of the PDE, but respects its underlying physical structure.

The Hamiltonian energy results provide additional evidence of the structure-preserving capability of the method. The nearly constant discrete energy over time indicates that the discrete-gradient treatment of the local cubic nonlinearity and the symmetric handling of the nonlocal interaction successfully prevent artificial energy dissipation or growth. This is particularly important for long-time simulations, where small non-physical energy errors can accumulate and distort the numerical dynamics. Energy-stable and invariant-preserving schemes have been shown to be essential in computational quantum mechanics and Kohn–Sham-type models, especially when nonlinear interactions and high-order spatial discretizations are involved (Motamarri et al., 2013; Zhan et al., 2025).

The adaptive local error distribution showed that the algorithm effectively identified regions requiring refinement. Before adaptation, high-error regions were concentrated near areas of strong spatial variation, coefficient transitions, and nonlocal interaction effects. After refinement, the error became more balanced across the computational domain. This behavior agrees with the established principle of adaptive finite element methods: refinement should be concentrated where the local contribution to the global error is largest, rather than being distributed uniformly across the entire domain (Dörfler, 1996; Sun et al., 2007). Such localization is especially valuable in fractional and nonlocal models, where uniform refinement be computationally expensive due to dense operator coupling (Ainsworth & Glusa, 2017; Du et al., 2013).

The computational scaling results indicate that the proposed adaptive method is more efficient than a purely uniform refinement strategy. This improvement is expected because the algorithm balances spatial refinement, spectral enrichment, and time-step control based on certified error information. Modern adaptive finite element theory has shown that adaptivity can lead not only to convergence but also to quasi-optimal computational runtime when combined with appropriate iterative solvers and stopping criteria (Bringmann et al., 2024; Bringmann et al., 2025). The present results support this view by showing that the proposed method reduces unnecessary degrees of freedom and avoids excessive spectral enrichment when the dominant error source lies elsewhere. The use of spectral–finite element methodology is also justified by previous developments in spectral element and hp-adaptive methods. Spectral element methods combine high-order accuracy with element-wise flexibility, while hp-adaptive strategies can respond to both local singularities and smooth regions through mesh refinement and polynomial enrichment (Mavriplis, 1994; Valenciano & Owens, 2000; Xu et al., 2018). In the present method, this philosophy is extended to fractional quantum systems by adding spectral control of the fractional operator itself. This extension is important because fractional operators cannot be treated as ordinary local differential operators, and their approximation requires careful control of both low-frequency and high-frequency components (Hafeez & Krawczuk, 2024; Tang et al., 2020).

The results also highlight the importance of treating variable coefficients and nonlocal kernels explicitly within the error estimator. In heterogeneous media, coefficient approximation errors become significant, especially near sharp transitions or regions with strong spatial variability. Similarly, nonlocal kernel approximation and quadrature errors can influence the solution because the nonlocal term couples distant points in the domain. Operational matrix and nonlocal heat-conduction studies have also emphasized that fractional and nonlocal problems often require specialized algebraic or integral approximation strategies to maintain accuracy and stability (Khalil et al., 2019; Khalil et al., 2021). From a broader computational perspective, the proposed method is positioned within a growing research direction that seeks reliable numerical solvers for fractional, nonlocal, and quantum-inspired PDEs. Recent work on quantum algorithms for PDEs and variational quantum solvers suggests that high-dimensional and nonlinear PDEs eventually benefit from quantum or hybrid computational strategies (Nguyen, 2025; Sarma et al., 2024). However, for current deterministic simulations, certified adaptive spectral– finite element methods remain more directly implementable and mathematically transparent. The present framework could therefore serve as a classical certified benchmark for future comparisons with quantum or hybrid quantum-classical solvers.

8. Conclusion

The CSPA-SFEM method is an innovative numerical algorithm developed for use in solving Non-Local Quantum Fractional PDEs having variable coefficients by overcoming the major difficulties associated with the computation of numerical solutions to equations containing fractional and non-local operators, whilst allowing the use of natural and simple operator function meshes (i.e., two-dimensional polygons, three-dimensional tetrahedra) onto which fractional or non-local differential operators can be implemented. The CSPA-SFEM methodology provides high accuracy for computing solutions of non-local and fractional hyperbolic PDE problems by combining the spatially dense data operator structures that result from solving the original PDE of the problem under both implicit and explicit conditions and solving for the eigenvalues (determined from boundary value problems) of each operator in the spatial direction using a dual-pronged approach: First, by solving for the "coarse" solution in the traditional spectral approximation approach using global and local FEMs; Second, by using adaptive mesh refinement in conjunction with spectral quality enhancement to supply the necessary local mesh refinement needed in the spectral solutions to accurately approximate the discrete solution of the governing system according to the governing spatial, temporal, coefficient, quadrature, and non-local kernel boundaries would be determined; thereby providing a certified mechanism by which to achieve a higher degree of assurance regarding the numerical quality of the computed solution. In support of the application of CSPA-SFEM for accurately calculating the solutions of long-term simulations of quantum-type systems, numerical simulations demonstrate the use of CSPA-SFEM to provide stable and accurate solutions with good accounting for the fundamental physical characteristics of the governing system of equations. In addition, results show that both the average conservation laws (mass and energ y conservation) associated with two-and three-dimensional CSPA-SFEM simulations demonstrate that negligible drift in the mass or Hamiltonian energy arises during the course of the simulations. These results are very important to provide reliable numerical solutions to problems that exhibit long time histories where man-made or artificial dissipation of energy or numerical accumulation of energy at the boundary conditions would lead to unphysical or unreliable problem results.