1 Introduction↩︎

Relativistic hydrodynamics (RHD) provides a fundamental description of fluid motion at velocities approaching the speed of light and plays a central role in astrophysics, high-energy physics, and cosmology.  The governing equations feature strong nonlinearities and, in many practically relevant regimes, extremely large Lorentz factors; as a result, analytical solutions are rare and high-fidelity numerical simulation becomes indispensable for elucidating the dynamics of relativistic flows.  Compared with the nonrelativistic Euler equations, RHD exhibits a more intricate and nonlinearly coupled relationship between conservative and primitive variables, which poses additional challenges for designing robust, accurate, and efficient numerical methods.

Early numerical investigations date back to [1], [2], where finite-difference discretizations augmented by artificial viscosity were employed in Lagrangian or Eulerian coordinates. Over the past decades, a wide range of high-resolution and high-order methods have been developed for RHD, including finite-volume, finite-difference, and discontinuous Galerkin (DG) approaches; see, e.g., [3]–[7]. To further enhance the resolution of discontinuities and multiscale structures, adaptive mesh refinement and adaptive moving-mesh techniques have also been incorporated into RHD simulations [8], [9].

A central difficulty in RHD computations is the strict enforcement of intrinsic physical constraints: positivity of the rest-mass density and pressure, and the subluminal constraint on the fluid velocity (i.e., \(|{\boldsymbol{v}}|<1\) with the light speed \(c=1\)).  These constraints define the admissible state set of the system.  Numerical violations (negative density/pressure or superluminal velocity) are not only nonphysical but may also cause immediate breakdown of the computation, particularly in extreme regimes characterized by large Lorentz factors, near-vacuum conditions, or strong shocks.  In practice, such failures are sometimes handled in an ad hoc manner by restarting with more diffusive schemes and/or smaller CFL numbers until admissibility is recovered [8], [10].  This lack of rigor and predictability motivates the development of high-order methods that guarantee admissibility by construction, commonly referred to as physical-constraints-preserving (PCP) or, more generally, invariant-domain-preserving (IDP) schemes.

In the broader setting of hyperbolic conservation laws, the design of provably PCP/IDP schemes has advanced substantially, with two major limiter paradigms emerging.  The first class consists of scaling-type limiters, pioneered by Zhang and Shu, and successfully applied to scalar conservation laws [11], the nonrelativistic Euler equations [12], [13], and the compressible Navier–Stokes equations [14].  The second class enforces admissibility through flux correction, i.e., nonlinear blending of high-order and low-order numerical fluxes.  Representative examples include parametrized maximum-principle/positivity-preserving flux limiters for scalar laws [15], [16], convection-dominated diffusion [17], and the compressible Euler equations [18].  In addition, Hu, Adams, and Shu [19] proposed a sufficient positivity-preserving flux limiting condition for the Euler equations via a convex decomposition of cell averages, an approach later generalized to special RHD by Wu and Tang [5].

The flux-correction philosophy can be traced back to the Flux-Corrected Transport (FCT) method of Boris and Book [20]–[22].  Zalesak extended FCT to fully multidimensional settings [23] and developed a structured-grid limiter for scalar conservation laws [24], laying foundations for later IDP/PCP developments.  Recent work has further clarified the close connection between Zalesak-type FCT limiters and parametrized flux limiters; see, e.g., [15], [25].  Kuzmin and collaborators systematically extended FCT to multidimensional and high-order discretizations [26], [27] and proposed monolithic convex limiting techniques for hyperbolic systems [28].  Other notable FCT-type variants include the Point-Average-Moment PolynomiAl-Interpreted (PAMPA) schemes [29]–[31].  Another influential line is the convex limiting framework of Guermond, Popov, and Tomas [32]–[34], developed for discretization-independent IDP mechanisms and applied across continuous finite elements and other discretization frameworks.  For comprehensive discussions of FCT/IDP principles and algorithms, we refer to [25], [35] and the references therein.

Applying PCP/IDP ideas to RHD is nontrivial.  Although the equations are conservative, neither the fluxes nor the wave speeds admit simple closed-form expressions purely in terms of conservative variables, since they depend on primitive variables obtained through a nonlinear recovery procedure; see Section 2.  Despite these difficulties, major progress has been achieved for special RHD.  Wu and Tang [5] derived the first explicit equivalent characterization of the admissible state set, later generalized to general equations of state (EOS) [36], and developed high-order PCP finite-difference WENO schemes with flux-corrected limiters.  This foundation has supported extensions to unstructured finite-volume methods [37], discontinuous Galerkin (DG) methods [36], [38], and invariant-region/minimum-entropy-principle approaches [39], [40].  Related PCP frameworks have also been established for general relativistic hydrodynamics [41] and coupled with the oscillation-eliminating approach [42] in [43], and PCP techniques have been extended to relativistic MHD [44], [45] as well as nonrelativistic MHD [46]–[48].  For further background, we refer to the reviews [49], [50], the textbook [51], and recent studies [52]–[54].

More recently, the geometric quasilinearization (GQL) framework was introduced in [55], inspired by PCP analyses for (relativistic) MHD and related systems [44]–[48].  GQL converts nonlinear admissibility constraints into equivalent families of linear inequalities parameterized by auxiliary variables, offering a geometric route to PCP design.  However, for special RHD a key practical bottleneck remains: enforcing the nonlinear constraint \(q({\boldsymbol{U}})>0\) (encoding positive pressure and subluminal velocity) in a way that is both sharp and efficient.  Achieving high resolution requires the flux-corrected limiting parameter \(\theta\in[0,1]\) to stay as close to \(1\) as possible, while guaranteeing admissibility at every interface and Runge–Kutta stage.  Even within the GQL framework [55], the limiting step naturally leads to a worst-case optimization over an auxiliary variable \(v_{*}\); existing strategies often rely on discrete sampling or iterative root-finding, whose cost and implementation complexity grow rapidly in multiple space dimensions and may become overly restrictive when the solution approaches the boundary of the admissible set.

To address these issues, we propose a new, efficient PCP flux-limiting framework for special RHD that avoids iterative searches over auxiliary variables in the limiting step.  Building on the GQL characterization and the classical Zalesak-type FCT mechanism, our main contributions are:

• GQL–FCT scalar-style PCP limiting. We integrate the GQL family of linear inequalities into an FCT update to construct a scalar-style limiter that acts directly on conservative variables.  In contrast to classical characteristic FCT limiters for hyperbolic systems [23], [24], our GQL-based PCP limiter does not require characteristic decomposition in the limiting step.  Furthermore, unlike the FCT approach using GQL with a specific convex decomposition for PAMPA schemes [29], the present limiter avoids such a prior convex decomposition; our limiting approach is thus applicable to other high-order numerical frameworks, such as finite difference and finite volume schemes.

• Eigenvalue-based parameter estimation. We resolve the computational bottleneck associated with the nonlinear constraint.  By introducing a rational (stereographic-type) parameterization of the GQL normal vector, we transform the worst-case minimization over GQL auxiliary variables into a generalized Rayleigh-quotient maximization.  Consequently, the limiter parameters are obtained from only a few small symmetric eigenvalue problems (\(2\times2\) in 1D; \((d{+}1)\times(d{+}1)\) in \(d\) dimensions), yielding an explicit, non-iterative procedure.  For 2D and 3D, we further propose relaxed variants that substantially reduce the number of eigenvalue evaluations while retaining robustness (from 16 to 6 in 2D and from 64 to 9 in 3D).

• High-order PCP WENO schemes and numerical validation. We couple the new limiter with standard fifth-order finite-difference WENO reconstruction and strong-stability-preserving Runge–Kutta time integration.  A series of demanding one- to three-dimensional benchmarks, including ultra-relativistic Riemann problems, multidimensional Riemann problems, shock–bubble interaction, and relativistic jet simulations, demonstrate robustness and high-order accuracy.  Comparisons with the classical Wu–Tang PCP limiter [5] indicate that the new estimator provides better resolution.

From the viewpoint of invariant-domain enforcement, convex limiting and its monolithic variants provide a general mechanism to restrict high-order (anti-diffusive) corrections so that the updated degrees of freedom remain in a prescribed convex admissible set; see, e.g., the convex limiting framework of Guermond et al. [33], [34] and the monolithic convex limiting approach of Kuzmin [28].  In contrast, we exploit the specific geometry of the special-RHD admissible set in conservative variables: besides \(D>0\), the nonlinear constraint \(q({\boldsymbol{U}})=E-\sqrt{D^2+|{\boldsymbol{m}}|^2}>0\) is represented exactly as the intersection of linear half-space inequalities \({\boldsymbol{U}}\cdot{\boldsymbol{n}}({\boldsymbol{v}}_*)>0\) for all \({\boldsymbol{v}}_*\in \mathbb{B}_1({\boldsymbol{0}})\).  This makes it possible to embed constraint enforcement into a scalar-style Zalesak-type FCT update by applying lower-bound limiting to the projected quantities \({\boldsymbol{U}}\cdot{\boldsymbol{n}}({\boldsymbol{v}}_*)\), acting directly on conservative variables and avoiding any problem-specific convex decomposition of states in the limiting step.

The remainder of this paper is organized as follows.  Section 2 introduces the governing equations and the GQL characterization of the admissible state set.  Section 3 outlines the base finite-difference WENO discretization.  Section 4 details the construction of the GQL–FCT PCP flux limiter, the eigenvalue-based parameter estimation, and the proof of the PCP property.  Numerical validations are presented in Section 5, followed by conclusions in Section 6.

2 Relativistic Euler system and admissible state set↩︎

The \(d\)-dimensional special relativistic hydrodynamics system can be formulated as a system of hyperbolic conservation laws: \[\label{eq:RHD3D} \frac{\partial {\boldsymbol{U}}}{\partial t} + \sum_{i=1}^{d} \frac{\partial {\boldsymbol{F}}_i({\boldsymbol{U}})}{\partial x_i} = {\boldsymbol{0}},\tag{1}\] where the conservative vector \({\boldsymbol{U}}\) and flux function \({\boldsymbol{F}}_i\) are defined as \[\begin{align} &{\boldsymbol{U}} = \left( D, {\boldsymbol{m}}^\top, E \right)^\top = \left( \rho \gamma, \rho h \gamma^2 {\boldsymbol{v}}^\top, \rho h \gamma^2 - p \right)^\top, \tag{2} \\ &{\boldsymbol{F}}_i = \left( D v_i, v_i {\boldsymbol{m}}^\top + p {\boldsymbol{e}}_i^\top, m_i \right)^\top = \left( \rho \gamma v_i, \rho h \gamma^2 v_i {\boldsymbol{v}}^\top + p {\boldsymbol{e}}_i^\top, \rho h \gamma^2 v_i \right)^\top. \tag{3} \end{align}\] Here, \(D,\;{\boldsymbol{m}}\), and \(E\) correspond to the mass density, momentum density vector, and energy density, respectively, and are collectively referred to as conservative variables. The rest-mass density \(\rho\), fluid velocity \({\boldsymbol{v}} = (v_1, v_2, \dots, v_d)^\top\), and pressure \(p\) are often termed primitive variables. The Lorentz factor \(\gamma = (1 - |{\boldsymbol{v}}|^2)^{-\frac{1}{2}},\) with \({|\boldsymbol{v}|}\) denoting the magnitude of the velocity vector, and \({\boldsymbol{e}}_i\) stands for the \(i\)-th column of the \(d \times d\) identity matrix.

The RHD system 1 is closed by an EOS, which relates the specific enthalpy \(h\) to pressure \(p\) and rest-mass density \(\rho\). To ensure relativistic causality (local sound speed \(c_s<1\)), the EOS must satisfy the inequality [36]: \[h\;\left(\frac{1}{\rho} - \frac{\partial h}{\partial p}(p,\rho)\right) < \frac{\partial h}{\partial \rho}(p,\rho) < 0.\] In this work, we adopt the ideal EOS: \[\label{ID-EOS} h = 1 + \frac{\Gamma p}{(\Gamma -1) \rho} ,\tag{4}\] where \(\Gamma \in (1,2]\) is the adiabatic index. The admissible state set of the RHD system 1 with the ideal EOS 4 is defined as: \[\label{eq:3constraints} \mathcal{G} = \left\{ {\boldsymbol{U}} = (D, {\boldsymbol{m}}^\top, E)^\top \in \mathbb{R}^{d+2} \mid \rho({\boldsymbol{U}}) > 0,\;p({\boldsymbol{U}}) > 0,\;{\boldsymbol{v}}({\boldsymbol{U}}) \in \mathbb{B}_1({\boldsymbol{0}})\subseteq\mathbb{R}^d\right\},\tag{5}\] where \(\mathbb{B}_1({\boldsymbol{0}})\) denotes the \(d\)-dimensional unit ball centered at the origin. An explicit equivalent characterization of \(\mathcal{G}\) was derived in [5]: \[\label{1stequivG} \mathcal{G}^{(1)} = \left\{ {\boldsymbol{U}} = (D, {\boldsymbol{m}}^\top, E)^\top \in \mathbb{R}^{d+2} \mid D > 0,\;q({\boldsymbol{U}}) := E - \sqrt{D^2 + |{\boldsymbol{m}}|^2} > 0 \right\},\tag{6}\] which was further linearized via the GQL framework [55] as: \[\label{eq:1218} \mathcal{G}^{(2)}= \left\{ {\boldsymbol{U}} = (D,{\boldsymbol{m}}^\top,E)^\top \in \mathbb{R}^{d+2} \mid D > 0,\; {\boldsymbol{U}}\cdot {\boldsymbol{n}}_* > 0,\;\forall {\boldsymbol{v}_*} \in \mathbb{B}_1({\boldsymbol{0}}) \right\},\tag{7}\] where \({\boldsymbol{n}}_*:= \big( -\sqrt{1-|{\boldsymbol{v}}_*|^2}, -{\boldsymbol{v}}_*^\top, 1 \big)^\top\), and \({\boldsymbol{v}}_*\) is the auxiliary parameter introduced in the GQL framework to linearize the constraints. As demonstrated in Section 4, this linear representation \(\mathcal{G}^{(2)}\) offers substantial advantages over the original nonlinear form of \(\mathcal{G}\) in PCP analysis for RHD.

