arXiv:2603.23185v2 Announce Type: replace Abstract: We develop arbitrarily high-order, stationarity-preserving stabilized finite element methods for multidimensional nonlinear hyperbolic balance laws on Cartesian grids. We aim at approximating all the steady states of the problem at hand, including non-trivial genuinely multidimensional equilibria, with a level of accuracy higher than the nominal one of the underlying scheme. We formalize more precisely the meaning of stationarity preservation, providing some technical conditions for its realizability. We then recast the multidimensional global-flux quadrature of Barsukow et al. (Num. Meth. PDEs, 2025) as a local preprocessing of the physical fluxes that maps continuous polynomial vector fields to a local space with Raviart--Thomas-type structure. Both the Galerkin and SUPG formulations are recast in this setting. The resulting methods extend the stationarity-preserving finite-volume approach of Barsukow et al. (J. Comput. Phys., 2026) to high-order continuous finite elements and Barsukow et al. (Num. Meth. PDEs, 2025) to nonlinear balance laws. We analyze key properties of the proposed schemes, including local conservation and nodal superconvergence of the discrete steady kernel, and we discuss their relation to low-Mach-compliant discretizations. We apply the framework to the compressible Euler equations with gravity. A simple source-term reformulation yields machine-precision preservation of isothermal hydrostatic equilibria. Extensive numerical benchmarks, including moving equilibrium, near-equilibrium, and instability-dominated regimes, demonstrate clear improvements in robustness and accuracy over standard SUPG and reference finite-volume methods.