Tracking the boundary between absolute/convective instability using adjoint equations

Technical Summary: Tracking the Boundary Between Absolute/Convective Instability Using Adjoint Equations

Problem Statement
Determining the boundary between absolute instability (AI) and convective instability (CI) in open flows is a fundamental challenge in spatio-temporal stability theory. Conventionally, this boundary is defined by a stationary saddle point in the complex wavenumber plane where the group velocity vanishes ($d\omega/dk = 0$) and the absolute growth rate is zero ($\text{Im}(\omega) = 0$).
The standard numerical workflow involves a nested, computationally expensive process:

1. Solving temporal or spatial eigenvalue problems at numerous complex wavenumbers.
2. Associating eigenvalues with specific modal branches (a delicate task for non-normal problems where branches approach one another).
3. Approximating derivatives via finite differences to locate saddle points.
4. Embedding these saddle searches within an outer loop over physical parameters (e.g., Reynolds number) to find the zero-growth contour.
   This approach becomes costly for large-scale, multi-field discretizations and is sensitive to errors in modal branch association. Furthermore, standard continuation methods in a single physical parameter often fail when the solution manifold contains folds (turning points), requiring manual parameter switching.

Methodology
The authors propose a direct continuation method that treats the neutral saddle boundary as a one-dimensional solution manifold in a combined state–parameter space. The core of the methodology involves:

1. Adjoint Formulation of the Stationary Condition:
   Instead of differentiating eigenvalues numerically, the zero-group-velocity condition ($d\omega/dk = 0$) is expressed as an adjoint solvability residual. For a generalized eigenvalue problem $L(k, \omega; R)q = [A(k, R) - \omega B(k, R)]q = 0$, the condition is derived as:
   $$S \equiv (q^\dagger)^H C q = 0, \quad \text{where } C = \frac{\partial A}{\partial k} - \omega \frac{\partial B}{\partial k}$$
   Here, $q$ and $q^\dagger$ are the direct and adjoint eigenvectors, respectively.

2. Augmented Nonlinear System:
   The direct eigenproblem, adjoint eigenproblem, the adjoint saddle residual ($S=0$), and the neutral growth condition ($\text{Im}(\omega)=0$) are solved simultaneously. To handle the scale and phase freedoms of the eigenvectors, two complex gauge conditions are imposed (e.g., biorthogonality $(q^\dagger)^H B q = 1$ and a reference vector constraint). This creates a square Newton system for fixed parameters.

3. Scaled Pseudo-Arclength Continuation:
   To trace the boundary across two physical parameters ($R_1, R_2$), the fixed-parameter constraint is removed, resulting in an underdetermined system. The authors employ a scaled pseudo-arclength continuation method.

   • Scaling: A diagonal scaling matrix $S_x$ is constructed from the initial state to normalize the magnitudes of eigenvectors, eigenvalues, and parameters, ensuring the tangent vector is well-conditioned.
   • Predictor-Corrector: The method uses a scaled tangent vector to predict the next point and a damped Newton iteration to correct it.
   • Robustness: This formulation allows the algorithm to cross "folds" (turning points) in the projection of the solution manifold onto physical parameter space without switching the continuation variable, a limitation of natural continuation.

4. Jacobian Assembly:
   The Newton Jacobian is assembled semi-analytically. The large blocks corresponding to eigenvector variations are derived exactly from the linearized operators, while only the scalar columns (frequency, wavenumber, and control parameters) are approximated via centered finite differences.

Key Results and Validation
The method was validated and applied across three distinct levels of complexity:

1. Linear Ginzburg–Landau Equation (Scalar Limit):

   • The method was tested against an analytical neutral-saddle boundary ($\mu_c = U^2/4D$).
   • Result: The numerical continuation recovered the analytical boundary with machine-precision accuracy (errors $\sim 10^{-16}$), verifying that the adjoint formulation correctly reduces to the exact group-velocity condition in the scalar limit.

2. Planar Gaussian Wake (Non-Normal Matrix Problem):

   • Applied to an Orr–Sommerfeld problem with a mapped collocation discretization.
   • Cross-Validation: The continued boundary was compared against 31 independently formulated finite-difference saddle corrections. The relative difference in critical Reynolds numbers was approximately $10^{-8}$.
   • Efficiency: A comparison with a nested complex-wavenumber and parameter-plane scanning method showed that the direct adjoint continuation was 14.0 to 30.6 times faster than the scanning approach. The cost of the scanning method increased rapidly with grid refinement, whereas the continuation cost remained focused on the boundary itself.
   • Topology: Spatial root tracking confirmed that the computed stationary points correspond to the physical Briggs–Bers pinch points.

3. Viscoelastic Oldroyd–B Free-Surface Film (Multi-Field Coupled Problem):

   • Applied to a coupled eigenproblem involving hydrodynamic, constitutive, and interfacial variables.
   • Fold Crossing: The method successfully traced the neutral-saddle manifold through a genuine turning point (fold) in the Reynolds number projection. Natural continuation in Reynolds number failed at this point, whereas the pseudo-arclength method crossed the fold and revealed the post-fold branch.
   • Physical Insight: The results revealed a re-entrant CI–AI–CI boundary geometry. As the Deborah number ($De$) increased, the absolutely unstable region first expanded (AI region grew) and then contracted, demonstrating that elasticity does not act as a uniformly stabilizing or destabilizing parameter in this spatio-temporal context.

Significance and Claims
The paper claims that adjoint-augmented pseudo-arclength continuation offers a direct and computationally efficient alternative to the conventional nested saddle-search and parameter-plane reconstruction workflow.

• Computational Efficiency: By advancing directly along the one-dimensional neutral-saddle manifold, the method avoids the repeated spectral solves required by scanning methods. The computational cost is shown to be significantly lower, particularly as the required resolution of the boundary increases.
• Numerical Robustness: The formulation handles the "branch association" problem inherent in non-normal eigenvalue problems by solving the stationarity condition algebraically within the system, rather than relying on finite-difference derivatives of tracked branches.
• Geometric Capability: The use of scaled pseudo-arclength continuation allows for the automatic crossing of folds in the physical parameter space, a capability essential for capturing complex boundary geometries (such as re-entrant transitions) that are missed by single-parameter continuation.
• Generality: The framework is presented as a general tool for frequency-affine generalized eigenvalue problems, demonstrated here to work effectively on scalar, matrix, and multi-field coupled systems without requiring problem-specific reformulations of the continuation strategy.

The authors conclude that while the method assumes a locally simple finite eigenvalue and requires an initial saddle estimate, it provides a robust framework for computing connected neutral-saddle manifolds in multi-parameter spatio-temporal stability problems.