Different from the nonrelativistic case, there are no explicit expressions for either the fluxes \({\boldsymbol{F}}_i\) or the primitive quantities in terms of the conservative variables \({\boldsymbol{U}}\) for RHD. In the computation, we have to first recover the primitive variables from the conservative ones before evaluating the fluxes. For a given conservative vector \(\boldsymbol{U}=(D,\boldsymbol{m},E)^\top\in\mathcal{G}\), we recover primitive variables \((\rho,\boldsymbol{v},p)^\top\) for the ideal EOS 4 by solving a nonlinear algebraic equation [5] \[\label{RHDCon2Pri} \Phi_{\boldsymbol{U}}(p) := \frac{p}{\Gamma-1}-E+\frac{\|{\boldsymbol{m}}\|^2}{E+p}+D\sqrt{1-\frac{\|{\boldsymbol{m}}\|^2}{(E+p)^2}} = 0 \quad \text{for}\;\;p>0.\tag{8}\] Once \(p\) is determined from 8 , the velocity vector \({\boldsymbol{v}}\) and rest-mass density \(\rho\) are computed via: \[\boldsymbol{v}=\frac{\boldsymbol{m}}{E+p},\qquad \rho=D\sqrt{1-\|{\boldsymbol{v}}\|^2}.\] To solve 8 robustly, we suggest the provably convergent iterative algorithms for RHD primitive-variable recovery in [37], [56].

3 Review of finite difference WENO methods↩︎

The proposed GQL-based limiting framework is applicable to general high-order finite-difference and finite-volume schemes. Since WENO methods are particularly popular, our numerical experiments use a high-order finite-difference classical WENO-JS scheme [57] with \(\epsilon = 10^{-6}\) as a representative example to demonstrate the effectiveness of the framework; however, the same framework can be applied directly to other high-order discretizations as well, such as WENO-Z [58] and monotonicity-preserving [59] schemes. For completeness, we briefly recall the standard finite-difference WENO discretization for one-dimensional systems of conservation laws. Consider the \(d=1\) case of 1 , \[\label{eq:RHD1D} {\boldsymbol{U}}_t + {\boldsymbol{F}}({\boldsymbol{U}})_x = {\boldsymbol{0}},\tag{9}\] on a uniform grid \(x_i=i\Delta x\). A conservative flux-differencing semi-discretization takes the form \[\label{eq:semiWENO} \frac{d{\boldsymbol{U}}_i}{dt} = -\frac{1}{\Delta x}\Bigl(\widehat{\boldsymbol{F}}_{i+1/2}-\widehat{\boldsymbol{F}}_{i-1/2}\Bigr),\tag{10}\] where \(\widehat{\boldsymbol{F}}_{i+1/2}\) is a high-order accurate numerical flux at the interface \(x_{i+1/2}=(x_i+x_{i+1})/2\).

3.0.0.1 Global Lax–Friedrichs flux splitting

We employ the global Lax–Friedrichs (LF) splitting \[\label{eq:LFsplit} {\boldsymbol{F}}^{\pm}({\boldsymbol{U}})=\frac{1}{2}\Bigl({\boldsymbol{F}}({\boldsymbol{U}}) \pm \alpha{\boldsymbol{U}} \Bigr),\tag{11}\] where \(\alpha\) is the globally defined maximum wave speed, that is, the spectral radius of the Jacobian \({\boldsymbol{A}}({\boldsymbol{U}})=\partial{\boldsymbol{F}}/\partial{\boldsymbol{U}}\) over all grid points, i.e., \(\alpha = \max_j \rho\bigl({\boldsymbol{A}}({\boldsymbol{U}}_j)\bigr)\). Then \({\boldsymbol{F}}({\boldsymbol{U}})={\boldsymbol{F}}^{+}({\boldsymbol{U}})+{\boldsymbol{F}}^{-}({\boldsymbol{U}})\), and the positive/negative parts can be reconstructed using upwind-biased stencils.

3.0.0.2 Characteristic projection and WENO5 reconstruction

Let \({\boldsymbol{L}}_{i+1/2}\) and \({\boldsymbol{R}}_{i+1/2}={\boldsymbol{L}}_{i+1/2}^{-1}\) be the left/right eigenvector matrices of \({\boldsymbol{A}}({\boldsymbol{U}})\) evaluated at a suitable interface state (e.g., a local average of \(\{{\boldsymbol{U}}_j\}\)). Define the characteristic split fluxes \[{\boldsymbol{W}}_{j}^{\pm} := {\boldsymbol{L}}_{i+1/2}\,{\boldsymbol{F}}^{\pm}({\boldsymbol{U}}_j).\] For each component of \({\boldsymbol{W}}_j^{+}\), apply the standard fifth-order WENO reconstruction [57], [60] with a left-biased stencil to obtain the interface value \(\widehat{\boldsymbol{W}}^{+}_{i+1/2}\). Similarly, reconstruct \({\boldsymbol{W}}_j^{-}\) with a right-biased stencil to obtain \(\widehat{\boldsymbol{W}}^{-}_{i+1/2}\). Transform back to physical space, \[\widehat{\boldsymbol{F}}^{\pm}_{i+1/2} := {\boldsymbol{R}}_{i+1/2}\,\widehat{\boldsymbol{W}}^{\pm}_{i+1/2}, \qquad \widehat{\boldsymbol{F}}_{i+1/2} := \widehat{\boldsymbol{F}}^{+}_{i+1/2}+\widehat{\boldsymbol{F}}^{-}_{i+1/2}.\] In multiple spatial dimensions, the above procedure is applied dimension-by-dimension to each directional flux.

3.0.0.3 Time integration

The semi-discrete system 10 is advanced in time by the third-order SSP Runge–Kutta method. The PCP flux limiting procedure developed in Section 4 is applied at each Runge–Kutta stage.

4 A novel PCP flux limiter↩︎

In this section, we will discuss the PCP high-order schemes for one-, two-, and three-dimensional RHD with the ideal EOS 4 . The PCP property is achieved through our novel flux limiter, which is motivated by Zalesak’s limiter in Flux-Corrected Transport (FCT) schemes for scalar partial differential equations [23], [24], [35].

Recently, Wu and Shu proposed the Geometric Quasi Linearization (GQL) framework for PCP analysis, providing new insights into dealing with nonlinear hyperbolic systems, such as RHD 1 . The GQL theory reformulates the nonlinear constraints into equivalent linear ones by introducing additional free parameters. As observed and verified in 7 , we can judge whether the solution lies in \(\mathcal{G}\) by the sign of the scalar \({\boldsymbol{U}} \cdot {\boldsymbol{n}}_*\), which can be integrated with Zalesak’s limiter in scalar equations for RHD, without transforming the system into its local characteristic space or seeking a technical convex decomposition for a sufficient PCP condition.

4.1 One-dimensional case↩︎

We consider the following finite difference FCT scheme for the case \(d = 1\) in 1 .

In what follows, we denote:

• \(\widehat{\mathcal{F}}_{i+1/2}^L\): first-order low-order numerical flux (see 14 for the detailed definition).

• \(\widehat{\mathcal{F}}_{i+1/2}^H\): fifth-order characteristic WENO flux described in Section 3.

• \(\widehat{\mathcal{F}}_{i+1/2}^A := \widehat{\mathcal{F}}_{i+1/2}^H - \widehat{\mathcal{F}}_{i+1/2}^L\): anti-diffusive numerical flux.

For convenience, we omit the subscript of \(x_1\): \[\label{uLF1DRHD} {\boldsymbol{U}}_i^L := {\boldsymbol{U}}_i^n - \frac{\Delta t}{\Delta x}\left(\widehat{\mathcal{F}}_{i+{1}/{2}}^L-\widehat{\mathcal{F}}_{i-{1}/{2}}^L\right),\tag{12}\] \[\label{uPCP1D} {\boldsymbol{U}}_i^{n+1} = {\boldsymbol{U}}_i^L - \frac{\Delta t}{\Delta x}\left(\theta_{i+{1}/{2}}\widehat{\mathcal{F}}_{i+{1}/{2}}^A-\theta_{i-{1}/{2}}\widehat{\mathcal{F}}_{i-{1}/{2}}^A\right).\tag{13}\] In our implementation, the low-order flux is chosen as the first-order global LF (Rusanov) flux \[\label{eq:LF95loworder} \widehat{\mathcal{F}}_{i+{1}/{2}}^L=\frac{1}{2}\Bigl({\boldsymbol{F}}({\boldsymbol{U}}_i^n)+{\boldsymbol{F}}({\boldsymbol{U}}_{i+1}^n)-\alpha\,({\boldsymbol{U}}_{i+1}^n-{\boldsymbol{U}}_i^n)\Bigr)\tag{14}\] with the same LF parameter \(\alpha\) as in 11 . (At each SSP-RK stage, \({\boldsymbol{U}}^n\) is replaced by the corresponding stage solution.) The LF scheme 12 with flux 14 is known to preserve the admissible set for the RHD system under a suitable CFL condition; see, e.g., [5], [46]. Our goal is to find the explicit formula for the coefficients \(\theta_{i+{1}/{2}}\).

At each SSP-RK stage, we enforce the following \(\epsilon\)-admissible set in conservative variables: \[\label{eq:1219} \mathcal{G}_\epsilon^{(1)}=\left\{{\boldsymbol{U}}\in\mathbb{R}^3 \;\big|\;D({\boldsymbol{U}})\geq\epsilon,\;q({\boldsymbol{U}})=E-\sqrt{D^2+m^2}\geq\epsilon\right\}.\tag{15}\] We fix a small parameter \(\epsilon>0\) such that \({\boldsymbol{U}}_i^L\in \mathcal{G}_\epsilon^{(1)}\) for all \(i\). In computation, we choose it stage-by-stage as \(\epsilon=\min\bigl\{10^{-13},\;D_i^L\bigr\}\) or \(\epsilon=\min\bigl\{10^{-13},\;q({\boldsymbol{U}}_i^L)\bigr\}\) for each cell within the computational domain.

The \(q\)-constraint in 15 is nonlinear with respect to \({\boldsymbol{U}}\), which makes it difficult for the PCP analysis. We first find an equivalent linear form following the GQL framework [5], [55].

For the exact GQL characterization with zero lower bound, normalization of the normal vector is immaterial because the half-space inequality is homogeneous under positive rescaling. Proposition 1 states that, for the same \(\epsilon\), \(q({\boldsymbol{U}}) \ge \epsilon\) is exactly equivalent to \({\boldsymbol{U}}\cdot{\boldsymbol{n}}_*(v_*)\geq \epsilon\) for all \({ v}_*\in (-1,1)\) with the specific unnormalized normal \({\boldsymbol{n}}_*(v_*)\) in ?? , we thus do not rescale it; consequently the same lower bound \(\epsilon\) is used consistently in all projected inequalities. If one chooses to normalize \({\boldsymbol{n}}_*(v_*)\) to unit Euclidean length \((\widehat{\boldsymbol{n}}_* = {\boldsymbol{n}}_* / \sqrt{2})\), the equivalent constraint becomes \({\boldsymbol{U}}\cdot\widehat{\boldsymbol{n}}_*\geq\epsilon/\sqrt{2}\). The equivalence is preserved only if the constraint parameter is scaled by the reciprocal of the norm. Our choice of the unnormalized normal vector is purely for notational convenience, as it avoids the introduction of arbitrary scaling factors in the subsequent flux limiter derivation. This normalization property holds generally for the GQL framework in all spatial dimensions. Therefore, it suffices to apply the scalar lower-bound operation in Zalesak’s FCT limiter to the two linear constraints associated with \({\boldsymbol{n}}_D:=(1,0,0)^\top\) and \({\boldsymbol{n}}_*(v_*)\) in ?? . To make the description convenient and clear, we follow Zalesak’s notation [24] and redefine the related quantities for any fixed vector \({\boldsymbol{n}}\). \[\begin{align} P^{-}_{i} &= \min\left(0, -\widehat{\mathcal{F}}^{A}_{i+{1}/{2}}\cdot{\boldsymbol{n}}\right) + \min\left(0, \widehat{\mathcal{F}}^{A}_{i-{1}/{2}}\cdot{\boldsymbol{n}}\right).\\ Q^{-}_{i} &= \frac{\Delta x}{\Delta t}\left(\epsilon - {\boldsymbol{U}}^{L}_{i}\cdot{\boldsymbol{n}}\right),\quad {\boldsymbol{U}}^{L}_{i}\;\text{is defined in }\eqref{uLF1DRHD}.\\ R^{-}_{i} &= \begin{cases} \min\left(1, \frac{Q^{-}_{i}}{P^{-}_{i}}\right), & \text{if}\quad P^{-}_{i} < 0,\\ 1,& \text{otherwise}. \end{cases} \end{align}\]

