Analytical and Gegenbauer spectral solutions of damped wave equations with singular potentials in ℝn

Technical Summary: Analytical and Gegenbauer Spectral Solutions of Damped Wave Equations with Singular Potentials in $\mathbb{R}^n$

Problem Formulation
This study investigates the existence, uniqueness, and numerical approximation of solutions to a multi-dimensional damped wave equation defined on a rectangular domain $\Omega = [0, \pi]^n$. The governing equation is a mixed hyperbolic-parabolic operator incorporating a singular source term of the form $|x|^{-\alpha}$:
$$\frac{\partial^2 u}{\partial t^2} + \gamma \frac{\partial u}{\partial t} - \Delta u + \frac{u}{|x|^\alpha} = f(x, t)$$
The problem is formulated under Cauchy initial conditions ($u(x,0)=\phi(x)$, $u_t(x,0)=\psi(x)$) and non-homogeneous Neumann boundary conditions ($\partial u / \partial \nu = g_j$). The analysis focuses on the regime where the singularity parameter satisfies $0 < \alpha < n$ and the damping coefficient $\gamma > 0$. The presence of the singular potential introduces significant analytical challenges regarding integrability and regularity, while the mixed nature of the operator requires addressing both wave propagation and diffusive dissipation simultaneously.

Methodology
The paper employs a dual analytical and numerical approach:

1. Analytical Framework (Eigenfunction Expansion):

   • The solution is decomposed into a transient homogeneous part ($u_h$) and a steady-state particular part ($u_p$).
   • The particular solution is constructed using a Fourier series in time, reducing the time-dependent PDE to a system of coupled elliptic equations for each temporal harmonic $k$.
   • By introducing complex-valued spatial functions $U_k(x) = A_k(x) - iB_k(x)$, the system is reduced to a complex Helmholtz-type equation: $(L - \lambda_k)U_k = F_k$, where $L = -\Delta + |x|^{-\alpha}$ and $\lambda_k = k^2 - i\gamma k$.
   • The solution is explicitly represented as a double series using the eigenfunctions $\Phi_m(x) = \prod \cos(m_j x_j)$ of the Neumann Laplacian. The coefficients are determined via projection onto the orthogonal basis, accounting for both internal sources and boundary fluxes.

2. Numerical Framework (Gegenbauer Spectral Collocation):

   • A high-order spectral method is developed using Gegenbauer polynomials with parameter $\lambda=1$ (which reduces to Chebyshev polynomials for optimal node distribution).
   • Spatial discretization utilizes Chebyshev-Gauss-Lobatto (CGL) nodes on a tensor-product grid to handle the domain geometry.
   • The singular potential $|x|^{-\alpha}$ is regularized near the origin using $V_\epsilon(x) = (|x|^2 + \epsilon^2)^{-\alpha/2}$ to prevent numerical overflow while preserving asymptotic behavior.
   • Time integration is performed using a second-order Implicit-Explicit (IMEX) scheme, treating the stiff diffusion term implicitly and the hyperbolic terms explicitly to ensure stability without restrictive CFL conditions.

Key Contributions and Theoretical Results
The paper establishes rigorous mathematical properties of the solution:

• Convergence of Coefficients: It is proven that the expansion coefficients decay asymptotically as $O(|m|^{-n-\alpha})$, ensuring the series is well-defined despite the singularity, provided $0 < \alpha < n$.
• Stability and Non-Resonance: The damping coefficient $\gamma > 0$ is shown to be critical. It ensures the complex spectral denominator $|\mu_m - \lambda_k| \geq \gamma k > 0$, preventing division by zero and eliminating pure resonance phenomena.
• Regularity and Sobolev Spaces:
  • The solution $u$ is shown to be bounded and well-posed in $L^2(\Omega)$.
  • The spatial gradient $\nabla u$ remains square-integrable (i.e., $u \in H^1(\Omega)$) only under the stricter condition $0 < \alpha < n-2$. For $n-2 \leq \alpha < n$, the gradient exhibits a singularity at the origin.
• Global Stability: An explicit norm estimate is derived showing that the $L^2$-norm of the solution decays exponentially as $e^{-\gamma t/2}$, confirming global stability and energy dissipation.
• Error Analysis: Theoretical error bounds for the Gegenbauer collocation method are established. The spatial discretization error converges algebraically at a rate of $O(N^{-(n-\alpha-1/2)})$, reflecting the regularity loss induced by the singularity. For weak singularities, near-spectral convergence is observed.

Numerical Results
Comprehensive numerical experiments were conducted in one, two, and three dimensions.

• Validation: The simulations validate the theoretical convergence rates. For weak singularities ($\alpha = 0.3$), the method achieves high-order convergence (rates $> 3.5$). For strong singularities near the critical threshold ($\alpha \approx n - 0.5$), the convergence rate degrades gracefully to the theoretical algebraic limit.
• Energy Behavior: Numerical plots confirm the exponential decay of energy over time, consistent with the damping parameter $\gamma$.
• Dimensional Scaling: The results demonstrate that the method remains robust and stable across dimensions $n=1, 2, 3$, with computational costs scaling as $O(N^n)$.

Significance
The paper claims to bridge a gap in the literature by providing an explicit analytical representation for $n$-dimensional mixed hyperbolic-parabolic operators under singular potentials, a problem where standard existence theorems often fall short of providing constructive solutions. The work demonstrates that the damping coefficient acts as a fundamental stabilizer against resonance in the presence of singularities. Furthermore, the integrated analytical-numerical approach offers a robust model for characterizing wave propagation and diffusion in non-homogeneous media with localized singularities. The authors position these findings as directly applicable to physical fields such as thermoelasticity, fluid-structure interaction, and electromagnetic wave propagation in dissipative media.