A second-order-in-time scheme for the von Neumann equation with singular self-interaction and simulation of the onset of instability

1. Problem Statement

The paper addresses the numerical simulation of the Alber equation, a statistical model used in oceanography to describe the evolution of the second moments (autocorrelation) of unidirectional sea states. The equation is a specific case of the von Neumann equation with a singular self-interaction kernel (a delta function).

• Mathematical Formulation: The equation governs the inhomogeneous part $u(x, y, t)$ of the two-space autocorrelation $R(x, y, t) = \Gamma(x-y) + u(x, y, t)$, where $\Gamma$ is a known homogeneous background (autocorrelation function) and $u$ represents the inhomogeneity.
  $$i\partial_t u + p(\Delta_x - \Delta_y)u + q(V(x, t) - V(y, t))(\Gamma(x-y) + u) = 0$$
  where $V(x, t) = u(x, x, t)$ is the trace of the solution.
• Physical Context: The equation exhibits a bifurcation between a stable regime (nonlinear Landau damping) and an unstable regime (modulation instability). The unstable regime is critical for understanding the formation of rogue waves (extreme ocean events).
• Numerical Challenges:
  • Existing literature contains only one first-order-in-time numerical method for this specific singular case.
  • The equation involves a singular nonlinearity (trace term), making standard analysis tools for von Neumann/Wigner equations difficult to apply.
  • There are no global-in-time existence results for the focusing case ($pq > 0$), and solutions can grow by orders of magnitude, requiring schemes that preserve physical invariants to avoid numerical blow-up or instability.

2. Methodology

The authors propose a structure-preserving, linearly implicit, second-order-in-time scheme combined with fourth-order spatial discretization.

A. Time Discretization (Relaxation-Crank-Nicolson)

The scheme introduces an auxiliary variable $\phi(x, y, t) = u(x, x, t) - u(y, y, t)$ to linearize the nonlinear term. The time-stepping is a relaxation-Crank-Nicolson method:

1. Auxiliary Update: $\Phi^{n+1/2}$ is updated via linear extrapolation from previous steps.
2. Main Update: A Crank-Nicolson step is applied to $u$, treating the nonlinear term implicitly via the auxiliary variable $\Phi^{n+1/2}$.
   • This avoids solving a nonlinear system at each timestep (linearly implicit).
   • It preserves the symmetries of the original problem (e.g., $\Phi$ remains real-valued).
3. Initialization: The authors highlight that a "naive" first-order initialization for the auxiliary variable leads to a loss of temporal order for $\phi$. They propose an advanced second-order initialization to ensure global second-order convergence for both variables.

B. Spatial Discretization

• Method: Fourth-order central finite differences on a uniform grid with periodic boundary conditions.
• Implementation: The hyperbolic Laplacian $(\Delta_x - \Delta_y)$ is implemented using a sparse matrix (Kronecker product), and the pointwise multiplication is handled via diagonal matrices.
• Periodization: Special care is taken to periodize the background autocorrelation $\Gamma$ correctly to handle boundary interactions in the finite domain.

C. Theoretical Properties

The scheme is proven to satisfy:

• Discrete Balance Law: A discrete version of the $L^2$ growth law, ensuring the numerical solution respects the energy exchange between the inhomogeneity and the background.
• Invariant Preservation: The scheme exactly preserves discrete versions of two key invariants ($I_0$ and $I_1$) and approximately preserves others ($I_2, I_3$), which is crucial for long-time stability in unstable regimes.
• Consistency & Boundedness: The scheme is shown to be consistent with order $O(\tau^2 + h^4)$ and unconditionally $L^2$-bounded.

3. Key Contributions

1. New Numerical Scheme: The first second-order-in-time, structure-preserving scheme for the Alber equation with singular self-interaction.
2. Initialization Analysis: Identification and correction of order reduction in the auxiliary variable via a specific second-order initialization strategy.
3. Invariant Preservation: Demonstration that the scheme preserves discrete invariants, allowing for reliable simulation of solutions that grow by factors of 100 or more.
4. Discovery of Meta-Stable Patterns: Numerical discovery of a "space-time lattice" of coherent structures in the fully nonlinear unstable regime.
5. Monte Carlo Investigation: A systematic study revealing that the maximum amplitude of rogue events depends primarily on the background spectrum intensity, not the specific initial inhomogeneity.

4. Key Results

A. Validation and Convergence

• Using an exact soliton solution, the authors verified that the advanced initialization yields second-order convergence in time for both $u$ and $\phi$.
• The naive initialization yields second-order convergence for $u$ but only first-order for $\phi$.
• The error analysis shows that the bulk of the numerical error is due to a phase shift (phase velocity error), while the shape of the soliton is preserved with high accuracy.

B. Stability and Instability Dynamics

• Stable Regime: Confirmed nonlinear Landau damping where inhomogeneities disperse.
• Unstable Regime (Modulation Instability):
  • Linear vs. Nonlinear: Linear stability analysis accurately predicts the initial exponential growth rate but fails to predict the maximum amplitude or the formation of coherent structures.
  • Meta-Stable Pattern: In the fully nonlinear regime, the solution evolves into a "cone-supported lattice" of localized structures. This pattern is robust and its scaling (peak amplitude, spacing) depends on the background $\Gamma$, not the initial condition $u_0$.
  • Growth: Solutions can grow from small initial perturbations ($\sim 0.05$) to large amplitudes ($\sim 6$), representing a two-order-of-magnitude increase.

C. Amplification Factors and Rogue Waves

The authors introduce two metrics:

• Inhomogeneity Amplification Factor (IAF): Ratio of max $|u|$ to initial max $|u|$.
• Total Amplification Factor (TAF): Ratio of max $|u + \Gamma|$ to max $|\Gamma|$. This measures the physical extreme event relative to the background.

Monte Carlo Findings:

• Dependence on Background: The TAF depends quadratically on the background intensity parameter $C$.
• Thresholds:
  • Instability onset occurs around $C \approx 1$.
  • For $C$ slightly above the threshold, IAF is large, but TAF remains close to 1 (no significant rogue waves).
  • For $C \approx 1.4$, TAF $\approx 1.4$ (40% increase in variance).
  • Rogue wave conditions ($TAF > 2.84$) require very strong instability ($C \approx 2$), which is rare in observed sea states.
• Conclusion: Modulation instability creates recurrent hotspots of increased variance. While these do not guarantee a single rogue wave, they make extreme events significantly more likely than predicted by linear statistics.

5. Significance

• Oceanography: Provides a robust tool for simulating the onset of rogue waves. The findings suggest that while linear instability is necessary, the intensity of the background sea state is the primary driver for the formation of extreme events.
• Numerical Analysis: Demonstrates the necessity of structure-preserving schemes (invariant preservation) for singular nonlinear PDEs where solutions can grow unboundedly. It resolves the order-reduction issue common in relaxation schemes for this class of problems.
• Physical Insight: Challenges the assumption that linear instability alone predicts extreme events. It highlights the role of nonlinear saturation and the formation of meta-stable coherent structures in determining the actual risk of rogue waves.

The code for the solver and Monte Carlo analysis is made available via a GitHub repository, facilitating further research in this domain.