Step 1: Enforce \(D({\boldsymbol{U}})\geq\epsilon\). This step is identical to the scalar FCT lower-bound limiting. We take \({\boldsymbol{n}} = {\boldsymbol{n}}_D\) (so that \({\boldsymbol{U}}\cdot{\boldsymbol{n}}_D=D({\boldsymbol{U}})\)) and use the parameter \(\epsilon=\min\bigl\{10^{-13},\;D_i^L\bigr\}\). The coefficients of the density flux limiter are computed by \[\label{alphaD1D} \theta_{i+\frac{1}{2}}^D = \begin{cases} R_{i}^- & \text{if}\quad \widehat{\mathcal{F}}_{i+{1}/{2}}^{A}\cdot{\boldsymbol{n}}_D \geq 0, \\ R_{i+1}^- & \text{if}\quad \widehat{\mathcal{F}}_{i+{1}/{2}}^{A}\cdot{\boldsymbol{n}}_D < 0. \end{cases}\tag{16}\] The density flux limiter \(\theta_{i+{1}/{2}}^D\) can be exactly computed in each interface because \(\widehat{\mathcal{F}}_{i+{1}/{2}}^{A,D}\), representing the first component of the anti-diffusive numerical flux in the \(x\) direction, is a scalar.

Step 2: Enforce \({\boldsymbol{U}}\cdot{\boldsymbol{n}_*}\geq \epsilon\iff q({\boldsymbol{U}})=E-\sqrt{D^2+m^2}\geq\epsilon\). Similar to Step 1, we obtain the formulations of the flux limiter for \(q({\boldsymbol{U}})\) \[\label{1DRHDPressure} \theta_{i+\frac{1}{2}}^q = \begin{cases} R_{i}^- & \text{if}\quad \widehat{\mathcal{F}}_{i+{1}/{2}}^{A}\cdot{\boldsymbol{n}}_* \geq 0, \\ R_{i+1}^- & \text{if}\quad \widehat{\mathcal{F}}_{i+{1}/{2}}^{A}\cdot{\boldsymbol{n}}_* < 0. \end{cases}\tag{17}\] In this step, \({\boldsymbol{n}} = {\boldsymbol{n}}_*(v_*)=\left(-\sqrt{1-{ v}_*^2},\;-{v}_*,\;1\right)^\top\) with \({v}_*\in (-1,1)\), and \(\epsilon=\min\bigl\{10^{-13},\;q({\boldsymbol{U}}_i^L)\bigr\}\) is used. The main difference from Step 1 is the presence of the free parameter \({v}_*\). Since \(\widehat{\mathcal{F}}_{i+{1}/{2}}^{A}\cdot{\boldsymbol{n}}_*\) in Step 2 depends on \({v}_*\), its sign cannot be determined once and for all at a fixed interface; hence the sign-based coefficient in 17 cannot be used directly as a uniform limiter for all GQL half-spaces. We therefore solve the following worst-case optimization problem \[\begin{align} (P_1): \mathcal{L}_{i} &= \min\limits_{{v}_*\in (-1,1)}\frac{\frac{\Delta x}{\Delta t}\left({\epsilon} - {\boldsymbol{U}}^{L}_{i}\cdot{\boldsymbol{n}}_*\right)}{\min\left(0, -\widehat{\mathcal{F}}^{A}_{i+{1}/{2}}\cdot{\boldsymbol{n}}_*\right) + \min\left(0, \widehat{\mathcal{F}}^{A}_{i-{1}/{2}}\cdot{\boldsymbol{n}}_*\right)} \nonumber \\ &= \min\limits_{{v}_*\in (-1,1)}\frac{\frac{\Delta x}{\Delta t}\left({\boldsymbol{U}}^{L}_{i}\cdot{\boldsymbol{n}}_*-{\epsilon}\right)}{\max\left(0, \widehat{\mathcal{F}}^{A}_{i+{1}/{2}}\cdot{\boldsymbol{n}}_*\right) + \max\left(0, -\widehat{\mathcal{F}}^{A}_{i-{1}/{2}}\cdot{\boldsymbol{n}}_*\right)} \nonumber \\ &= \frac{\Delta x}{\Delta t}\left(\max\limits_{{ v}_*\in (-1,1)} \frac{\max\left(0, \widehat{\mathcal{F}}^{A}_{i+{1}/{2}}\cdot{\boldsymbol{n}}_*\right) + \max\left(0, -\widehat{\mathcal{F}}^{A}_{i-{1}/{2}}\cdot{\boldsymbol{n}}_*\right)}{{\boldsymbol{U}}^{L}_{i}\cdot{\boldsymbol{n}}_*-{\epsilon}}\right)^{-1}. \label{1DbeforeQuadForm} \end{align}\tag{18}\] As in the scalar FCT coefficient \(R_i^-\), the actual cell-wise coefficient is capped at one. Equivalently, if \(M_i\) denotes the denominator maximum in the last line of 18 , then \(\mathcal{L}_i=1\) when \(M_i=0\), and otherwise \[\mathcal{L}_i=\min\left\{1,\frac{\Delta x}{\Delta t\,M_i}\right\}.\] The uncapped denominator ratio is introduced below only for deriving the quadratic maximization. \[\label{1Doptfun} \mathcal{L}_i({\boldsymbol{n}}_*) := \frac{\max\left(0, \widehat{\mathcal{F}}^{A}_{i+{1}/{2}}\cdot{\boldsymbol{n}}_*\right) + \max\left(0, -\widehat{\mathcal{F}}^{A}_{i-{1}/{2}}\cdot{\boldsymbol{n}}_*\right)}{{\boldsymbol{U}}^{L}_{i}\cdot{\boldsymbol{n}}_*-{\epsilon}}.\tag{19}\]

Since the normal vector \({\boldsymbol{n}}_*(v_*)\) contains the irrational term \(\sqrt{1-|v_{*}|^{2}}\), the minimization problem 18 is algebraically cumbersome to handle directly. To rigorously transform this into a quadratic form, we introduce a stereographic-type rational parameterization of the unit interval. Let \(u_*\) be an auxiliary variable; we define the mapping: \[\label{1Dvstar95para} v_{*} := \frac{2u_{*}}{1+u_{*}^{2}} {~~ withu_* \in (-1,1) \quad \implies \quad {\boldsymbol{n}}_*(u_*) = \left(-\frac{1-u_{*}^{2}}{1+u_{*}^{2}}, -\frac{2u_{*}}{1+u_{*}^{2}}, 1\right)^\top.}\tag{20}\] For the open admissible range, the restriction \(|u_*|<1\) is in one-to-one correspondence for \(v_*\in(-1,1)\), with inverse \(u_*={v_*}/(1+\sqrt{1-v_*^2})\). To simplify the notation and make our presentation clearer, we introduce the following matrix and vector: \[\mathbf{M} = \begin{pmatrix} 1 & 0 & 1 \\ 0 & 2 & 0 \\ -1 & 0 & 1 \end{pmatrix}, \quad \mathbf{S} = \begin{pmatrix} \epsilon \\ 0 \\ \epsilon \end{pmatrix}.\] Then we obtain \[{\boldsymbol{U}}_{i}^L\cdot{\boldsymbol{n}}_*-\epsilon:=\frac{au_{*}^2-bu_{*}+c}{1+u_{*}^2}, \qquad \widehat{\mathcal{F}}_{i\pm{1}/{2}}^A\cdot{\boldsymbol{n}}_*:=\frac{d^{\pm}u_{*}^2-e^{\pm}u_{*}+f^{\pm}}{1+u_{*}^2}, \label{1DUFcdotn}\tag{21}\] where the definitions of the quantities \(a\), \(b\), \(c\), \(d^{\pm}\), \(e^{\pm}\), and \(f^{\pm}\) are listed in Table 1.

Using 21 , we can solve \((P_1)\) by explicitly expressing the optimization function \({\mathcal{L}}_{i}({\boldsymbol{n}}_*)\) as the quadratic form with an extra constraint \(|{ u}_*|<1\). We first introduce some matrix notations, transforming 19 into the matrix-vector formulation. Let \(\star \in \{+,-\}\) and define the sign function \[\tt{sign}(\star) := \begin{cases} 1 & \text{if}\;\star = +,\\ -1 & \text{if}\;\star = -. \end{cases}\] Then we define \[A^{\star}:= {\tt sign}(\star) \begin{pmatrix} d^{\star} & -\frac{e^{\star}}{2} \\ -\frac{e^{\star}}{2} & f^{\star} \end{pmatrix},\qquad B:= \begin{pmatrix} a & -\frac{b}{2} \\ -\frac{b}{2} & c \end{pmatrix},\qquad {\boldsymbol{w}}:= \begin{pmatrix} u_* \\ 1 \end{pmatrix}.\] Since \({\boldsymbol{w}}^\top B{\boldsymbol{w}}={\boldsymbol{U}}^{L}_{i}\cdot{\boldsymbol{n}}_*-\epsilon\geq 0\), the matrix \(B\) is symmetrically positive definite except for the case \({b}^2=4ac\iff q({\boldsymbol{U}}_{i}^L)=\epsilon\). In this borderline case, any anti-diffusive correction may immediately violate \(q({\boldsymbol{U}})\ge\epsilon\), and we simply take \(\mathcal{L}_{i}=0\) (i.e., keep the low-order update). Otherwise, we rewrite \[B = R^\top R,\qquad R = \begin{pmatrix} \sqrt{a} & -\frac{b}{2\sqrt{a}} \\ 0 & \sqrt{c-\frac{{b}^2}{4a}} \end{pmatrix}.\] We can then derive the explicit formulae for matrices \(\widehat{A}^{\star}:=R^{-\top}A^{\star}R^{-1}\): \[\label{Astar} \widehat{A}^{\star} = {\tt sign}(\star) \begin{pmatrix} \frac{d^{\star}}{a} & \frac{bd^{\star}-ae^{\star}}{a\sqrt{4ac-{b}^2}} \\ \frac{bd^{\star}-ae^{\star}}{a\sqrt{4ac-{b}^2}} & \frac{4a^2f^{\star}-2a{b}{e^{\star}}+d^{\star}{b}^2}{a\left(4ac-{b}^2\right)} \end{pmatrix}.\tag{22}\] Based on the above matrix-vector notations, we continue to solve \((P_1)\).

\[\begin{align} {\mathcal{L}}_{i}&=\frac{\Delta x}{\Delta t}\left(\max\limits_{{\boldsymbol{w}}:=({u}_*,1)^\top\in\mathbb{R}^2 \atop {|{ u}_{*}|<1}} \frac{\max\left(0,{\boldsymbol{w}}^\top A^+{\boldsymbol{w}}\right)+\max\left(0,{\boldsymbol{w}}^\top A^-{\boldsymbol{w}}\right)}{{\boldsymbol{w}}^\top B{\boldsymbol{w}}}\right)^{-1}\\ &=\frac{\Delta x}{\Delta t}\left(\max\limits_{{\boldsymbol{z}}:=R{\boldsymbol{w}}\;\in\;\mathbb{R}^2 \atop {|{ u}_{*}|<1}} \frac{\max\left(0,{\boldsymbol{z}}^\top \widehat{A}^+{\boldsymbol{z}}\right)+\max\left(0,{\boldsymbol{z}}^\top \widehat{A}^-{\boldsymbol{z}}\right)}{{\boldsymbol{z}}^\top {\boldsymbol{z}}}\right)^{-1}\\ &=\frac{\Delta x}{\Delta t}\left(\max\limits_{{\boldsymbol{z}}\in\mathbb{R}^2,\;\Vert {\boldsymbol{z}} \Vert_2=1 \atop {|{ u}_{*}|<1}} \max\left(0,{\boldsymbol{z}}^\top \widehat{A}^+{\boldsymbol{z}}\right)+\max\left(0,{\boldsymbol{z}}^\top \widehat{A}^-{\boldsymbol{z}}\right)\right)^{-1}. \end{align}\] The solution to the Rayleigh quotient \(\max\limits_{{\boldsymbol{z}}\in\mathbb{R}^2,\;\Vert {\boldsymbol{z}} \Vert_2=1}{\boldsymbol{z}}^\top\widehat{A}^{\star}{\boldsymbol{z}}\) is the maximal eigenvalue of the matrix \(\widehat{A}^{\star}\). However, this is not a standard Rayleigh quotient, and there are two main differences. On the one hand, we need an extra check for the feasibility of the associated optimizer \({\boldsymbol{w}}=({u}_*,1)^\top\), i.e., \(|{u}_{*}|<1\) or not. On the other hand, for the case \(d = 1\) in 1 , the optimization function 19 is a sum of two quotients composed of the maximum operation, which is clearly different. We adopt the following lemma to simplify the maximum operation.

