Jet Marching Method: High-Order PDE Solver

• Jet Marching Method is a high-order numerical solver for Hamilton–Jacobi PDEs, explicitly propagating both function values and their derivatives (jets) at each mesh point.
• It uses compact Hermite interpolants and Dijkstra-style update schemes to achieve high-order convergence for both function and gradient accuracy.
• Applications range from accurate eikonal and quasipotential computations to precise level-set reinitialization and jet acoustics, demonstrating robust efficiency.

The Jet Marching Method (JMM) encompasses a family of high-order numerical solvers for Hamilton–Jacobi-type partial differential equations (PDEs), with leading applications in solving the eikonal equation, quasipotential computation in stochastic differential equations, accurate level-set reinitialization, and parabolic stability analyses in jet acoustics. JMMs are characterized by the explicit evolution (“marching”) of not only the field value but also its derivatives—the so-called jet—at each mesh point, employing Hermite or related interpolations and Dijkstra-like label-setting algorithms for optimal efficiency and accuracy. These methods yield higher-order convergence for both function and gradient, and in certain formulations, systematically extend to second derivatives, facilitating computation of quantities such as WKB amplitudes and geometric spreading factors (Potter et al., 2020, Paskal et al., 2021, Salac, 2011, Padilla-Montero et al., 7 Aug 2025).

1. Mathematical Formulation and Jet Structure

At the core of JMMs is the explicit representation and propagation of jets: for a scalar field $\tau(x)$ (e.g., solution of the eikonal equation), the total first jet at $x$ is

$j^1\tau(x) = \{\tau(x), \nabla\tau(x)\}.$

This stands in contrast to classical methods which propagate only $\tau(x)$. For higher accuracy and applications such as WKB analysis and level-set methods, second jets—including the Hessian—are also evolved:

$$j^2\tau(x) = \{\tau(x), \nabla\tau(x), \nabla\nabla\tau(x)\}.$$

On a discrete mesh (structured or unstructured), each node $x_i$ thus stores, e.g., $T_i \approx \tau(x_i)$, $G_i \approx \nabla\tau(x_i)$, and, if needed, second derivatives. These nodal jets enable high-order Hermite interpolation over cells, essential for the method's accuracy and compactness (Potter et al., 2020, Salac, 2011).

2. Semi-Lagrangian Update Schemes and Interpolant Construction

JMMs employ compact high-order semi-Lagrangian update stencils, leveraging Hermite-type piecewise-polynomial interpolants on the local mesh neighborhood. In the eikonal and quasipotential settings, the update at a target node $\hat{y}$ proceeds by:

• Constructing a local Hermite interpolant $\mathcal{T}$ for $\tau$ (or $V$ in quasipotential problems) over a simplex or facet spanned by neighboring Valid nodes, fully exploiting the available jets.
• Parameterizing characteristics (e.g., trial rays, minimum action paths) via high-degree polynomials (cubic for Hermite), with endpoint constraints on position and tangent, ensuring locally optimal path representation.
• Minimizing an action or Fermat-type cost functional—often using Simpson's rule or higher quadrature—over the interpolated characteristic, constraining to admissible base points and tangent directions.
• Setting the candidate value and gradient at $\hat{y}$ to the action minimum and the tangent of the minimizing characteristic, i.e., $T(\hat{y}) \leftarrow F_\text{min}$, $$\nabla T(\hat{y}) \leftarrow s(\hat{y}) \gamma'(L)/\|\gamma'(L)\|$$ (Potter et al., 2020, Paskal et al., 2021).

These local updates, involving small-scale nonlinear minimizations (often efficiently solved via Newton-type iterations), are crucial for ensuring both accuracy and computational tractability.

3. Algorithmic Structure: Label-Setting, Heap Ordering, and Causality

JMMs are fundamentally Dijkstra-style, label-setting solvers. Each mesh node is assigned a state (Far, Trial, Valid) and a value (function and derivatives), with a priority queue (typically a min-heap) of Trial nodes keyed on their tentative values. The main loop is:

1. Initialize mesh (e.g., Dirichlet data on boundary $\Gamma_h$); set initial jets (analytically or via interpolation).
2. While the heap is non-empty:
   • Pop the node with minimum tentative value; mark as Valid.
   • For each neighbor stencil containing this node, and for all base subsets:
     • If the base is Valid, perform the local semi-Lagrangian update.
     • If the new candidate value is superior (e.g., smaller for arrival times), update Trial status and heap key.
3. For methods supporting higher derivative marching, a secondary sweep updates Hessians or related second-jet quantities (Potter et al., 2020, Salac, 2011).

Causality and monotonicity are carefully enforced, especially at points with limited accepted neighbors, to prevent information corruption and spurious oscillations.

4. Accuracy, Convergence, and Special Cases

Quantitative error results delineate the advantage of JMMs. The Hermite interpolants constructed from jets at $O(h)$-separated nodes are $O(h^{4-p})$ in the $p$th derivative ($p \leq 1$), and with fifth-order local cost approximations, the per-update error is $O(h^5)$. Globally, after accounting for graph depth ($O(1/h)$) and dynamic programming error propagation, the convergence rate is:

