Multi-domain spectral approach for Zakharov-Kuznetsov equations in 3D with cylindrical symmetry

1. Problem Statement

The paper addresses the numerical simulation of the three-dimensional generalized Zakharov-Kuznetsov (gZK) equation:
$$u_t + (\Delta u + u^p)_x = 0, \quad (x, y, z) \in \mathbb{R}^3, t \in \mathbb{R}$$
where $p > 1$. The specific focus is on the $L^2$-critical case where the nonlinearity power is $p = 7/3$.

Key Challenges:

• Dimensionality and Symmetry: The equation models weakly magnetized ion-acoustic waves in plasma. While the full problem is 3D, solitary waves travel along a preferred axis ($x$) with cylindrical symmetry in the transverse plane ($y, z$). This reduces the problem to 2D in coordinates $(x, \rho)$, where $\rho = \sqrt{y^2+z^2}$.
• Fractional Nonlinearity: The critical power $p=7/3$ is fractional. This introduces significant numerical difficulties, particularly regarding the definition of $u^{7/3}$ for negative values and the preservation of real-valued solutions.
• Singular Terms: In cylindrical coordinates, the Laplacian introduces a singular term $\frac{1}{\rho}u_\rho$, which becomes undefined at the symmetry axis ($\rho=0$). Standard spectral methods (like Fourier) fail or become inefficient near this singularity.
• Blow-up Dynamics: The $L^2$-critical regime is the threshold between global existence (dispersion) and finite-time blow-up (wave collapse). Accurately tracking the solution profile as it approaches singularity requires extremely high resolution and stability.

2. Methodology

The authors propose a novel Multi-Domain Spectral Approach combined with an implicit time-stepping scheme.

A. Spatial Discretization

The computational domain is split based on the behavior of the solution and the coordinate singularity:

1. Fourier Spectral Method (x-direction):

   • The $x$-dependence is treated on a large torus $T_{2\pi L}$ assuming periodicity for localized solutions.
   • A Fast Fourier Transform (FFT) is used for spatial derivatives.

2. Multi-Domain Chebyshev Collocation (ρ-direction):
   To handle the singularity at $\rho=0$ and the need for high resolution near the core, the radial domain is partitioned into two sub-domains with different independent variables:

   • Domain I ($\rho < \rho_0$, typically $\rho_0=1$):
     • Uses the transformation $s = \rho^2$. This removes the singularity because functions regular at the axis are analytic in $\rho^2$.
     • The equation is rewritten in terms of $s$, yielding a less singular form: $u_t + (u_{xx} + 4s u_{ss} + 4u_s + u^p)_x = 0$.
     • Discretized using Chebyshev polynomials in the transformed variable $l \in [-1, 1]$.
   • Domain II ($\rho_0 < \rho < \rho_1$):
     • Uses the standard variable $\rho$. Since $\rho$ is bounded away from zero, the $\frac{1}{\rho}$ term is non-singular.
     • Discretized using Chebyshev polynomials.
   • Matching Conditions: The solutions in Domain I ($u_I$) and Domain II ($u_{II}$) are coupled at $\rho = \rho_0$ by enforcing $C^1$ continuity (continuity of the function and its first derivative). This is implemented via a $\tau$-method, replacing specific rows in the differentiation matrices with the boundary conditions.

B. Time Integration

• Implicit Runge-Kutta (IRK4): Due to the stiffness introduced by the third-order derivative in $x$ and the spectral discretization in $\rho$, explicit schemes are inefficient. The authors use a 4th-order, 2-stage implicit Gauss-Legendre Runge-Kutta method (Hammer-Hollingsworth).
• Nonlinear Solver: The implicit step results in a nonlinear system solved via a simplified Newton iteration. The linear systems are solved by inverting block matrices derived from the differentiation operators.
• Handling Fractional Powers: To manage the fractional nonlinearity $u^{7/3}$ and avoid complex numbers arising from FFT artifacts, the authors enforce real-valuedness after every inverse FFT and compute powers as $|u|^{7/3}\text{sign}(u)$.

3. Key Contributions

1. Novel Numerical Framework: Development of a hybrid spectral method (Fourier in $x$, Multi-domain Chebyshev in $\rho$) specifically designed for 3D dispersive equations with cylindrical symmetry and singular terms.
2. Efficient Handling of Critical Fractional Nonlinearity: Successful implementation of a stable scheme for the $p=7/3$ case, overcoming numerical instabilities typically associated with fractional powers and singular coordinates.
3. High-Resolution Blow-up Tracking: The method achieves spectral accuracy (exponential convergence) and maintains $C^1$ continuity across domain boundaries, allowing for precise tracking of the solution profile up to the point of singularity.
4. Threshold Determination: The framework provides a robust tool to numerically verify the theoretical threshold for global existence vs. blow-up in the $L^2$-critical regime.

4. Key Results

The authors applied their method to the 3D critical ZK equation ($p=7/3$) and obtained the following results:

• Ground State Construction: They numerically constructed the ground state soliton $Q$ (radial positive solution).

  • Calculated Mass: $\|Q\|_{L^2} \approx 7.98$ (Mass squared $\approx 63.78$).
  • Calculated Energy: $E[Q] \approx 0$ (consistent with theoretical predictions for the critical case).
  • The solution was verified to be stable under time evolution with errors on the order of $10^{-10}$.

• Stability Threshold Verification:

  • Sub-critical ($\|u_0\|_{L^2} < \|Q\|_{L^2}$): Solutions with initial mass slightly below the ground state (e.g., $\lambda = 0.99$ for $u = \lambda Q$) were observed to disperse and flatten over time.
  • Super-critical ($\|u_0\|_{L^2} > \|Q\|_{L^2}$): Solutions with initial mass slightly above the ground state (e.g., $\lambda = 1.01$) showed initial growth. While finite domain effects caused temporary stabilization for very small perturbations, larger perturbations ($\lambda = 1.1$) led to finite-time blow-up in the $L^\infty$ norm.
  • Gaussian Data: Similar thresholds were observed for Gaussian initial data. Mass below $\approx 5.7$ (for specific width parameters) led to dispersion; mass above led to blow-up.

• Blow-up Dynamics:

  • The simulations successfully tracked the formation of the singularity.
  • Analysis of the radiation pattern during blow-up revealed a cone of radiation propagating in the negative $x$-direction with an opening angle of approximately $60^\circ$ (half-angle $\approx 29.7^\circ$), consistent with previous theoretical and 2D results.

5. Significance

• Theoretical Validation: The study provides strong numerical evidence that the ground state mass serves as the sharp threshold for global existence versus finite-time blow-up in the 3D $L^2$-critical ZK equation, a regime that was previously under-explored numerically due to computational difficulties.
• Physical Relevance: The findings are crucial for understanding wave collapse and energy localization in magnetized plasmas (laboratory, space, and astrophysical contexts), where the balance between nonlinearity and dispersion is critical.
• Methodological Advancement: The multi-domain spectral approach presented here offers a template for solving other high-dimensional nonlinear dispersive equations with symmetries and singularities, particularly those involving fractional powers. It bridges the gap between analytical theory and high-fidelity simulation in complex physical regimes.