We can apply Lemma 1 to the Rayleigh quotient to solve \({\mathcal{L}}_{i}\) immediately by taking \(K=2\). \[{\mathcal{L}}_{i} =\frac{\Delta x}{\Delta t}\left( \max\left\{0,\; \max\limits_{{\boldsymbol{z}}\in\mathbb{R}^2,\;\Vert {\boldsymbol{z}} \Vert_2=1 \atop {|{ u}_{*}|<1}}{\boldsymbol{z}}^\top \widehat{A}^+{\boldsymbol{z}},\;\max\limits_{{\boldsymbol{z}}\in\mathbb{R}^2,\;\Vert {\boldsymbol{z}} \Vert_2=1 \atop {|{ u}_{*}|<1}}{\boldsymbol{z}}^\top \widehat{A}^-{\boldsymbol{z}},\; \max\limits_{{\boldsymbol{z}}\in\mathbb{R}^2,\;\Vert {\boldsymbol{z}} \Vert_2=1 \atop {|{ u}_{*}|<1}}{\boldsymbol{z}}^\top \left(\widehat{A}^++\widehat{A}^-\right){\boldsymbol{z}}\right\} \right)^{-1}.\] The remaining step is to solve the Rayleigh quotient with the additional constraint \(|{u}_{*}|<1\): \[\widetilde{\mathcal{L}}_{i}:=\max\limits_{{\boldsymbol{z}}\in\mathbb{R}^2,\;\Vert {\boldsymbol{z}} \Vert_2=1 \atop {|{ u}_{*}|<1}}{\boldsymbol{z}}^\top \widehat{A}\;{\boldsymbol{z}},\] which is the same as the standard Rayleigh quotient except for a further validation of the associated optimizer \({\boldsymbol{w}}=({u}_*,1)^\top\). If the optimizer \({\boldsymbol{w}}\) associated with the maximal eigenvalue satisfies \(|{u}_{*}|<1\), then the solution is consistent with the standard Rayleigh quotient. Otherwise, we compare the remaining feasible eigenvalues with the boundary values. Let \(\widehat{\lambda}_1\) and \(\widehat{\lambda}_2\) be the eigenvalues of the model matrix \(\widehat{A}\), and \(\widehat{\lambda}_1 > \widehat{\lambda}_2\). Then we have \[\widetilde{\mathcal{L}}_{i}= \begin{cases} \max\left\{0,\;\widehat{\lambda}_1\right\} & \text{if}\;\left|{u}_{*,1}\right|<1\\ \max\left\{0,\;\widehat{\lambda}_2,\;\frac{d}{a},\;\frac{f}{c}\right\} & \text{if}\;\left|{u}_{*,1}\right|\geq 1\;\text{and}\;\left|{u}_{*,2}\right|<1\\ \max\left\{0,\;\frac{d}{a},\;\frac{f}{c}\right\} & \text{otherwise} \end{cases},\] where \(u_{*,1}\) and \(u_{*,2}\) are the first component of the optimizers corresponding to \(\widehat{\lambda}_1\) and \(\widehat{\lambda}_2\), respectively. The definitions of \(d\), \(e\), and \(f\) are consistent with \(\widehat{A}\). The eigenvalues of the matrices \(\widehat{A}^{\star}\) (\(\star\in\{+,-\}\)) and \(\widehat{A}^+ + \widehat{A}^-\) are listed below, denoted by \(\widehat{\lambda}_{k}^{\star}\) and \(\widehat{\lambda}_{k}^{*}\) (\(k = 1,2\)): \[\begin{align} \widehat{\lambda}_k^{\star} &= \frac{(-1)^{k-1}\sigma^{\star}+{\tt sign}(\star)(2af^{\star}-{b}{e^{\star}}+2cd^{\star})}{4ac-{b}^2},\tag{23}\\ \widehat{\lambda}_k^{*} &= \frac{(-1)^{k-1}\sigma^{*}+2af^{*}-{b}{e^{*}}+2cd^{*}}{4ac-{b}^2},\tag{24} \end{align}\] where \[\begin{align} d^{*}&:=d^+-d^-,\qquad e^{*}:=e^+ - e^-,\qquad f^{*}:=f^+-f^-,\nonumber\\ \sigma^{\widehat{\star}}&:=2\sqrt{(af^{\widehat{\star}}-cd^{\widehat{\star}})^2-(af^{\widehat{\star}}+cd^{\widehat{\star}}){b}{e^{\widehat{\star}}}+ac{\left(e^{\widehat{\star}}\right)}^2+{b}^2d^{\widehat{\star}} f^{\widehat{\star}}},\quad {\widehat{\star}\in\{\star, *\}}.\label{sigma1D} \end{align}\tag{25}\] The non-negativity of the radicand in 25 is a fundamental prerequisite for the well-posedness of \(\sigma^{\widehat{\star}}\). If \(b=0\), the radicand reduces to \((af^{\widehat{\star}}-cd^{\widehat{\star}})^2+ac\big(e^{\widehat{\star}}\big)^2\ge0\). Thus, in the quadratic-form representation below, we only need to consider \(b\ne0\). We rewrite it as a symmetric quadratic form to reveal its underlying mathematical structure: \[\sigma^{\widehat{\star}} = \sqrt{\mathbf{r}^\top \mathbf{Q} \mathbf{r}}, \quad \text{where} \quad \mathbf{r} = \left( \frac{d^{\widehat{\star}}}{a}, \frac{e^{\widehat{\star}}}{b}, \frac{f^{\widehat{\star}}}{c} \right)^\top,\] and the symmetric matrix \(\mathbf{Q}\) is defined as \[\mathbf{Q} = 2ac \begin{pmatrix} 2ac & -b^2 & b^2-2ac \\ -b^2 & 2b^2 & -b^2 \\ b^2-2ac & -b^2 & 2ac \end{pmatrix}.\] We prove that \(\mathbf{Q}\) is symmetric positive semi-definite for all admissible states by verifying that all the eigenvalues of \({\boldsymbol{Q}}\) are non-negative. The three eigenvalues of \({\boldsymbol{Q}}\) \[\lambda_{\boldsymbol{Q}}^{(1)} = 0,\qquad \lambda_{\boldsymbol{Q}}^{(2)} = 6ab^2c,\qquad \lambda_{\boldsymbol{Q}}^{(3)} = 2ac(4ac-b^2)\] are all non-negative, following the definitions of \(a,b,c\) and the admissible state \({\boldsymbol{U}}_i^L\in\mathcal{G}_{\epsilon}^{(1)}\). Therefore, the quadratic form \(\langle \mathbf{r}, \mathbf{Q} \mathbf{r} \rangle\) is non-negative for all \(\mathbf{r}\), and \(\sigma^{\widehat{\star}}\) is well-defined.

The key components in computing the \(q\)-limiting factor are the eigenvalues of matrices \(\widehat{A}^{\star}\), defined in 22 . It is necessary to explicitly clarify the situation where \(\widehat{A}^{\star}\) is positive definite, positive semi-definite, or indefinite. We summarize it in the following proposition to improve the readability and rigor.

Before ending this subsection, we summarize the one-dimensional PCP property.

4.2 Multi-dimensional case↩︎

The discussions about the case \(d = 1\) in 1 can be easily generalized to \(d = 2\) and \(d = 3\). This is because the GQL tools conveniently transform the dimension of computational variables \({\boldsymbol{U}}\) and \({\widehat{\mathcal{F}}}\) into that of the introduced auxiliary ones. For uniform rectangular cells, we only need to apply the one-dimensional flux limiter in each spatial dimension separately in the finite difference framework. To make the notation clear, we use \(y\) and \(z\) to represent the other two spatial dimensions, rather than \(x_2\) and \(x_3\) in 1 . We analyze the two-dimensional case and provide a brief summary of the three-dimensional case.

The finite difference FCT scheme for the two-dimensional RHD \[{\boldsymbol{U}}_t+{\boldsymbol{F}}({\boldsymbol{U}})_x+{\boldsymbol{G}}({\boldsymbol{U}})_y = 0,\] \[{\boldsymbol{U}} = \begin{pmatrix} D \\ m_x \\ m_y \\ E \end{pmatrix}, \qquad {\boldsymbol{F}}({\boldsymbol{U}}) = \begin{pmatrix} D v_x \\ m_xv_x+p \\ m_yv_x \\ m_x \end{pmatrix}, \qquad {\boldsymbol{G}}({\boldsymbol{U}}) = \begin{pmatrix} D v_y \\ m_xv_y \\ m_yv_y+p \\ m_y \end{pmatrix}\] can be formulated as \[\label{uLF2DRHD} {\boldsymbol{U}}_{ij}^{L} := {\boldsymbol{U}}_{ij}^n -\frac{\Delta t}{\Delta x}\left(\widehat{\mathcal{F}}_{i+\frac{1}{2},j}^L-\widehat{\mathcal{F}}_{i-\frac{1}{2},j}^L\right) -\frac{\Delta t}{\Delta y}\left(\widehat{\mathcal{G}}_{i,j+\frac{1}{2}}^L-\widehat{\mathcal{G}}_{i,j-\frac{1}{2}}^L\right),\tag{28}\] \[\label{2DFE} {\boldsymbol{U}}_{ij}^{n+1} = {\boldsymbol{U}}_{ij}^L -\frac{\Delta t}{\Delta x}\left(\theta_{i+\frac{1}{2},j}\widehat{\mathcal{F}}_{i+\frac{1}{2},j}^A-\theta_{i-\frac{1}{2},j}\widehat{\mathcal{F}}_{i-\frac{1}{2},j}^A\right) -\frac{\Delta t}{\Delta y}\left(\theta_{i,j+\frac{1}{2}}\widehat{\mathcal{G}}_{i,j+\frac{1}{2}}^A-\theta_{i,j-\frac{1}{2}}\widehat{\mathcal{G}}_{i,j-\frac{1}{2}}^A\right).\tag{29}\]

The density limiter is similar to the one-dimensional case, but it requires computation in both the \(x\) and \(y\) directions. Following Zalesak’s notation [24], for a fixed vector \({\boldsymbol{n}}\), we define \[\begin{align} P^{-}_{ij} &= \Delta y\left[\min\left(0, -\widehat{\mathcal{F}}^{A}_{i+\frac{1}{2},j}\cdot{\boldsymbol{n}}\right) + \min\left(0, \widehat{\mathcal{F}}^{A}_{i-\frac{1}{2},j}\cdot{\boldsymbol{n}}\right)\right] +\Delta x\left[\min\left(0, -\widehat{\mathcal{G}}^{A}_{i,j+\frac{1}{2}}\cdot{\boldsymbol{n}}\right) + \min\left(0, \widehat{\mathcal{G}}^{A}_{i,j-\frac{1}{2}}\cdot{\boldsymbol{n}}\right)\right].\\ Q^{-}_{ij} &= \frac{\Delta x\Delta y}{\Delta t}\left(\epsilon - {\boldsymbol{U}}^{L}_{ij}\cdot{\boldsymbol{n}}\right),\quad {\boldsymbol{U}}^{L}_{ij}\;\text{is defined in }\eqref{uLF2DRHD}.\\ R^{-}_{ij} &= \begin{cases} \min\left(1, \frac{Q^{-}_{ij}}{P^{-}_{ij}}\right), & \text{if}\quad P^{-}_{ij} < 0,\\ 1,& \text{otherwise}. \end{cases} \end{align}\] Taking \({\boldsymbol{n}} = {\boldsymbol{n}}_D=(1,0,0,0)^\top\) (so that \({\boldsymbol{U}}\cdot{\boldsymbol{n}}_D=D({\boldsymbol{U}})\)) and using \(\epsilon = \min\bigl\{10^{-13},D({\boldsymbol{U}}_{ij}^L)\bigr\}\), the coefficients of the density flux limiter are \[\begin{align} \theta_{i+\frac{1}{2},j}^D = \begin{cases} R_{ij}^- & \text{if}\;\widehat{\mathcal{F}}_{i+\frac{1}{2},j}^{A}\cdot{\boldsymbol{n}}_D \geq 0, \\ R_{i+1,j}^- & \text{if}\;\widehat{\mathcal{F}}_{i+\frac{1}{2},j}^{A}\cdot{\boldsymbol{n}}_D < 0. \end{cases}, \qquad\qquad \theta_{i,j+\frac{1}{2}}^D = \begin{cases} R_{ij}^- & \text{if}\;\widehat{\mathcal{G}}_{i,j+\frac{1}{2}}^{A}\cdot{\boldsymbol{n}}_D \geq 0, \\ R_{i,j+1}^- & \text{if}\;\widehat{\mathcal{G}}_{i,j+\frac{1}{2}}^{A}\cdot{\boldsymbol{n}}_D < 0. \end{cases}, \end{align}\] where \({\boldsymbol{n}}_D:=(1,0,0,0)^\top\) here.

The limiter for \(q({\boldsymbol{U}})=E-\sqrt{D^2+|{\boldsymbol{m}}|^2}\geq\epsilon\) is a little different. The main difference lies in the dimension of the GQL auxiliary variable \({\boldsymbol{v}}_*\in\mathbb{B}_1({\boldsymbol{0}})\subseteq\mathbb{R}^2\). The GQL representation of the admissible state set in the two-dimensional case is \[\mathcal{G}_\epsilon^{\tt 2D}=\left\{{\boldsymbol{U}}\in\mathbb{R}^4 \mid D({\boldsymbol{U}})\geq\epsilon,\;{\boldsymbol{U}}\cdot{\boldsymbol{n}_*}\geq \epsilon,\;{\boldsymbol{n}_*}=\left(-\sqrt{1-|{\boldsymbol{v}}_*|^2},\;-{\boldsymbol{v}}_*^\top,\;1\right)^\top,\;\forall\;{\boldsymbol{v}}_*\in\mathbb{B}_1({\boldsymbol{0}})\subseteq\mathbb{R}^2\right\}.\] As in the one-dimensional case, we fix a small parameter \(\epsilon>0\) such that the low-order update satisfies \({\boldsymbol{U}}_{ij}^L\in \mathcal{G}_\epsilon^{\tt 2D}\) for all cells. In computation, at each SSP-RK stage we set \(\epsilon=\min\bigl\{10^{-13},\;q({\boldsymbol{U}}_{ij}^L)\bigr\}\). We generalize the stereographic-type transformation in 20 to higher dimensions (a rational parameterization of the unit ball) and obtain \[\label{2Dvstar95para} {\boldsymbol{v}}_*:=\frac{2{\boldsymbol{u}}_{*}}{1+\|{\boldsymbol{u}}_*\|_2^2} {~~ with\|{\boldsymbol{u}}_*\|_2<1 \quad \implies \quad {\boldsymbol{n}}_*({\boldsymbol{u}}_*)=\left(-\frac{1-\|{\boldsymbol{u}}_*\|_2^2}{1+\|{\boldsymbol{u}}_*\|_2^2},\;-\frac{2{\boldsymbol{u}}_{*}^\top}{1+\|{\boldsymbol{u}}_*\|_2^2},\;1\right)^\top,}\tag{30}\] which has the same form as the one-dimensional case because we apply the same parameterization for all velocity components. Therefore, it is also applicable to the three-dimensional RHD. Next comes the key step. We need to solve the following optimization problem \[(P_2): \mathcal{L}_{ij} = \min\limits_{{\boldsymbol{v}}_*\in\mathbb{B}_1({\boldsymbol{0}})\subseteq\mathbb{R}^2} \frac{Q^{-}_{ij}}{P^{-}_{ij}}.\]