Quantity	Local Error	Global Error	Special Case (linear/constant $s$)
Function (T)	$O(h^5)$	$O(h^3)$	$O(h^3)$ (can be $O(h^4)$ locally)
Gradient ($\nabla T$)	$O(h^4)$	$O(h^2)$	$O(h^3)$

In cases where the slowness function $s(x)$ (eikonal) or underlying field is linear or constant, the local cubic interpolation precisely recovers the true characteristic, eliminating the typical one-order loss, and resulting in empirically cubic convergence for both field and gradient (Potter et al., 2020). These properties extend to the computation of minimum action paths in anisotropic quasipotential settings, with cubic-Hermite updates yielding global $O(h^2)$ accuracy for the quasipotential and $O(h^{1.8})$ in its gradient (Paskal et al., 2021).

5. Extension to Second Derivatives and WKB Amplitude Computation

After a pass marching the 1-jet, JMMs can reconstruct second derivatives (Hessian, or in application-specific notation, $T_{ij}$) via upwind “cell-marching” procedures. For example, central differences of gradient components on cell edges yield midpoint Hessians with $O(h^2)$ accuracy; bilinear extrapolation and averaging enforce $C^1$ continuity at corners.

This capability is central to WKB approximations of high-frequency wave propagation, where the amplitude $A$ in $u \approx A e^{i \omega \tau}$ evolves according to geometric spreading $J$ computed using second derivatives:

$$J(\hat{y}) = |1 + \epsilon \cdot (\Delta T(x) - t^T \nabla s(x))| J(x), \quad \epsilon \approx \int_0^L c(\gamma(s)) ds \approx L\,c\left(\frac{x+\hat{y}}{2}\right).$$

The amplitude formula is

$$A(\hat{y}) = \frac{e^{i \pi/4}}{2 \sqrt{2\pi\, \omega}} \sqrt{\frac{c(\hat{y})}{J(\hat{y})}},$$

allowing a complete solver for the WKB system alongside the phase (Potter et al., 2020).

6. Implementation Aspects and Computational Efficiency

JMMs exhibit $O(N \log N)$ complexity for $N$ mesh nodes. Each is “settled” once; at each step, a small number of neighbor stencils is tested. The nonlinear minimizations per update converge rapidly (in a few Newton steps). Heap operations typically account for a minor fraction of computational time.

• Structured and unstructured meshes are supported, with the mesh represented via adjacency lists.
• Dirichlet boundary conditions are initialized by direct assignment or Hermite interpolation; special treatments (e.g., rarefaction fans, specular reflection) accommodate physical singularities.
• For higher-dimensional or highly anisotropic problems (notably in efficient jet marching for quasipotential SDEs), tailored anisotropic stencils based on angular bins substantially reduce neighbor counts (O(1) per node) without loss of accuracy, following the methodology of Mirebeau (Paskal et al., 2021).
• Storage per node includes field value, gradient vector, and optional Hessian or auxiliary fields.
• Reproducible codebases exist in MATLAB/C++ (see (Potter et al., 2020)).

7. Applications and Representative Results

JMMs have been successfully applied in varied domains:

• Eikonal equation: Achieving third-order convergence for arrival time and gradient, with sharp shockline localization and accurate WKB prefactor computation, as demonstrated on constant, linear, and nonlinear slowness fields (Potter et al., 2020).
• Quasipotential computation for SDEs: Enabling precise estimation of invariant measures, mean escape times, and most probable exit paths in metastable systems; critical for applications in stochastic dynamics (Paskal et al., 2021).
• Level-set reinitialization: The Augmented Fast Marching Method (AFMM) achieves second-order accurate distance functions and gradients, and first-order curvature, vastly outperforming classical FMM, and preserving unconditional stability and smoothness (Salac, 2011).
• Jet acoustics and wavepacket modeling: The Plane-Marching Parabolized Stability Equations (PM–PSE), a jet-marching method in this context, compute spatial evolution of wavepacket modes in complex jets with direct application to prediction and analysis of mixing noise and instability structures in supersonic twin jets (Padilla-Montero et al., 7 Aug 2025).

A summary of representative convergence results:

Problem/Scenario	Quantity	Observed Order	Reference
Eikonal, $s \equiv 1$	$T$, $\nabla T$	$\sim$3	(Potter et al., 2020)
Quasipotential, 2D SDE	$V$, $\nabla V$	2, $\sim$1.8	(Paskal et al., 2021)
Level-set, 2D/3D geometry	$\phi$, $\nabla\phi$, $\kappa$	2, 2, 1	(Salac, 2011)
Jet acoustics (PM–PSE)	Wavepacket modes	High-fidelity structure, agreement with PIV	(Padilla-Montero et al., 7 Aug 2025)

In all cases, JMMs offer a robust, high-order, and efficient alternative to classical single-field update schemes, with a particularly significant impact where accurate gradient and curvature information are essential.

---

References:

• "Jet Marching Methods for Solving the Eikonal Equation" (Potter et al., 2020)
• "An efficient jet marcher for computing the quasipotential for 2D SDEs" (Paskal et al., 2021)
• "The Augmented Fast Marching Method for Level Set Reinitialization" (Salac, 2011)
• "Plane-marching PSE wavepacket models for perfectly-expanded twin jets" (Padilla-Montero et al., 7 Aug 2025)