Let \({\boldsymbol{w}}:=({\boldsymbol{u}}_*^\top,1)^\top\in\mathbb{R}^3\). Define \(\widehat{A}_x^{\star}:=\Delta y\widetilde{A}_x^{\star}\), \(\widehat{A}_y^{\star}:=\Delta x\widetilde{A}_y^{\star}\), and the matrix \[\label{2DQuadCoeff} \widetilde{A}_s^{\star}:=R^{-\top}A_s^{\star} R^{-1} = {\tt sign}(\star) \begin{pmatrix} \frac{d_s^{\star}}{a}\mathcal{I}_{2\times2} & \frac{d_s^{\star}{\boldsymbol{b}}-a{\boldsymbol{e}}_s^{\star}}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} \\ \frac{d_s^{\star}{\boldsymbol{b}}^\top-a({\boldsymbol{e}}_s^{\star})^\top}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} & \frac{4a^2f_s^{\star}-2a{\boldsymbol{b}}\cdot{\boldsymbol{e}}_s^{\star}+d_s^{\star}|{\boldsymbol{b}}|^2}{a\left(4ac-|{\boldsymbol{b}}|^2\right)} \end{pmatrix},\qquad s = x,y,\quad \star\in\{+,-\}\tag{31}\] with the coefficients listed in Table 2. To simplify the notation and make our presentation clearer, we introduce the following matrix and vector: \[\mathbf{M}_d = \begin{pmatrix} 1 & {\boldsymbol{0}}_d^\top & 1 \\ {\boldsymbol{0}}_d & 2\mathcal{I}_{d\times d} & {\boldsymbol{0}}_d \\ -1 & {\boldsymbol{0}}_d^\top & 1 \end{pmatrix}, \quad \mathbf{S}_d = \begin{pmatrix} \epsilon \\ {\boldsymbol{0}}_d \\ \epsilon \end{pmatrix},\] where \(d\) is consistent with the space dimension. Following the idea in the one-dimensional case and applying Lemma 1 for \(K = 4\), we have

\[\begin{align} \widehat{\mathcal{L}}_{ij}:&=\max\limits_{{\boldsymbol{w}}\in\mathbb{R}^3 \atop {\|{\boldsymbol{u}}_*\|_2<1} } \frac{\Delta y\left[\max\left(0,{\boldsymbol{w}}^\top A_x^+{\boldsymbol{w}}\right)+\max\left(0,{\boldsymbol{w}}^\top A_x^-{\boldsymbol{w}}\right)\right] +\Delta x\left[\max\left(0,{\boldsymbol{w}}^\top A_y^+{\boldsymbol{w}}\right)+\max\left(0,{\boldsymbol{w}}^\top A_y^-{\boldsymbol{w}}\right)\right]}{{\boldsymbol{w}}^\top B{\boldsymbol{w}}}\nonumber \\ &=\max\limits_{{\boldsymbol{z}}:=R{\boldsymbol{w}}\;\in\;\mathbb{R}^3 \atop {\|{\boldsymbol{u}}_*\|_2<1} } \frac{\Delta y\left[\max\left(0,{\boldsymbol{z}}^\top \widetilde{A}_x^+{\boldsymbol{z}}\right)+\max\left(0,{\boldsymbol{z}}^\top \widetilde{A}_x^-{\boldsymbol{z}}\right)\right] +\Delta x\left[\max\left(0,{\boldsymbol{z}}^\top \widetilde{A}_y^+{\boldsymbol{z}}\right)+\max\left(0,{\boldsymbol{z}}^\top \widetilde{A}_y^-{\boldsymbol{z}}\right)\right]}{{\boldsymbol{z}}^\top {\boldsymbol{z}}} \nonumber \\ &=\max\limits_{{\boldsymbol{z}}\in\mathbb{R}^3,\;\Vert {\boldsymbol{z}} \Vert_2=1,\;{\|{\boldsymbol{u}}_*\|_2<1} } \Big\{\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_x^+{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_x^-{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_y^+{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_y^-{\boldsymbol{z}}\right)\Big\},\nonumber \\ &=\max\Big\{0,\; \max\limits_{{\boldsymbol{z}}\in\mathbb{R}^3,\;\Vert {\boldsymbol{z}} \Vert_2=1,\;{\|{\boldsymbol{u}}_*\|_2<1} }{\boldsymbol{z}}^\top \widehat{A}_x^+{\boldsymbol{z}},\;\ldots,\; \max\limits_{{\boldsymbol{z}}\in\mathbb{R}^3,\;\Vert {\boldsymbol{z}} \Vert_2=1,\;{\|{\boldsymbol{u}}_*\|_2<1} }{\boldsymbol{z}}^\top \big(\widehat{A}_x^++\widehat{A}_x^-+\widehat{A}_y^++\widehat{A}_y^-\big) {\boldsymbol{z}}\Big\}.\label{16mat2D} \end{align}\tag{32}\]

From 32 , we see that the solution to \(\widehat{\mathcal{L}}_{ij}\) is the maximum of \(2^4=16\) maximal eigenvalues of associated matrices. Let \({\boldsymbol{s}}:=(s_1,s_2,s_3,s_4)\in\{0,1\}^4\) be an index vector, and \(\widehat{A}_{\boldsymbol{s}}:=s_1\widehat{A}_x^++s_2\widehat{A}_x^-+s_3\widehat{A}_y^++s_4\widehat{A}_y^-\) contains all 16 possible combinations that make a difference in solving 32 . The basic model of \(\widehat{A}_{\boldsymbol{s}}\) can be expressed as \[\label{basicmodel} \widehat{A}_{\boldsymbol{s}}= \begin{pmatrix} \frac{d}{a}\mathcal{I}_{2\times2} & \frac{d{\boldsymbol{b}}-a{\boldsymbol{e}}}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} \\ \frac{d{\boldsymbol{b}}^\top-a{\boldsymbol{e}}^\top}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} & \frac{4a^2f-2a{\boldsymbol{b}}\cdot{\boldsymbol{e}}+|{\boldsymbol{b}}|^2d}{a(4ac-|{\boldsymbol{b}}|^2)} \end{pmatrix} = \begin{pmatrix} \frac{d}{a} & 0 & \frac{b_1d-ae_1}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} \\ 0 & \frac{d}{a} & \frac{b_2d-ae_2}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} \\ \frac{b_1d-ae_1}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} & \frac{b_2d-ae_2}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} & \frac{4a^2f-2a{\boldsymbol{b}}\cdot{\boldsymbol{e}}+|{\boldsymbol{b}}|^2d}{a(4ac-|{\boldsymbol{b}}|^2)} \end{pmatrix},\tag{33}\] where \({\boldsymbol{b}}=(b_1,b_2)^\top,\;{\boldsymbol{e}}=(e_1,e_2)^\top.\) Three eigenvalues of \(\widehat{A}_{\boldsymbol{s}}\) are \[\label{basicmodel:evalue} \widehat{\lambda}_{\boldsymbol{s}}^{(k)} = \frac{(-1)^{k-1}\sigma_{\boldsymbol{s}}+2af-{\boldsymbol{b}}\cdot{\boldsymbol{e}}+2cd}{4ac-|{\boldsymbol{b}}|^2}\;\;(k=1,2),\quad \widehat{\lambda}_{\boldsymbol{s}}^{(3)} = \frac{d}{a},\tag{34}\] where \[\label{sigma2D} \sigma_{\boldsymbol{s}}:=\sqrt{4\Big((af-cd)^2-(af+cd)({\boldsymbol{b}}\cdot{\boldsymbol{e}})+ac|{\boldsymbol{e}}|^2+|{\boldsymbol{b}}|^2df\Big)-(b_1e_2-b_2e_1)^2}.\tag{35}\] The non-negativity of the radicand in 35 follows directly from the symmetric block structure of 33 . Let \[\Delta:=4ac-|{\boldsymbol{b}}|^2>0,\qquad {\boldsymbol{h}}:=\frac{d{\boldsymbol{b}}-a{\boldsymbol{e}}}{a\sqrt{\Delta}},\qquad g:=\frac{4a^2f-2a{\boldsymbol{b}}\cdot{\boldsymbol{e}}+|{\boldsymbol{b}}|^2d}{a\Delta}.\] Then 33 has the eigenvalue \(d/a\) in the subspace orthogonal to \({\boldsymbol{h}}\) (with the obvious multiplicity-two degeneracy if \({\boldsymbol{h}}={\boldsymbol{0}}\)), and its two remaining eigenvalues are those of the real symmetric \(2\times2\) matrix \[\begin{pmatrix} d/a & \|{\boldsymbol{h}}\|_2 \\ \|{\boldsymbol{h}}\|_2 & g \end{pmatrix}.\] Consequently, \[\sigma_{\boldsymbol{s}}^2=\frac{\Delta^2}{4}\left(\left(g-\frac{d}{a}\right)^2+4\|{\boldsymbol{h}}\|_2^2\right)\ge0,\] which is algebraically equivalent to the radicand in 35 . Hence \(\sigma_{\boldsymbol{s}}\) is real and the eigenvalue formula 34 is well-defined for all admissible states. Since \(\sigma_{\boldsymbol{s}}\geq0\), the maximal eigenvalue of \(\widehat{A}_{\boldsymbol{s}}\) is \[\widehat{\lambda}_{\boldsymbol{s}}:=\max\left\{\widehat{\lambda}_{\boldsymbol{s}}^{(1)},\;\widehat{\lambda}_{\boldsymbol{s}}^{(3)}\right\}.\] We are also interested in the optimizer \({\boldsymbol{w}} = ({\boldsymbol{u}}_*^\top, 1)^\top\in\mathbb{R}^3\) when dealing with the additional constraint \(\|{\boldsymbol{u}}_*\|_2<1\). It can be found by multiplying the eigenvector \({\boldsymbol{z}}\) associated with the maximal eigenvalue of \(\widehat{A}_{\boldsymbol{s}}\) since we apply the substitution \({\boldsymbol{z}}=R{\boldsymbol{w}}\) to express \(\widehat{\mathcal{L}}_{ij}\) as a Rayleigh quotient. \[\begin{align} {\boldsymbol{w}}_{\boldsymbol{s}}^{(1)} &= \frac{1}{\sqrt{c-\frac{|{\boldsymbol{b}}|^2}{4a}}}\left(\frac{{\boldsymbol{b}}^\top}{2a}+\frac{(d{\boldsymbol{b}}^\top-a{\boldsymbol{e}}^\top)\widetilde{\sigma}_{\boldsymbol{s}}^{(1)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2},\;1\right)^\top,\tag{36}\\ {\boldsymbol{w}}_{\boldsymbol{s}}^{(2)} &= \frac{1}{\sqrt{c-\frac{|{\boldsymbol{b}}|^2}{4a}}}\left(\frac{{\boldsymbol{b}}^\top}{2a}-\frac{(d{\boldsymbol{b}}^\top-a{\boldsymbol{e}}^\top)\widetilde{\sigma}_{\boldsymbol{s}}^{(2)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2},\;1\right)^\top,\tag{37}\\ {\boldsymbol{w}}_{\boldsymbol{s}}^{(3)} &= \frac{1}{\sqrt{a}}\left(-\frac{b_2d-ae_2}{b_1d-ae_1},\;1,\;0\right)^\top,\tag{38} \end{align}\] where \[\label{basicemodel:evec95sig} \widetilde{\sigma}_{\boldsymbol{s}}^{(k)} = a\sigma_{\boldsymbol{s}}+(-1)^k\left(|{\boldsymbol{b}}|^2d+2a^2f-2acd-a{\boldsymbol{b}}\cdot{\boldsymbol{e}}\right),\quad k = 1,2.\tag{39}\] The maximal eigenvalue may not be a feasible optimal solution if its corresponding eigenvector \({\boldsymbol{w}} = ({\boldsymbol{u}}_*^\top, 1)^\top\) does not satisfy \(\|{\boldsymbol{u}}_*\|_2<1\). We list this quantity for feasibility judgment for 36 , 37 , and 38 : \[\begin{align} \left|{\boldsymbol{u}}_{\boldsymbol{s}}^{(1)}\right|^2&=\left(\frac{b_1}{2a}+\frac{(b_1d-ae_1)\widetilde{\sigma}_{\boldsymbol{s}}^{(1)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}\right)^2+\left(\frac{b_2}{2a}+\frac{(b_2d-ae_2)\widetilde{\sigma}_{\boldsymbol{s}}^{(1)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}\right)^2,\tag{40}\\ \left|{\boldsymbol{u}}_{\boldsymbol{s}}^{(2)}\right|^2&=\left(\frac{b_1}{2a}-\frac{(b_1d-ae_1)\widetilde{\sigma}_{\boldsymbol{s}}^{(2)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}\right)^2+\left(\frac{b_2}{2a}-\frac{(b_2d-ae_2)\widetilde{\sigma}_{\boldsymbol{s}}^{(2)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}\right)^2,\tag{41}\\ \left|{\boldsymbol{u}}_{\boldsymbol{s}}^{(3)}\right|^2&=\left(-\frac{b_2d-ae_2}{b_1d-ae_1}\right)^2+1^2\geq1,\;\text{and the last component of {\boldsymbol{w}}_{\boldsymbol{s}}^{(3)} is not 1}.\tag{42} \end{align}\] All the above formulations correspond to the basic model of \(\widehat{A}_{\boldsymbol{s}}\). For a given index vector \({\boldsymbol{s}}_0:=(s_1,s_2,s_3,s_4)^\top\), we can obtain \(\widehat{A}_{{\boldsymbol{s}}_0}\) by the corresponding choices of \(d\), \(e_1\), \(e_2\), and \(f\). To simplify the notation and make our presentation clearer, we introduce the base coefficient vectors: \[{\boldsymbol{V}}_d:= \begin{pmatrix} \Delta y\;d_x^+\\ -\Delta y\;d_x^-\\ \Delta x\;d_y^+\\ -\Delta x\;d_y^- \end{pmatrix}, \quad {\boldsymbol{V}}_{e_1}:= \begin{pmatrix} \Delta y\;e_x^{+,1}\\ -\Delta y\;e_x^{-,1}\\ \Delta x\;e_y^{+,1}\\ -\Delta x\;e_y^{-,1} \end{pmatrix}, \quad {\boldsymbol{V}}_{e_2}:= \begin{pmatrix} \Delta y\;e_x^{+,2}\\ -\Delta y\;e_x^{-,2}\\ \Delta x\;e_y^{+,2}\\ -\Delta x\;e_y^{-,2} \end{pmatrix}, \quad {\boldsymbol{V}}_f:= \begin{pmatrix} \Delta y\;f_x^+\\ -\Delta y\;f_x^-\\ \Delta x\;f_y^+\\ -\Delta x\;f_y^- \end{pmatrix}.\] Then the coefficients \(d,e_1,e_2,f\) associated with the index vector \({\boldsymbol{s}}_0\) are given by \[c_{\star} = {\boldsymbol{s}}_0\cdot{\boldsymbol{V}}_{c_{\star}} = \sum\limits_{i=1}^4s_i\;({\boldsymbol{V}}_{c_{\star}})_i,\qquad c_{\star}\in\{d,e_1,e_2,f\}.\] This compact representation unifies all 16 possible sign combinations.

Before the end, a special case should be verified. If the eigenvectors 36 and 37 are not well defined, i.e., \(a{\boldsymbol{e}}-d{\boldsymbol{b}}={\boldsymbol{0}}\), we should check whether it is consistent with the general formulation. In this special case, the basic model \(\widehat{A}_{\boldsymbol{s}}\) degenerates into \[\label{basicmodel:special} \widehat{A}_{\boldsymbol{s}}= \begin{pmatrix} \frac{d}{a} & 0 & 0 \\ 0 & \frac{d}{a} & 0 \\ 0 & 0 & \frac{4a^2f-2a{\boldsymbol{b}}\cdot{\boldsymbol{e}}+|{\boldsymbol{b}}|^2d}{a(4ac-|{\boldsymbol{b}}|^2)} \end{pmatrix}.\tag{43}\] Three eigenvalues of \(\widehat{A}_{\boldsymbol{s}}\) are \[\label{basicmodel:specialevalue} \widehat{\lambda}_{\boldsymbol{s}}^{(1)} = \frac{4af-{\boldsymbol{b}}\cdot{\boldsymbol{e}}}{4ac-|{\boldsymbol{b}}|^2},\quad \widehat{\lambda}_{\boldsymbol{s}}^{(2)} = \frac{d}{a},\quad \widehat{\lambda}_{\boldsymbol{s}}^{(3)} = \frac{d}{a},\tag{44}\] and the maximal eigenvalue is given by \[\widehat{\lambda}_{\boldsymbol{s}} = \begin{cases} \widehat{\lambda}_{\boldsymbol{s}}^{(1)} = \frac{4af-{\boldsymbol{b}}\cdot{\boldsymbol{e}}}{4ac-|{\boldsymbol{b}}|^2} & \text{if}\;af>cd\\ \widehat{\lambda}_{\boldsymbol{s}}^{(2)} = \widehat{\lambda}_{\boldsymbol{s}}^{(3)} = \frac{d}{a} & \text{otherwise} \end{cases}.\] The associated (scaled) eigenvectors are \[\begin{align} {\boldsymbol{w}}_{\boldsymbol{s}}^{(1)} &= \frac{1}{\sqrt{c-\frac{|{\boldsymbol{b}}|^2}{4a}}}\left(\frac{{\boldsymbol{b}}^\top}{2a},\;1\right)^\top,\tag{45}\\ {\boldsymbol{w}}_{\boldsymbol{s}}^{(2)} &= \frac{1}{\sqrt{a}}\left(0,\;1,\;0\right)^\top,\tag{46}\\ {\boldsymbol{w}}_{\boldsymbol{s}}^{(3)} &= \frac{1}{\sqrt{a}}\left(1,\;0,\;0\right)^\top.\tag{47} \end{align}\] We also list the quantity for feasibility judgment for 45 , 46 , and 47 : \[\begin{align} \left|{\boldsymbol{u}}_{\boldsymbol{s}}^{(1)}\right|^2 &= \left(\frac{b_1}{2a}\right)^2+\left(\frac{b_2}{2a}\right)^2=\frac{|{\boldsymbol{b}}|^2}{4a^2},\tag{48}\\ \left|{\boldsymbol{u}}_{\boldsymbol{s}}^{(2)}\right|^2 &= 1 \geq 1 \;\text{and the last component of {\boldsymbol{w}}_{\boldsymbol{s}}^{(2)} is not 1},\tag{49}\\ \left|{\boldsymbol{u}}_{\boldsymbol{s}}^{(3)}\right|^2 &= 1 \geq 1 \;\text{and the last component of {\boldsymbol{w}}_{\boldsymbol{s}}^{(3)} is not 1}.\tag{50} \end{align}\] However, computing 16 eigenvalues may be too expensive and complex, and there exist \(2^6=64\) cases if we extend 32 to the three-dimensional case, which are impractical in application. One relaxation approach is to divide \(\widehat{\mathcal{L}}_{ij}\) into two parts, corresponding to two one-dimensional cases that we have already solved. Note that \[\begin{align} \widehat{\mathcal{L}}_{ij}&=\max\limits_{{\boldsymbol{z}}\in\mathbb{R}^3,\;\Vert {\boldsymbol{z}} \Vert_2=1 \atop {\|{\boldsymbol{u}}_*\|_2<1} } \Big\{\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_x^+{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_x^-{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_y^+{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_y^-{\boldsymbol{z}}\right)\Big\}\\ &\leq \max\limits_{{\boldsymbol{z}}\in\mathbb{R}^3,\;\Vert {\boldsymbol{z}} \Vert_2=1 \atop {\|{\boldsymbol{u}}_*\|_2<1} } \Big\{\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_x^+{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_x^-{\boldsymbol{z}}\right)\Big\}+ \max\limits_{{\boldsymbol{z}}\in\mathbb{R}^3,\;\Vert {\boldsymbol{z}} \Vert_2=1 \atop {\|{\boldsymbol{u}}_*\|_2<1} } \Big\{\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_y^+{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_y^-{\boldsymbol{z}}\right)\Big\}. \end{align}\] For convenience, we denote \[\widehat{\mathcal{L}}_{ij}^{x}:=\max\limits_{{\boldsymbol{z}}\in\mathbb{R}^3,\;\Vert {\boldsymbol{z}} \Vert_2=1,\;{\|{\boldsymbol{u}}_*\|_2<1} } \Big\{\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_x^+{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_x^-{\boldsymbol{z}}\right)\Big\} =\max\left\{0,\;\widehat{\lambda}_x^+,\;\widehat{\lambda}_x^-,\;\widehat{\lambda}_x^*\right\},\] \[\widehat{\mathcal{L}}_{ij}^{y}:=\max\limits_{{\boldsymbol{z}}\in\mathbb{R}^3,\;\Vert {\boldsymbol{z}} \Vert_2=1,\;{\|{\boldsymbol{u}}_*\|_2<1} } \Big\{\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_y^+{\boldsymbol{z}}\right) +\max\left(0,{\boldsymbol{z}}^\top \widehat{A}_y^-{\boldsymbol{z}}\right)\Big\} =\max\left\{0,\;\widehat{\lambda}_y^+,\;\widehat{\lambda}_y^-,\;\widehat{\lambda}_y^*\right\},\] where \(\lambda_s^+\), \(\lambda_s^-\), and \(\lambda_s^*\) represent the maximal eigenvalues of \(\widehat{A}_s^+\), \(\widehat{A}_s^-\), and \(\widehat{A}_s^++\widehat{A}_s^-\) (\(s=x\) or \(y\)), respectively. We then obtain the final relaxation formulation \[\label{relax2D} \widehat{\mathcal{L}}_{ij}^{\mathcal{R}}:=\widehat{\mathcal{L}}_{ij}^{x}+\widehat{\mathcal{L}}_{ij}^{y}= \max\left\{0,\;\widehat{\lambda}_x^+,\;\widehat{\lambda}_x^-,\;\widehat{\lambda}_x^*\right\}+\max\left\{0,\;\widehat{\lambda}_y^+,\;\widehat{\lambda}_y^-,\;\widehat{\lambda}_y^*\right\},\tag{51}\] which is significantly more efficient than 32 . Here \(\widehat{\mathcal{L}}\) denotes the maximal denominator ratio in the optimization problem, not the limiter coefficient itself. Thus, for the relaxed 2D estimator, the corresponding cell-wise \(q\)-limiting factor is \[\mathcal{L}_{ij}^{q,\mathcal{R}}:= \begin{cases} 1, & \widehat{\mathcal{L}}_{ij}^{\mathcal{R}}=0,\\[1mm] \min\left\{1,\dfrac{\Delta x\Delta y}{\Delta t}\left(\widehat{\mathcal{L}}_{ij}^{\mathcal{R}}\right)^{-1}\right\}, & \widehat{\mathcal{L}}_{ij}^{\mathcal{R}}>0, \end{cases}\] with the exact estimator obtained by replacing \(\widehat{\mathcal{L}}_{ij}^{\mathcal{R}}\) by the full 16-matrix maximum in 32 . The interface factors are then assigned by the usual adjacent-cell minimum, e.g., \(\theta^q_{i+1/2,j}=\min\{\mathcal{L}_{ij}^{q},\mathcal{L}_{i+1,j}^{q}\}\) and analogously in the \(y\) direction.

Finally, we briefly explain how to implement our flux limiter in the three-dimensional RHD, which is very similar to the two-dimensional case, and only the different parts will be listed. The FCT finite difference scheme for three-dimensional uniform rectangular meshes can be formulated by \[\label{uLF3DRHD} {\boldsymbol{U}}_{ijk}^{L} := {\boldsymbol{U}}_{ijk}^n -\frac{\Delta t}{\Delta x}\left(\widehat{\mathcal{F}}_{i+\frac{1}{2},j,k}^L-\widehat{\mathcal{F}}_{i-\frac{1}{2},j,k}^L\right) -\frac{\Delta t}{\Delta y}\left(\widehat{\mathcal{G}}_{i,j+\frac{1}{2},k}^L-\widehat{\mathcal{G}}_{i,j-\frac{1}{2},k}^L\right) -\frac{\Delta t}{\Delta z}\left(\widehat{\mathcal{H}}_{i,j,k+\frac{1}{2}}^L-\widehat{\mathcal{H}}_{i,j,k-\frac{1}{2}}^L\right),\tag{52}\] \[\begin{align} {\boldsymbol{U}}_{ijk}^{n+1} = {\boldsymbol{U}}_{ijk}^L &-\frac{\Delta t}{\Delta x}\left(\theta_{i+\frac{1}{2},j,k}\widehat{\mathcal{F}}_{i+\frac{1}{2},j,k}^A-\theta_{i-\frac{1}{2},j,k}\widehat{\mathcal{F}}_{i-\frac{1}{2},j,k}^A\right)\\ &-\frac{\Delta t}{\Delta y}\left(\theta_{i,j+\frac{1}{2},k}\widehat{\mathcal{G}}_{i,j+\frac{1}{2},k}^A-\theta_{i,j-\frac{1}{2},k}\widehat{\mathcal{G}}_{i,j-\frac{1}{2},k}^A\right) -\frac{\Delta t}{\Delta z}\left(\theta_{i,j,k+\frac{1}{2}}\widehat{\mathcal{H}}_{i,j,k+\frac{1}{2}}^A-\theta_{i,j,k-\frac{1}{2}}\widehat{\mathcal{H}}_{i,j,k-\frac{1}{2}}^A\right), \end{align}\] where \({\widehat{\mathcal{F}}}\), \({\widehat{\mathcal{G}}}\), and \({\widehat{\mathcal{H}}}\) represent the numerical fluxes in \(x\), \(y\), and \(z\) direction, respectively. The three-dimensional case can share the same framework as the two-dimensional one, except for the definitions of the following quantities related to the spatial dimension \[\begin{align} P^{-}_{ijk} &= \Delta y\Delta z\left[\min\left(0, -\widehat{\mathcal{F}}^{A}_{i+\frac{1}{2},j,k}\cdot{\boldsymbol{n}}_*\right) + \min\left(0, \widehat{\mathcal{F}}^{A}_{i-\frac{1}{2},j,k}\cdot{\boldsymbol{n}}_*\right)\right]\\ &+\Delta x\Delta z\left[\min\left(0, -\widehat{\mathcal{G}}^{A}_{i,j+\frac{1}{2},k}\cdot{\boldsymbol{n}}_*\right) + \min\left(0, \widehat{\mathcal{G}}^{A}_{i,j-\frac{1}{2},k}\cdot{\boldsymbol{n}}_*\right)\right]\\ &+\Delta x\Delta y\left[\min\left(0, -\widehat{\mathcal{H}}^{A}_{i,j,k+\frac{1}{2}}\cdot{\boldsymbol{n}}_*\right) + \min\left(0, \widehat{\mathcal{H}}^{A}_{i,j,k-\frac{1}{2}}\cdot{\boldsymbol{n}}_*\right)\right].\\ Q^{-}_{ijk} &= \frac{\Delta x\Delta y\Delta z}{\Delta t}\left(\epsilon - {\boldsymbol{U}}^{L}_{ijk}\cdot{\boldsymbol{n}}_*\right),\quad {\boldsymbol{U}}^{L}_{ijk}\;\text{is defined in }\eqref{uLF3DRHD}.\\ R^{-}_{ijk} &= \begin{cases} \min\left(1, \frac{Q^{-}_{ijk}}{P^{-}_{ijk}}\right), & \text{if}\quad P^{-}_{ijk} < 0,\\ 1,& \text{otherwise}. \end{cases} \end{align}\]

The formulations of the model matrix \(\widehat{A}_{\boldsymbol{s}}\) and its associated eigenpairs share the same form as 33 , 34 , 36 , 37 , and 38 but only differ in the size of \({\boldsymbol{b}}\) and \({\boldsymbol{e}}\). \[\widehat{A}_{\boldsymbol{s}}= \begin{pmatrix} \frac{d}{a}\mathcal{I}_{3\times 3} & \frac{d{\boldsymbol{b}}-a{\boldsymbol{e}}}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} \\ \frac{d{\boldsymbol{b}}^\top-a{\boldsymbol{e}}^\top}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} & \frac{4a^2f-2a{\boldsymbol{b}}\cdot{\boldsymbol{e}}+|{\boldsymbol{b}}|^2d}{a(4ac-|{\boldsymbol{b}}|^2)} \end{pmatrix} = \begin{pmatrix} \frac{d}{a} & 0 & 0 & \frac{b_1d-ae_1}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} \\ 0 & \frac{d}{a} & 0 & \frac{b_2d-ae_2}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} \\ 0 & 0 & \frac{d}{a} & \frac{b_3d-ae_3}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} \\ \frac{b_1d-ae_1}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} & \frac{b_2d-ae_2}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} & \frac{b_3d-ae_3}{a\sqrt{4ac-|{\boldsymbol{b}}|^2}} & \frac{4a^2f-2a{\boldsymbol{b}}\cdot{\boldsymbol{e}}+|{\boldsymbol{b}}|^2d}{a(4ac-|{\boldsymbol{b}}|^2)} \end{pmatrix},\] where \({\boldsymbol{b}}=(b_1,b_2,b_3)^\top,\;{\boldsymbol{e}}=(e_1,e_2,e_3)^\top.\)

There is one more repeated eigenvalue compared with 34 : \[\widehat{\lambda}_{\boldsymbol{s}}^{(k)} = \frac{(-1)^{k-1}\sigma_{\boldsymbol{s}}+2af-{\boldsymbol{b}}\cdot{\boldsymbol{e}}+2cd}{4ac-|{\boldsymbol{b}}|^2}\;\;(k=1,2),\quad \widehat{\lambda}_{\boldsymbol{s}}^{(3)} = \widehat{\lambda}_{\boldsymbol{s}}^{(4)} = \frac{d}{a},\] where \[\begin{align} &\sigma_{\boldsymbol{s}}:=\sqrt{4\Big((af-cd)^2-(af+cd)({\boldsymbol{b}}\cdot{\boldsymbol{e}})+ac|{\boldsymbol{e}}|^2+|{\boldsymbol{b}}|^2df\Big)-{|{\boldsymbol{b}}\times{\boldsymbol{e}}|^2}},\label{sigma3D}\\ &{{\boldsymbol{b}}\times{\boldsymbol{e}}=\left(b_2e_3-b_3e_2,\;b_3e_1-b_1e_3,\;b_1e_2-b_2e_1\right)^\top.\nonumber} \end{align}\tag{53}\] Similarly, the radicand of \(\sigma_{\boldsymbol{s}}\) in 53 is non-negative by the same block-structure argument. With \(\Delta:=4ac-|{\boldsymbol{b}}|^2>0\), \({\boldsymbol{h}}:=(d{\boldsymbol{b}}-a{\boldsymbol{e}})/(a\sqrt{\Delta})\), and \(g:=(4a^2f-2a{\boldsymbol{b}}\cdot{\boldsymbol{e}}+|{\boldsymbol{b}}|^2d)/(a\Delta)\), the \(4\times4\) matrix above has the repeated eigenvalue \(d/a\) on the two-dimensional subspace orthogonal to \({\boldsymbol{h}}\), while the two non-repeated eigenvalues are determined by the same real symmetric \(2\times2\) block with entries \(d/a\), \(\|{\boldsymbol{h}}\|_2\), and \(g\). Thus \[\sigma_{\boldsymbol{s}}^2=\frac{\Delta^2}{4}\left(\left(g-\frac{d}{a}\right)^2+4\|{\boldsymbol{h}}\|_2^2\right)\ge0,\] which is equivalent to the radicand in 53 . This guarantees the well-posedness of the 3D eigenvalue formulas. The associated eigenvectors are \[\begin{align} {\boldsymbol{w}}_{\boldsymbol{s}}^{(1)} &= \frac{1}{\sqrt{c-\frac{|{\boldsymbol{b}}|^2}{4a}}}\left( \frac{b_1}{2a}+\frac{(b_1d-ae_1)\widetilde{\sigma}_{\boldsymbol{s}}^{(1)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}, \frac{b_2}{2a}+\frac{(b_2d-ae_2)\widetilde{\sigma}_{\boldsymbol{s}}^{(1)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}, \frac{b_3}{2a}+\frac{(b_3d-ae_3)\widetilde{\sigma}_{\boldsymbol{s}}^{(1)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}, 1\right)^\top,\\ {\boldsymbol{w}}_{\boldsymbol{s}}^{(2)} &= \frac{1}{\sqrt{c-\frac{|{\boldsymbol{b}}|^2}{4a}}}\left( \frac{b_1}{2a}-\frac{(b_1d-ae_1)\widetilde{\sigma}_{\boldsymbol{s}}^{(2)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}, \frac{b_2}{2a}-\frac{(b_2d-ae_2)\widetilde{\sigma}_{\boldsymbol{s}}^{(2)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}, \frac{b_3}{2a}-\frac{(b_3d-ae_3)\widetilde{\sigma}_{\boldsymbol{s}}^{(2)}}{2a\|a{\boldsymbol{e}}-d{\boldsymbol{b}}\|_2^2}, 1\right)^\top,\\ {\boldsymbol{w}}_{\boldsymbol{s}}^{(3)} &= \frac{1}{\sqrt{a}}\left(-\frac{b_2d-ae_2}{b_1d-ae_1},1,0,0\right)^\top,\\ {\boldsymbol{w}}_{\boldsymbol{s}}^{(4)} &= \frac{1}{\sqrt{a}}\left(-\frac{b_3d-ae_3}{b_1d-ae_1},0,1,0\right)^\top, \end{align}\] where \(\widetilde{\sigma}_{\boldsymbol{s}}^{(1)}\) and \(\widetilde{\sigma}_{\boldsymbol{s}}^{(2)}\) share the same definition in 39 . The judgment for the feasibility of eigenvectors and the special case of \(\widehat{A}_{\boldsymbol{s}}\) are very similar to the two-dimensional case, and hence we omit them here.

Since solving \(2^6=64\) Rayleigh quotients is impractical, we suggest using the relaxation approach similar to 51 , which is also adopted and shown in the numerical results. \[\label{relax3D} \widehat{\mathcal{L}}_{ijk}^{\mathcal{R}}:= \sum_{s\in\{x,y,z\}}\max\Bigl\{0,\;\widehat{\lambda}_{\boldsymbol{s}}^{+},\;\widehat{\lambda}_{\boldsymbol{s}}^{-},\;\widehat{\lambda}_{\boldsymbol{s}}^{*}\Bigr\},\tag{54}\] where \(\widehat{\lambda}_{\boldsymbol{s}}^{+}\), \(\widehat{\lambda}_{\boldsymbol{s}}^{-}\) and \(\widehat{\lambda}_{\boldsymbol{s}}^{*}\) denote the maximal eigenvalues of \(\widehat{A}_{\boldsymbol{s}}^{+}\), \(\widehat{A}_{\boldsymbol{s}}^{-}\) and \(\widehat{A}_{\boldsymbol{s}}^{+}+\widehat{A}_{\boldsymbol{s}}^{-}\), respectively. The relaxed estimator 54 reduces the number of eigenvalue evaluations from \(2^6=64\) to \(3\times 3=9\) per cell, and the corresponding 3D cell-wise factor is obtained analogously as \(\mathcal{L}_{ijk}^{q,\mathcal{R}}=1\) when \(\widehat{\mathcal{L}}_{ijk}^{\mathcal{R}}=0\), and otherwise \(\mathcal{L}_{ijk}^{q,\mathcal{R}}=\min\{1,(\Delta x\Delta y\Delta z)/(\Delta t\,\widehat{\mathcal{L}}_{ijk}^{\mathcal{R}})\}\).

4.3 Computational cost and implementation remarks↩︎

The new estimator reduces the determination of the \(q\)-limiting factor to a small number of symmetric eigenvalue problems of size \((d{+}1)\times(d{+}1)\). In 1D, the maximization is reduced to three \(2\times2\) eigenvalue problems, and in 2D and 3D, the exact estimators require \(2^4=16\) and \(2^6=64\) eigenvalue evaluations per cell, respectively. The relaxed variants 51 and 54 significantly reduce this number (see Table 3 for details) while remaining robust in all tests reported in Section 5.

In addition, the symmetric matrices arising from the GQL parameterization have a special low-rank structure, and their maximal eigenvalues (and eigenvectors when needed) can be evaluated by closed-form formulas; see 23 –24 in 1D and 34 and 36 in 2D/3D. This avoids iterative root finding and makes the limiter well suited for large-scale simulations.

As noted in Remark 4, the relaxed estimators are conservative lower bounds on the exact limiting factor, which introduces a small amount of additional numerical dissipation in theory. However, extensive numerical benchmarks in Section 5 demonstrate that this dissipation is negligible or acceptable in practice. This minimal accuracy loss is more than offset by the dramatic reduction in computational cost.

This completes the description of our PCP flux limiting framework. To facilitate reproducibility, 7 provides complete pseudocode, and Section 5 lists all computational parameters, including grids, CFL numbers, final times, etc., used in the experiments.

5 Numerical results↩︎

In this section, we present a series of numerical experiments to validate the performance of the proposed flux-limited high-order accurate method for the one-dimensional (1D), two-dimensional (2D), and three-dimensional (3D) RHD. The primary objective is to verify that the implemented flux limiter preserves the designed spatial accuracy when coupled with the classical fifth-order WENO (WENO5) reconstruction. This is examined using smooth problems where the formal convergence order can be measured. A key focus of the tests is to demonstrate the robustness granted by the flux limiter’s ability to enforce the admissible states of the relativistic system. Several challenging benchmark problems, including strong shocks, ultra-relativistic Riemann problems, and multi-dimensional jets, are simulated. Without the limiting procedure, the standard WENO5 would break down due to the generation of nonphysical solutions. Furthermore, comparison results will be shown between our novel flux limiter and the Wu-Tang limiter [5]. These experiments, therefore, underscore the necessity of the PCP mechanism and simultaneously assess the resolution in capturing complex wave structures. For time integration, we employ the third-order strong-stability-preserving Runge–Kutta method. The CFL number is set to 0.4, and the ideal equation of state 4 is used in all cases, with the adiabatic index set to \(\Gamma=5/3\), unless otherwise specified.

For discontinuous solutions such as Riemann problems, formal convergence order does not apply, and the error is concentrated in a narrow region around shocks and contact discontinuities. The primary goal of numerical methods for such problems is to produce sharp, non-oscillatory shock profiles while preserving physical constraints.

6 Conclusions↩︎

We have developed a robust and efficient physical-constraint-preserving (PCP) flux-limiting framework for high-order schemes, exemplified by finite-difference WENO methods, for special relativistic hydrodynamics. The admissible state set in conservative variables is strictly enforced through the equivalent conditions \(D>0\) and \(q({\boldsymbol{U}})=E-\sqrt{D^2+|{\boldsymbol{m}}|^2}>0\), corresponding to positive rest-mass density, positive pressure, and subluminal velocity.

The central idea is to leverage the geometric quasilinearization (GQL) characterization, which represents the nonlinear constraint \(q({\boldsymbol{U}})>0\) exactly as the intersection of linear half-space inequalities \({\boldsymbol{U}}\cdot{\boldsymbol{n}}_*({\boldsymbol{v}}_*)>0\) for all auxiliary variables \({\boldsymbol{v}}_*\in \mathbb{B}_1({\boldsymbol{0}})\). This formulation allows us to embed the RHD constraints into a scalar-style Zalesak-type FCT limiter: by projecting the flux-corrected update onto the relevant normals, the constraint enforcement reduces to lower-bound limiting of scalar quantities while the limiter operates directly on conservative variables.

A critical innovation is the explicit, non-iterative determination of the limiting parameters via a rational stereographic parameterization of the GQL normal vector. This technique transforms the required worst-case minimization over auxiliary variables into a generalized Rayleigh-quotient formulation, allowing the optimal parameters to be obtained by solving small symmetric eigenvalue problems (\(2\times2\) in 1D; \((d+1)\times(d+1)\) in \(d\) dimensions). In multidimensions, relaxed variants further reduce the number of eigenvalue evaluations while preserving a rigorous PCP guarantee. Consequently, the proposed parameter estimator avoids expensive iterative searches over \({\boldsymbol{v}}_*\) and replaces them with closed-form eigenvalue evaluations, resulting in a lightweight limiter suitable for large-scale multidimensional simulations.

Extensive numerical experiments in one to three space dimensions confirm that the resulting schemes robustly enforce physical admissibility, retain fifth-order accuracy for smooth solutions, and sharply resolve strong discontinuities. Future work will focus on less restrictive (yet still explicit) estimators for the \(q\)-limiting factor, extensions to curvilinear or unstructured meshes and more general equations of state, and applications of the GQL–FCT strategy to other hyperbolic systems with nonlinear invariant domains.

Data availability↩︎

The data and code that support the findings of this study are available from the corresponding author upon reasonable request.

Declaration of competing interest↩︎

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement↩︎

This work was partially supported by Science Challenge Project (No. TZ2025007) and the Shenzhen Science and Technology Program (Grant Nos. JCYJ20250604144300001 and RCJC20221008092757098).

7 Complete Pseudocode for the proposed PCP flux Limiter↩︎

In this appendix, we provide complete pseudocode for our physical-constraint-preserving flux-limiting framework. All notation matches the main text.

7.1 Notation Summary↩︎

• \({\boldsymbol{U}}^n\): Conservative state at time \(t^n\)

• \(\widehat{\mathcal{F}}^H\): High-order WENO5 numerical flux

• \(\widehat{\mathcal{F}}^L\): Low-order global Lax-Friedrichs (Rusanov) flux

• \(\widehat{\mathcal{F}}^A = \widehat{\mathcal{F}}^H - \widehat{\mathcal{F}}^L\): Anti-diffusive flux

• \(\Delta x, \Delta y, \Delta z, \Delta t\): Mesh size and time step

• \(\alpha\): Global maximum wave speed

7.2 Main Algorithm: PCP flux Limiter↩︎

Algorithm A.1 GQL-Based PCP Flux-Limiting Framework (1D/2D/3D)

Input: Conservative state \({\boldsymbol{U}}^n\), high-order flux \(\widehat{\mathcal{F}}^H\), low-order flux \(\widehat{\mathcal{F}}^L\), mesh parameters \(\Delta x, \Delta y, \Delta z\), time step \(\Delta t\), CFL number

Output: Updated state \({\boldsymbol{U}}^{n+1}\) satisfying \({\boldsymbol{U}}^{n+1} \in \mathcal{G}_\epsilon\)

1. Step 1: Compute low-order update (guaranteed to be admissible)

   Compute low-order state \({\boldsymbol{U}}^L\) using first-order LF flux: \[{\boldsymbol{U}}^L = {\boldsymbol{U}}^n - \frac{\Delta t}{\Delta x} \left( \widehat{\mathcal{F}}^L_{i+1/2,j,k} - \widehat{\mathcal{F}}^L_{i-1/2,j,k} \right) - \frac{\Delta t}{\Delta y} \left( \widehat{\mathcal{G}}^L_{i,j+1/2,k} - \widehat{\mathcal{G}}^L_{i,j-1/2,k} \right) - \frac{\Delta t}{\Delta z} \left( \widehat{\mathcal{H}}^L_{i,j,k+1/2} - \widehat{\mathcal{H}}^L_{i,j,k-1/2} \right)\] (For 1D, omit \(y,z\) terms; for 2D, omit \(z\) term.)

2. Step 2: Compute anti-diffusive flux

   \(\widehat{\mathcal{F}}^A_{i+1/2} \leftarrow \widehat{\mathcal{F}}^H_{i+1/2} - \widehat{\mathcal{F}}^L_{i+1/2}\)

   (Similarly for \(\widehat{\mathcal{G}}^A, \widehat{\mathcal{H}}^A\) in 2D/3D.)

3. Step 3: Compute density limiter \(\theta^D\) (enforces \(D \geq \epsilon\))

   1. For each cell \(i\) (and \(j,k\) in 2D/3D):

      • \(\epsilon_D \leftarrow \min(10^{-13}, D({\boldsymbol{U}}^L_i))\)

      • Compute \(P^-_i, Q^-_i, R^-_i\) using Zalesak’s notation with \({\boldsymbol{n}}_D = (1, 0, \dots, 0)^\top\)

   2. For each interface \(i+1/2\):

      • If \(\widehat{\mathcal{F}}^A_{i+1/2} \cdot {\boldsymbol{n}}_D \geq 0\): \(\theta^D_{i+1/2} \leftarrow R^-_i\)

      • Else: \(\theta^D_{i+1/2} \leftarrow R^-_{i+1}\)

   (Similarly for \(y,z\) directions in 2D/3D.)

4. Step 4: Compute q-limiter \(\theta^q\) (enforces \(q(U) \geq \epsilon\))

   • If dimension == 1: \(\theta^q \leftarrow \texttt{Compute1DLimiter}({\boldsymbol{U}}^L, \widehat{\mathcal{F}}^A, \Delta x, \Delta t)\)

   • Else: \(\theta^q \leftarrow \texttt{ComputeRelaxedLimiter}({\boldsymbol{U}}^L, \widehat{\mathcal{F}}^A, \widehat{\mathcal{G}}^A, \widehat{\mathcal{H}}^A, \Delta x, \Delta y, \Delta z, \Delta t)\)

5. Step 5: Apply final limiter and update state

   • For each interface: \(\theta_{i+1/2} \leftarrow \min(\theta^D_{i+1/2}, \theta^q_{i+1/2})\)

   • Compute flux-corrected update: \[{\boldsymbol{U}}^{n+1} = {\boldsymbol{U}}^L - \frac{\Delta t}{\Delta x} \left( \theta_{i+1/2} \widehat{\mathcal{F}}^A_{i+1/2} - \theta_{i-1/2} \widehat{\mathcal{F}}^A_{i-1/2} \right)\]

Return: \({\boldsymbol{U}}^{n+1}\)

7.3 Subroutine: Exact 1D q-Limiter↩︎

Algorithm A.2 Compute1DLimiter(\({\boldsymbol{U}}^L, \widehat{\mathcal{F}}^A, \Delta x, \Delta t\))

Input: Low-order state \({\boldsymbol{U}}^L\), anti-diffusive flux \(\widehat{\mathcal{F}}^A\), mesh parameters

Output: q-limiter \(\theta^q\)

1. For each cell \(i\):

   1. Step 1: Extract coefficients

      \({\boldsymbol{U}}^L_i = (D^L, m^L, E^L)^\top\)

      \(\epsilon_q \leftarrow \min(10^{-13}, q({\boldsymbol{U}}^L_i))\)

      \(a \leftarrow E^L + D^L - \epsilon_q\), \(b \leftarrow 2m^L\), \(c \leftarrow E^L - D^L - \epsilon_q\)

      \(d^\pm \leftarrow \widehat{\mathcal{F}}^A_{i\pm1/2} \cdot (1, 0, 1)^\top\), \(e^\pm \leftarrow 2\widehat{\mathcal{F}}^A_{i\pm1/2} \cdot (0, 1, 0)^\top\), \(f^\pm \leftarrow \widehat{\mathcal{F}}^A_{i\pm1/2} \cdot (-1, 0, 1)^\top\)

   2. Step 2: Rational stereographic parameterization and Rayleigh quotient

      Define symmetric matrices: \[B = \begin{pmatrix} a & -b/2 \\ -b/2 & c \end{pmatrix}, \quad {A^\star={\tt sign}(\star)\begin{pmatrix} d^\star & -e^\star/2 \\ -e^\star/2 & f^\star \end{pmatrix},\quad \star\in\{+,-\}}\] Here \({\tt sign}(+)=1\) and \({\tt sign}(-)=-1\).

      Compute Cholesky decomposition \(B = R^\top R\)

      \({\widetilde{A}^\star \leftarrow R^{-\top} A^\star R^{-1},\quad \star\in\{+,-\}}\)

   3. Step 3: Compute maximal eigenvalues

      \(\lambda^+_{\text{max}} \leftarrow \texttt{MaxEigenvalue}(\widetilde{A}^+)\)

      \(\lambda^-_{\text{max}} \leftarrow \texttt{MaxEigenvalue}(\widetilde{A}^-)\)

      \(\lambda^*_{\text{max}} \leftarrow \texttt{MaxEigenvalue}(\widetilde{A}^+ + \widetilde{A}^-)\)

   4. Step 4: Check feasibility of eigenvectors (\(|u_*| < 1\))

      For each eigenvalue \(\lambda\):

      • Compute the corresponding eigenvector \(z\)

      • \({\boldsymbol{w}} \leftarrow R^{-1} {\boldsymbol{z}}\)

      • If the last component of \({\boldsymbol{w}}\) is nonzero, set \(u_* \leftarrow w_1/w_2\); otherwise mark this eigenvector as infeasible

      • If \(|u_*| \geq 1\): Reject this eigenvalue (infeasible)

   5. Step 5: Compute cell-wise limiting factor

      \(M_i \leftarrow \max(0, \lambda^+_{\text{max, feasible}}, \lambda^-_{\text{max, feasible}}, \lambda^*_{\text{max, feasible}},\text{ boundary candidates})\)

      \(\mathcal{L}_i \leftarrow 1\) if \(M_i=0\); otherwise \(\mathcal{L}_i\leftarrow \min\left\{1,\dfrac{\Delta x}{\Delta t\,M_i}\right\}\).

2. Step 6: Compute interface limiting factors

   For each interface \(i+1/2\): \(\theta^q_{i+1/2} \leftarrow \min(\mathcal{L}_i, \mathcal{L}_{i+1})\)

Return: \(\theta^q\)

7.4 Subroutine: Relaxed 2D/3D q-Limiter↩︎

Algorithm A.3 ComputeRelaxedLimiter(\({\boldsymbol{U}}^L, \widehat{\mathcal{F}}^A, \widehat{\mathcal{G}}^A, \widehat{\mathcal{H}}^A, \Delta x, \Delta y, \Delta z, \Delta t\))

Input: Low-order state \({\boldsymbol{U}}^L\), anti-diffusive fluxes in all directions, mesh parameters

Output: Relaxed q-limiter \(\theta^q\)

Note: ComputeDirectionalDenominator returns the weighted maximal denominator estimate in 51 or 54 ; the conversion to the actual limiter coefficient is performed only after summing the directional denominator estimates.

1. For each cell \(i,j,k\):

   Step 1: Extract coefficients for each spatial direction

   Extract \((a, {\boldsymbol{b}}, c)\) from \({\boldsymbol{U}}^L_{ijk}\) (vector \({\boldsymbol{b}}\) for momentum)

   Extract \((d^\pm_x, {\boldsymbol{e}}^\pm_x, f^\pm_x)\) from \(\widehat{\mathcal{F}}^A_{i\pm1/2,j,k}\)

   Extract \((d^\pm_y, {\boldsymbol{e}}^\pm_y, f^\pm_y)\) from \(\widehat{\mathcal{G}}^A_{i,j\pm1/2,k}\)

   If dimension == 3:

   • Extract \((d^\pm_z, {\boldsymbol{e}}^\pm_z, f^\pm_z)\) from \(\widehat{\mathcal{H}}^A_{i,j,k\pm1/2}\)

   • \(M^x_{ijk} \leftarrow \texttt{ComputeDirectionalDenominator}(a, {\boldsymbol{b}}, c, d^\pm_x, {\boldsymbol{e}}^\pm_x, f^\pm_x, \Delta y \Delta z)\)

   • \(M^y_{ijk} \leftarrow \texttt{ComputeDirectionalDenominator}(a, {\boldsymbol{b}}, c, d^\pm_y, {\boldsymbol{e}}^\pm_y, f^\pm_y, \Delta x \Delta z)\)

   • \(M^z_{ijk} \leftarrow \texttt{ComputeDirectionalDenominator}(a, {\boldsymbol{b}}, c, d^\pm_z, {\boldsymbol{e}}^\pm_z, f^\pm_z, \Delta x \Delta y)\)

   • \(M_{ijk} \leftarrow M^x_{ijk} + M^y_{ijk} + M^z_{ijk}\)

   • \(\mathcal{L}_{ijk}\leftarrow 1\) if \(M_{ijk}=0\); otherwise \(\mathcal{L}_{ijk}\leftarrow\min\left\{1,\dfrac{\Delta x\Delta y\Delta z}{\Delta t\,M_{ijk}}\right\}\)

   Else:

   • \(M^x_{ij} \leftarrow \texttt{ComputeDirectionalDenominator}(a, {\boldsymbol{b}}, c, d^\pm_x, {\boldsymbol{e}}^\pm_x, f^\pm_x, \Delta y)\)

   • \(M^y_{ij} \leftarrow \texttt{ComputeDirectionalDenominator}(a, {\boldsymbol{b}}, c, d^\pm_y, {\boldsymbol{e}}^\pm_y, f^\pm_y, \Delta x)\)

   • \(M_{ij} \leftarrow M^x_{ij} + M^y_{ij}\)

   • \(\mathcal{L}_{ij}\leftarrow 1\) if \(M_{ij}=0\); otherwise \(\mathcal{L}_{ij}\leftarrow\min\left\{1,\dfrac{\Delta x\Delta y}{\Delta t\,M_{ij}}\right\}\)

2. Step 2: Compute interface limiting factors

   For each interface: \(\theta^q \leftarrow \min(\mathcal{L}_{\text{left}}, \mathcal{L}_{\text{right}})\)

Return: \(\theta^q\)

7.5 Auxiliary Function: Maximal Eigenvalue↩︎

Algorithm A.4 MaxEigenvalue(\(A\))

Input: Symmetric matrix \(A\) (size 2x2, 3x3, or 4x4)

Output: Maximal eigenvalue \(\lambda_{\text{max}}\)

• Compute eigenvalues using the explicit formulae from the main text

• \(\lambda_{\text{max}} \leftarrow \text{maximum of all eigenvalues}\)

Return: \(\lambda_{\text{max}}\)

References↩︎