Third-Order Relativistic Hydrodynamics Beyond Müller–Israel–Stewart: A Rational Extended Thermodynamics Perspective

1. Introduction

Relativistic dissipative fluid dynamics has become one of the principal theoretical frameworks for describing nonequilibrium transport phenomena in systems ranging from relativistic heavy-ion collisions to astrophysical plasmas and the early Universe. In ultra-relativistic nucleus–nucleus collisions, the formation of a strongly interacting quark–gluon plasma (QGP) has been convincingly established through experiments conducted at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). Hydrodynamic simulations have played a central role in interpreting these experimental observations, demonstrating that the produced matter behaves as an almost perfect fluid possessing an exceptionally small shear viscosity-to-entropy density ratio [1,2,3]. The remarkable success of relativistic hydrodynamics in reproducing collective flow observables has motivated continuous efforts to improve its theoretical foundations and extend its range of validity far from local thermodynamic equilibrium.

The earliest relativistic formulations of dissipative fluid dynamics were proposed independently by Eckart [4] and by Landau and Lifshitz [5]. These first-order theories express dissipative fluxes as algebraic functions of thermodynamic gradients, thereby extending the classical Navier–Stokes equations to relativistic systems. Although successful near equilibrium, first-order relativistic theories are now known to violate causality and exhibit unstable equilibrium states because the governing equations are parabolic rather than hyperbolic [6,7]. These deficiencies severely restrict their applicability to rapidly evolving relativistic systems, particularly those characterized by large gradients and short relaxation times.

To overcome these shortcomings, Israel and Stewart developed a second-order theory based on Grad’s 14-moment approximation to relativistic kinetic theory [8]. By introducing finite relaxation times for dissipative fluxes, the Israel–Stewart formulation restores causality, ensures hyperbolicity under suitable conditions, and has become the standard theoretical framework employed in modern simulations of relativistic heavy-ion collisions. Subsequent developments have substantially refined this approach through systematic derivations from the Boltzmann equation using the Chapman–Enskog expansion, moments methods, and effective field theory techniques [9,10,11,12]. These advances have considerably improved the quantitative description of nonequilibrium transport while maintaining consistency with relativistic kinetic theory.

Despite their widespread success, second-order theories possess important limitations. As the Knudsen number and inverse Reynolds number increase, nonlinear dissipative interactions become increasingly significant, and second-order truncations may fail to accurately capture the evolution of highly nonequilibrium systems. Numerical investigations have demonstrated that higher-order corrections become particularly relevant during the early stages of heavy-ion collisions, where strong gradients and rapid longitudinal expansion dominate the dynamics [10,11,13]. These observations have motivated the development of third-order relativistic hydrodynamic theories capable of extending the domain of validity beyond the second-order approximation.

One of the most influential third-order formulations was proposed by El, Xu, and Greiner [13], who derived evolution equations including cubic contributions to the entropy current directly from relativistic kinetic theory. Later, Jaiswal [10] obtained a complete third-order theory through a Chapman–Enskog expansion of the Boltzmann equation in the relaxation-time approximation. More recently, Denicol and collaborators developed systematic higher-order hydrodynamic equations derived from relativistic kinetic theory that improve agreement with exact solutions of the Boltzmann equation under a wide range of nonequilibrium conditions [9,11]. These kinetic-theory-based formulations have demonstrated substantial improvements in describing pressure anisotropies, entropy production, and dissipative relaxation during rapidly expanding relativistic flows.

An alternative route to relativistic dissipative fluid dynamics is provided by Rational Extended Thermodynamics (RET), which constructs continuum theories directly from balance laws, thermodynamic principles, and mathematical consistency requirements rather than explicit kinetic descriptions. Originating from the pioneering work of Müller and later systematically developed by Müller and Ruggeri, RET enlarges the set of thermodynamic state variables by promoting dissipative fluxes to independent fields satisfying additional balance equations [14,15,16]. The resulting governing equations possess a symmetric hyperbolic structure, satisfy the entropy principle exactly, admit finite signal propagation speeds, and avoid the causality and stability problems associated with classical irreversible thermodynamics. RET has achieved considerable success in describing rarefied gases, polyatomic fluids, chemically reacting mixtures, and high-frequency transport phenomena across a broad range of applications [17,18].

Although RET has been extensively investigated for second-order relativistic fluids, comparatively little attention has been devoted to systematic third-order formulations and their quantitative comparison with contemporary kinetic-theory-based hydrodynamics. In particular, it remains unclear to what extent entropy-based continuum theories reproduce the nonlinear transport mechanisms predicted by higher-order relativistic kinetic theory. Establishing this correspondence is essential for assessing whether RET can serve as a mathematically rigorous continuum alternative to Boltzmann-based closures in strongly nonequilibrium relativistic systems.

The objective of the present work is to address this gap by performing a comprehensive comparison between a third-order Rational Extended Thermodynamics formulation and established third-order relativistic hydrodynamics derived from kinetic theory. The analysis focuses on Bjorken flow, which provides an analytically tractable and physically relevant benchmark for longitudinally expanding relativistic matter [19]. Using identical initial conditions and transport parameters, we compare the evolution of normalized shear stress, pressure anisotropy, effective pressures, entropy production, and dissipative relaxation over a range of shear viscosity-to-entropy density ratios representative of strongly interacting QGP matter. Quantitative error metrics are introduced to evaluate the level of agreement between the two formulations and to identify regimes in which higher-order nonlinear effects become significant.

The principal contribution of this paper is the demonstration that third-order RET reproduces the dominant dissipative dynamics predicted by relativistic kinetic theory while simultaneously preserving the rigorous mathematical structure imposed by the entropy principle and symmetric hyperbolicity. The results provide new evidence that entropy-based continuum theories are capable of accurately describing higher-order relativistic transport phenomena without requiring explicit kinetic closures. Consequently, the proposed framework offers a promising foundation for future developments in multidimensional relativistic fluid dynamics, chemically reacting gas mixtures, relativistic magnetohydrodynamics, and Bayesian identification of higher-order transport coefficients.

The remainder of this paper is organized as follows. Section 2 presents the governing equations of the third-order Rational Extended Thermodynamics formulation together with the associated entropy structure and constitutive relations. Section 2 discusses the correspondence between RET and third-order kinetic theory. Section 3 describes the numerical implementation for relativistic Bjorken flow and the computational methodology. Section 4 presents the numerical comparisons together with quantitative error analyses and physical interpretation of the higher-order transport mechanisms. Finally, Section 5 summarizes the principal findings and outlines several directions for future research.

2. Third-Order Rational Extended Thermodynamics Framework

Relativistic dissipative fluid dynamics seeks to provide a macroscopic description of nonequilibrium transport processes while preserving the fundamental principles of thermodynamics, covariance, and causality. Although relativistic kinetic theory offers a systematic route for deriving transport equations from the Boltzmann equation, continuum formulations based directly on thermodynamic principles possess several important advantages, including mathematical generality, independence from specific microscopic interaction models, and a natural framework for describing systems in which kinetic closures are unavailable or difficult to obtain. Rational Extended Thermodynamics (RET) represents one such continuum theory and has been extensively developed over the past several decades as a mathematically rigorous extension of classical irreversible thermodynamics [14,15,16].

The central idea of RET is to enlarge the thermodynamic state space by promoting dissipative quantities, such as the heat flux and viscous stresses, from constitutive variables to independent thermodynamic fields satisfying their own balance equations. Consequently, the governing equations form a closed system of first-order hyperbolic partial differential equations whose constitutive structure is determined by the entropy principle rather than phenomenological assumptions. Unlike classical Navier–Stokes theory, in which dissipative fluxes respond instantaneously to thermodynamic gradients, RET predicts finite relaxation times and finite propagation speeds for thermal and mechanical disturbances. These properties ensure consistency with special relativity and eliminate the acausal behavior associated with first-order relativistic theories [6,7].

A defining characteristic of Rational Extended Thermodynamics is that the entropy density and entropy flux are regarded as constitutive functions of an enlarged set of state variables. The entropy inequality is then exploited to derive admissible constitutive relations and evolution equations through Liu’s multiplier method or equivalent procedures based on Lagrange multipliers [16,20]. This approach guarantees compatibility with the Second Law of Thermodynamics while simultaneously establishing the symmetric hyperbolic structure required for well-posed initial-value problems.

The present work extends the RET formalism to third order in dissipative quantities in order to investigate nonlinear transport phenomena beyond the conventional Israel–Stewart approximation. Third-order contributions become increasingly important when the Knudsen number or inverse Reynolds number is no longer sufficiently small, as occurs during the early evolution of relativistic heavy-ion collisions and other rapidly expanding nonequilibrium systems [10,11,13]. The resulting formulation retains the exact entropy structure of RET while incorporating higher-order nonlinear couplings that improve agreement with relativistic kinetic theory.

Throughout this work we adopt the metric signature

$$g^{\mu \nu }=\mathrm{diag}(1,-1,-1,-1),$$

(1)

and employ natural units in which

$$c=\hslash =k_{\mathrm{B}}=1.$$

(2)

Greek indices run over spacetime coordinates $(0,1,2,3)$ and repeated indices imply Einstein summation. The hydrodynamic four-velocity satisfies the normalization condition

$$u^{\mu }{u}_{\mu }=1,$$

(3)

while the spatial projection tensor is defined by

$${\Delta }^{\mu \nu }=g^{\mu \nu }-u^{\mu }{u}^{\nu },$$

(4)

which projects tensors onto the three-dimensional space orthogonal to the fluid velocity.

The remainder of this section derives the governing balance equations, establishes the entropy structure of the third-order RET formulation, develops the associated constitutive relations, and demonstrates the correspondence between the resulting transport equations and contemporary third-order relativistic hydrodynamics derived from kinetic theory.

2.1. Balance Laws for Relativistic Dissipative Fluids

The foundation of Rational Extended Thermodynamics is provided by the local conservation laws expressing the balance of particle number, energy, momentum, and entropy. In the absence of external body forces and source terms, particle conservation is governed by

$${\partial }_{\mu }{N}^{\mu }=0,$$

(5)

where

$$N^{\mu }=n{u}^{\mu }+V^{\mu }$$

(6)

is the particle four-current, n denotes the particle number density measured in the local rest frame, and $V^{\mu }$ represents the particle diffusion current satisfying

$$u_{\mu }{V}^{\mu }=0.$$

(7)

For the baryon-free systems considered in relativistic heavy-ion collisions, diffusion effects are negligible, implying

$$V^{\mu }=0,$$

(8)

so that

$$N^{\mu }=n{u}^{\mu }.$$

(9)

The local conservation of energy and momentum is expressed through

$${\partial }_{\mu }{T}^{\mu \nu }=0,$$

(10)

where the energy-momentum tensor is decomposed according to

$$T^{\mu \nu }=\epsilon u^{\mu }{u}^{\nu }-(p+\mathrm{\Pi }){\Delta }^{\mu \nu }+{\pi }^{\mu \nu }+q^{\mu }{u}^{\nu }+q^{\nu }{u}^{\mu }.$$

(11)

Here,

• $\epsilon$ denotes the local energy density,
• p is the equilibrium pressure,
• $\mathrm{\Pi }$ is the bulk viscous pressure,
• $q^{\mu }$ is the heat-flux four-vector,
• ${\pi }^{\mu \nu }$ is the symmetric traceless shear stress tensor.

The dissipative quantities satisfy the orthogonality conditions

$$\begin{array}{cc}\hfill u_{\mu }{q}^{\mu }& =0,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill u_{\mu }{\pi }^{\mu \nu }& =0,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\pi }_{\mu }^{\mu }& =0,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\pi }^{\mu \nu }& ={\pi }^{\nu \mu }.\hfill \end{array}$$

(15)

Projecting Equation (10) parallel and orthogonal to the fluid four-velocity yields the evolution equations for the energy density and momentum,

$$\begin{array}{cc}\hfill D\epsilon +(\epsilon +p+\mathrm{\Pi })\theta -{\pi }^{\mu \nu }{\sigma }_{\mu \nu }+{\nabla }_{\mu }{q}^{\mu }+2q^{\mu }{a}_{\mu }& =0,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill (\epsilon +p+\mathrm{\Pi })a^{\alpha }-{\nabla }^{\alpha }(p+\mathrm{\Pi })+{\Delta }_{\nu }^{\alpha }{\partial }_{\mu }{\pi }^{\mu \nu }+{\Delta }_{\nu }^{\alpha }D{q}^{\nu }+q^{\alpha }\theta +q^{\mu }{\sigma }_{\mu }^{\alpha }& =0,\hfill \end{array}$$

(17)

where

$$\begin{array}{cc}\hfill D& =u^{\mu }{\partial }_{\mu },\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\nabla }^{\mu }& ={\Delta }^{\mu \nu }{\partial }_{\nu },\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \theta & ={\partial }_{\mu }{u}^{\mu },\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill a^{\mu }& =D{u}^{\mu },\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\sigma }^{\mu \nu }& ={\Delta }^{\mu \nu \alpha \beta }{\partial }_{\alpha }{u}_{\beta }\hfill \end{array}$$

(22)

represent the comoving derivative, spatial gradient, expansion scalar, four-acceleration, and shear tensor, respectively.

In classical irreversible thermodynamics the dissipative fluxes are determined algebraically by constitutive laws such as Fourier’s and Navier–Stokes laws. Rational Extended Thermodynamics adopts a fundamentally different viewpoint by considering $\mathrm{\Pi }$, $q^{\mu }$, and ${\pi }^{\mu \nu }$ as independent thermodynamic variables possessing their own evolution equations. Consequently, the complete governing system consists of the conservation laws (5)–(10) together with additional balance equations for the dissipative fields, whose constitutive structure follows directly from the entropy principle.

This enlarged state space constitutes the defining feature of Rational Extended Thermodynamics and provides the mathematical foundation for constructing causal, stable, and symmetric hyperbolic models of relativistic nonequilibrium transport.

2.2. Thermodynamic State Space and the Entropy Principle

The distinguishing feature of Rational Extended Thermodynamics (RET) is the enlargement of the thermodynamic state space to include dissipative fluxes as independent state variables. In contrast to classical irreversible thermodynamics, where constitutive quantities are assumed to respond instantaneously to thermodynamic gradients, RET assigns independent evolution equations to each dissipative field. Consequently, the constitutive structure is determined by the entropy principle together with the requirement that the governing equations form a symmetric hyperbolic system [15,16].

The local thermodynamic state is therefore characterized by the set of variables

$$\mathcal{Z}=\left\{n,\epsilon ,u^{\mu },\mathrm{\Pi },q^{\mu },{\pi }^{\mu \nu }\right\},$$

(23)

where n denotes the particle density, $\epsilon$ is the energy density, $u^{\mu }$ is the fluid four-velocity, $\mathrm{\Pi }$ represents the bulk viscous pressure, $q^{\mu }$ is the heat-flux four-vector, and ${\pi }^{\mu \nu }$ denotes the symmetric traceless shear stress tensor. These variables satisfy the orthogonality and tracelessness conditions introduced in Section 2.1. Within RET, the entropy four-current is regarded as a constitutive quantity depending on the complete thermodynamic state,

$$S^{\mu }=S^{\mu }\left(\mathcal{Z}\right),$$

(24)

rather than solely on equilibrium variables. The Second Law of Thermodynamics requires that the local entropy production be non-negative,

$${\partial }_{\mu }{S}^{\mu }\ge 0,$$

(25)

for every physically admissible thermodynamic process. For systems sufficiently close to local equilibrium, the entropy density may be expanded in powers of the dissipative variables. Retaining contributions up to third order yields

$$\begin{array}{cc}\hfill s=s_{\mathrm{eq}}& -{\displaystyle \frac{{\beta }_0}{2T}}{\mathrm{\Pi }}^2-{\displaystyle \frac{{\beta }_1}{2T}}{q}^{\mu }{q}_{\mu }-{\displaystyle \frac{{\beta }_2}{2T}}{\pi }^{\mu \nu }{\pi }_{\mu \nu }\hfill \\ \hfill & -{\displaystyle \frac{{\gamma }_1}{3T}}{\mathrm{\Pi }}^3-{\displaystyle \frac{{\gamma }_2}{T}}\mathrm{\Pi }{q}^{\mu }{q}_{\mu }-{\displaystyle \frac{{\gamma }_3}{T}}\mathrm{\Pi }{\pi }^{\mu \nu }{\pi }_{\mu \nu }-{\displaystyle \frac{{\gamma }_4}{T}}{\pi }_{\lambda }^{\mu }{\pi }_{\nu }^{\lambda }{\pi }_{\mu }^{\nu }+\mathcal{O}\left(4\right),\hfill \end{array}$$

(26)

where T denotes the absolute temperature, while the coefficients ${\beta }_i$ characterize second-order thermodynamic corrections and ${\gamma }_i$ represent third-order nonlinear couplings. These coefficients are material-dependent functions of the equilibrium state and constitute higher-order transport parameters within the RET framework. Equation (26) extends the quadratic entropy density employed in Israel–Stewart theory by incorporating cubic invariants constructed from the dissipative fields. These higher-order contributions become increasingly important when the Knudsen number or inverse Reynolds number is no longer asymptotically small, thereby providing a more accurate representation of strongly nonequilibrium transport processes. The associated entropy four-current is written as

$$S^{\mu }=s{u}^{\mu }+{\displaystyle \frac{q^{\mu }}{T}}+Q_{\left(2\right)}^{\mu }+Q_{\left(3\right)}^{\mu },$$

(27)

where the first term represents the equilibrium entropy transport, the second corresponds to the classical heat contribution, and $Q_{\left(2\right)}^{\mu }$ and $Q_{\left(3\right)}^{\mu }$ denote second- and third-order nonequilibrium corrections. A convenient representation is

$$\begin{array}{cc}\hfill Q_{\left(2\right)}^{\mu }& =-{\displaystyle \frac{{\alpha }_0}{T}}\mathrm{\Pi }{q}^{\mu }-{\displaystyle \frac{{\alpha }_1}{T}}{\pi }^{\mu \nu }{q}_{\nu },\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill Q_{\left(3\right)}^{\mu }& =-{\displaystyle \frac{{\delta }_1}{T}}{\mathrm{\Pi }}^2{q}^{\mu }-{\displaystyle \frac{{\delta }_2}{T}}\mathrm{\Pi }{\pi }^{\mu \nu }{q}_{\nu }-{\displaystyle \frac{{\delta }_3}{T}}{\pi }^{\mu \lambda }{\pi }_{\lambda \nu }{q}^{\nu },\hfill \end{array}$$

(29)

where ${\alpha }_i$ and ${\delta }_i$ are higher-order coupling coefficients describing nonlinear entropy transport. The entropy inequality is evaluated by substituting the balance equations derived in Section 2.1 into the divergence of Equation (27). After considerable algebra one obtains

$$\begin{array}{cc}\hfill {\partial }_{\mu }{S}^{\mu }=& -{\displaystyle \frac{\mathrm{\Pi }}{T}}\left[{\mathcal{X}}_{\mathrm{\Pi }}\right]-{\displaystyle \frac{q^{\mu }}{T}}\left[{\mathcal{X}}_{\mu }^{\left(q\right)}\right]-{\displaystyle \frac{{\pi }^{\mu \nu }}{T}}\left[{\mathcal{X}}_{\mu \nu }^{\left(\pi \right)}\right]+\mathcal{O}\left(4\right),\hfill \end{array}$$

(30)

where the quantities ${\mathcal{X}}_{\mathrm{\Pi }}$, ${\mathcal{X}}_{\mu }^{\left(q\right)}$, and ${\mathcal{X}}_{\mu \nu }^{\left(\pi \right)}$ represent generalized thermodynamic forces associated with the bulk pressure, heat flux, and shear stress, respectively.

The Second Law requires that Equation (30) remain non-negative for arbitrary admissible thermodynamic processes. Following the standard RET procedure, this requirement is satisfied by assuming linear relations between thermodynamic fluxes and forces,

$$\begin{array}{cc}\hfill {\mathcal{X}}_{\mathrm{\Pi }}& =-{\displaystyle \frac{\mathrm{\Pi }}{\zeta }},\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\mathcal{X}}_{\mu }^{\left(q\right)}& =-{\displaystyle \frac{q_{\mu }}{\kappa T}},\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\mathcal{X}}_{\mu \nu }^{\left(\pi \right)}& =-{\displaystyle \frac{{\pi }_{\mu \nu }}{2\eta }},\hfill \end{array}$$

(33)

where $\zeta$, $\kappa$, and $\eta$ denote the bulk viscosity, thermal conductivity, and shear viscosity, respectively. The entropy production therefore assumes the manifestly positive form

$${\partial }_{\mu }{S}^{\mu }={\displaystyle \frac{{\mathrm{\Pi }}^2}{\zeta T}}+{\displaystyle \frac{q^{\mu }{q}_{\mu }}{\kappa T^2}}+{\displaystyle \frac{{\pi }^{\mu \nu }{\pi }_{\mu \nu }}{2\eta T}}+\mathcal{O}\left(3\right)\ge 0,$$

(34)

provided that

$$\eta >0,\zeta >0,\kappa >0.$$

(35)

Equation (34) constitutes the central thermodynamic constraint of the RET formulation. Unlike phenomenological constitutive theories, the admissible evolution equations are not postulated but are derived directly from the entropy inequality. Consequently, thermodynamic consistency, causality, and nonlinear stability are embedded within the mathematical structure of the governing equations.

In the following subsection, these thermodynamic restrictions are employed to derive the third-order constitutive equations governing the evolution of the dissipative fields. The resulting balance laws extend the Israel–Stewart formulation by incorporating nonlinear entropy-generated couplings while preserving the symmetric hyperbolic structure characteristic of Rational Extended Thermodynamics.

2.3. Third-Order Constitutive Relations

Having established the enlarged thermodynamic state space and the entropy inequality, we now derive the constitutive equations governing the evolution of the dissipative variables. In Rational Extended Thermodynamics (RET), these equations are not introduced phenomenologically but emerge from the simultaneous requirements of local entropy production, covariance, frame indifference, and symmetric hyperbolicity [15,16].

Unlike classical irreversible thermodynamics, where constitutive laws relate dissipative fluxes algebraically to thermodynamic gradients, RET promotes the bulk viscous pressure, heat flux, and shear stress tensor to independent thermodynamic fields satisfying balance-type evolution equations. Consequently, irreversible processes evolve over finite relaxation times, thereby preserving causality and eliminating the instantaneous propagation of disturbances characteristic of parabolic transport theories.

The constitutive equations are obtained by requiring that the entropy production derived in Equation (34) remain positive for arbitrary thermodynamic processes. Following the standard RET procedure, the generalized thermodynamic forces are assumed to depend linearly on the dissipative variables while allowing nonlinear corrections generated by the third-order entropy expansion. The resulting evolution equation for the bulk viscous pressure is written as

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{\Pi }}D\mathrm{\Pi }+\mathrm{\Pi }=& -\zeta \theta -{\delta }_{\mathrm{\Pi }\mathrm{\Pi }}\mathrm{\Pi }\theta +{\lambda }_{\mathrm{\Pi }q}{\nabla }_{\mu }{q}^{\mu }+{\lambda }_{\mathrm{\Pi }\pi }{\pi }^{\mu \nu }{\sigma }_{\mu \nu }\hfill \\ \hfill & -{\chi }_1{\mathrm{\Pi }}^2-{\chi }_2{q}^{\mu }{q}_{\mu }-{\chi }_3{\pi }^{\mu \nu }{\pi }_{\mu \nu }+\mathcal{O}\left(4\right),\hfill \end{array}$$

(36)

where ${\tau }_{\mathrm{\Pi }}$ denotes the bulk relaxation time, $\zeta$ is the bulk viscosity, and the coefficients ${\delta }_{\mathrm{\Pi }\mathrm{\Pi }}$, ${\lambda }_{\mathrm{\Pi }q}$, ${\lambda }_{\mathrm{\Pi }\pi }$, and ${\chi }_i$ describe nonlinear third-order interactions generated by the entropy expansion. Similarly, the evolution equation for the heat-flux vector becomes

$$\begin{array}{cc}\hfill {\tau }_q{\Delta }_{\nu }^{\mu }D{q}^{\nu }+q^{\mu }=& -\kappa \left({\nabla }^{\mu }T+T{a}^{\mu }\right)-{\delta }_{qq}{q}^{\mu }\theta -{\lambda }_{q\mathrm{\Pi }}\mathrm{\Pi }{a}^{\mu }\hfill \\ \hfill & +{\lambda }_{q\pi }{\pi }^{\mu \nu }{a}_{\nu }-{\chi }_4\mathrm{\Pi }{q}^{\mu }-{\chi }_5{\pi }^{\mu \nu }{q}_{\nu }+\mathcal{O}\left(4\right),\hfill \end{array}$$

(37)

where ${\tau }_q$ denotes the thermal relaxation time and $\kappa$ represents the thermal conductivity. The shear stress tensor satisfies

$$\begin{array}{cc}\hfill {\tau }_{\pi }{\Delta }_{\alpha \beta }^{\mu \nu }D{\pi }^{\alpha \beta }+{\pi }^{\mu \nu }=& 2\eta {\sigma }^{\mu \nu }-{\delta }_{\pi \pi }{\pi }^{\mu \nu }\theta +2{\tau }_{\pi }{\pi }_{\lambda }^{\langle \mu }{\omega }^{\nu \rangle \lambda }\hfill \\ \hfill & -{\tau }_{\pi \pi }{\pi }_{\lambda }^{\langle \mu }{\sigma }^{\nu \rangle \lambda }+{\lambda }_{\pi \mathrm{\Pi }}\mathrm{\Pi }{\sigma }^{\mu \nu }+{\lambda }_{\pi q}{q}^{\langle \mu }{a}^{\nu \rangle }\hfill \\ \hfill & -{\chi }_6\mathrm{\Pi }{\pi }^{\mu \nu }-{\chi }_7{\pi }_{\lambda }^{\langle \mu }{\pi }^{\nu \rangle \lambda }-{\chi }_8{q}^{\langle \mu }{q}^{\nu \rangle }+\mathcal{O}\left(4\right),\hfill \end{array}$$

(38)

where

$${\omega }^{\mu \nu }={\displaystyle \frac{1}{2}}{\Delta }^{\mu \alpha }{\Delta }^{\nu \beta }\left({\partial }_{\alpha }{u}_{\beta }-{\partial }_{\beta }{u}_{\alpha }\right)$$

(39)

is the relativistic vorticity tensor and

$${\Delta }_{\alpha \beta }^{\mu \nu }={\displaystyle \frac{1}{2}}\left({\Delta }_{\alpha }^{\mu }{\Delta }_{\beta }^{\nu }+{\Delta }_{\beta }^{\mu }{\Delta }_{\alpha }^{\nu }\right)-{\displaystyle \frac{1}{3}}{\Delta }^{\mu \nu }{\Delta }_{\alpha \beta }$$

(40)

denotes the traceless projection operator. The first terms on the right-hand side of Equations (36)–(38) coincide with the familiar Israel–Stewart relaxation equations. The remaining terms represent genuine third-order contributions arising from the cubic entropy expansion introduced in Section 2.2. These nonlinear couplings describe interactions between different dissipative mechanisms and become increasingly important as the system evolves away from local equilibrium.

In particular, the coefficients ${\chi }_i$ quantify cubic interactions among bulk viscosity, heat conduction, and shear stress. Such terms are absent in first-order relativistic hydrodynamics and appear only partially in kinetic-theory-based third-order formulations. Within RET, however, they arise naturally as direct consequences of the entropy principle and therefore possess a clear thermodynamic interpretation.

An important property of the present formulation is that the evolution equations remain first-order in space and time. Consequently, the complete governing system preserves the symmetric hyperbolic structure characteristic of Rational Extended Thermodynamics. This guarantees finite characteristic propagation speeds, local well-posedness of the initial-value problem, and compatibility with relativistic causality. Near local thermodynamic equilibrium, where

$$\left|\mathrm{\Pi }\right|,\left|q^{\mu }\right|,\left|{\pi }^{\mu \nu }\right|\ll p,$$

(41)

all cubic contributions become asymptotically negligible. Equations (36)–(38) therefore reduce to the standard Israel–Stewart equations,

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{\Pi }}D\mathrm{\Pi }+\mathrm{\Pi }& =-\zeta \theta ,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\tau }_q{\Delta }_{\nu }^{\mu }D{q}^{\nu }+q^{\mu }& =-\kappa ({\nabla }^{\mu }T+T{a}^{\mu }),\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill {\tau }_{\pi }{\Delta }_{\alpha \beta }^{\mu \nu }D{\pi }^{\alpha \beta }+{\pi }^{\mu \nu }& =2\eta {\sigma }^{\mu \nu },\hfill \end{array}$$

(44)

demonstrating that the RET formulation is fully consistent with established second-order relativistic hydrodynamics in the near-equilibrium limit.

The third-order terms therefore constitute nonlinear corrections extending the range of validity of relativistic hydrodynamics toward regimes characterized by large Knudsen numbers and strong dissipative effects. These additional couplings play a central role in accurately describing the early stages of relativistic heavy-ion collisions, rapidly expanding astrophysical plasmas, and other strongly nonequilibrium relativistic flows. In the following subsection, the transport coefficients entering Equations (36)–(38) are discussed together with their physical interpretation and their correspondence with coefficients obtained from relativistic kinetic theory.

2.4. Symmetric Hyperbolic Structure

A defining feature of Rational Extended Thermodynamics is that the governing equations constitute a first-order symmetric hyperbolic system of balance laws. This mathematical structure guarantees finite propagation speeds, local well-posedness of the initial-value problem, uniqueness of sufficiently smooth solutions, and compatibility with relativistic causality [15,16,21,22,23]. Unlike classical Navier–Stokes theory, whose parabolic character leads to instantaneous propagation of disturbances, RET predicts finite characteristic wave speeds determined by the thermodynamic state of the medium.

The complete set of governing variables introduced in Section 2.1, Section 2.2 and Section 2.3 may be collected into the state vector

$$\mathbf{U}={\left(n,\epsilon ,u^{\mu },\mathrm{\Pi },q^{\mu },{\pi }^{\mu \nu }\right)}^T,$$

(45)

where each component satisfies a first-order balance equation. The governing system may therefore be written compactly as

$${\mathbf{A}}^0\left(\mathbf{U}\right){\partial }_t\mathbf{U}+{\mathbf{A}}^i\left(\mathbf{U}\right){\partial }_i\mathbf{U}=\mathbf{S}\left(\mathbf{U}\right),$$

(46)

where

• ${\mathbf{A}}^0$ is the temporal coefficient matrix,
• ${\mathbf{A}}^i$ denote the spatial flux Jacobians,
• $\mathbf{S}$ contains relaxation and nonlinear source terms.

For the system to be symmetric hyperbolic, the coefficient matrices must satisfy

$$\begin{array}{cc}\hfill {\left({\mathbf{A}}^0\right)}^T& ={\mathbf{A}}^0,\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\left({\mathbf{A}}^i\right)}^T& ={\mathbf{A}}^i,\hfill \end{array}$$

(48)

with

$${\mathbf{A}}^0>0,$$

(49)

that is, the temporal matrix is symmetric and positive definite. Within RET these properties follow directly from the convexity of the entropy density. Introducing the entropy potential

$$h=-h^{*}=-\rho s,$$

(50)

where $\rho$ denotes the rest-mass density and s is the specific entropy, the entropy variables are defined by

$$\mathbf{W}={\displaystyle \frac{\partial h}{\partial \mathbf{U}}}.$$

(51)

The Jacobian

$$\mathbf{H}={\displaystyle \frac{{\partial }^{2}h}{\partial {\mathbf{U}}^2}}$$

(52)

is the Hessian of the entropy potential. The strict convexity of the entropy density implies

$$\mathbf{H}>0,$$

(53)

which establishes a one-to-one mapping between conservative variables and entropy variables. Consequently,

$${\mathbf{A}}^0=\mathbf{H},$$

(54)

is automatically symmetric positive definite. The governing equations may therefore be rewritten in entropy variables as

$$\mathbf{H}{\partial }_t\mathbf{W}+{\mathbf{K}}^i{\partial }_i\mathbf{W}=\mathbf{R},$$

(55)

where

$${\left({\mathbf{K}}^i\right)}^T={\mathbf{K}}^i.$$

(56)

Equation (55) is the canonical Godunov form of a symmetric hyperbolic system. An immediate consequence is the existence of real characteristic speeds. Considering propagation in the direction of a unit normal vector $n_i$, one obtains the generalized eigenvalue problem

$$det\left(\lambda {\mathbf{A}}^0-n_i{\mathbf{A}}^i\right)=0.$$

(57)

Since ${\mathbf{A}}^0$ is positive definite and the flux matrices are symmetric, all eigenvalues

$${\lambda }_k\in \mathbb{R}$$

(58)

are real. The corresponding characteristic fields are complete, ensuring that disturbances propagate along finite characteristic surfaces rather than instantaneously throughout the domain. For relativistic fluids, the characteristic speeds satisfy

$$\left|{\lambda }_k\right|\le c,$$

(59)

where c denotes the speed of light (equal to unity in the present units). Consequently, the theory satisfies Einstein causality. Another important consequence of the symmetric hyperbolic structure is the local well-posedness of the initial-value problem. Given sufficiently smooth initial data,

$$\mathbf{U}(x,0)={\mathbf{U}}_0\left(x\right),$$

(60)

there exists a unique local solution depending continuously upon the prescribed initial conditions. This property follows directly from the classical theory of Friedrichs and Lax for first-order symmetric hyperbolic systems [22]. The third-order nonlinear couplings introduced in Section 2.3 modify only the relaxation source vector

$$\mathbf{S}={\mathbf{S}}^{\left(1\right)}+{\mathbf{S}}^{\left(2\right)}+{\mathbf{S}}^{\left(3\right)},$$

(61)

while leaving the principal part of the governing equations unchanged. Consequently, the inclusion of cubic entropy terms preserves the hyperbolic character of the system. This property represents a significant advantage over higher-order gradient expansions whose mathematical structure may become ill posed for sufficiently large gradients.

It is important to emphasize that the preservation of symmetric hyperbolicity is not imposed as an additional assumption but emerges naturally from the entropy principle. The entropy density simultaneously serves as a thermodynamic potential and as a mathematical generating function for the governing equations. This dual role is one of the defining characteristics of Rational Extended Thermodynamics and distinguishes it fundamentally from phenomenological extensions of relativistic Navier–Stokes theory.

The hyperbolic structure established above forms the mathematical foundation for the numerical simulations presented in Section 3. Because the governing equations possess finite characteristic speeds and a convex entropy function, they are particularly well suited for high-resolution finite-volume and discontinuous Galerkin discretizations that preserve both stability and thermodynamic consistency.

2.5. Transport Coefficients and Correspondence with Third-Order Kinetic Theory

The governing equations derived in the preceding subsections contain a collection of transport coefficients that characterize irreversible processes and nonlinear dissipative interactions. Within Rational Extended Thermodynamics (RET), these coefficients are not arbitrary phenomenological parameters but constitutive quantities constrained by the entropy principle, material frame indifference, and the requirement that the governing equations remain symmetric hyperbolic. Their values depend on the underlying equation of state and the microscopic properties of the medium, while their mathematical admissibility is determined by the convexity of the entropy function and the positivity of entropy production. The first class of coefficients consists of the classical transport parameters,

$${\mathcal{T}}_1=\left\{\eta ,\zeta ,\kappa \right\},$$

(62)

where $\eta$ denotes the shear viscosity, $\zeta$ the bulk viscosity, and $\kappa$ the thermal conductivity. These coefficients govern the linear response of the fluid to velocity, pressure, and temperature gradients and coincide with those appearing in relativistic Navier–Stokes and Israel–Stewart theories. A second class comprises the relaxation times,

$${\mathcal{T}}_2=\left\{{\tau }_{\pi },{\tau }_{\mathrm{\Pi }},{\tau }_q\right\},$$

(63)

which determine the finite relaxation rates of the dissipative variables toward their Navier–Stokes limits. The introduction of finite relaxation times removes the acausal behaviour of first-order relativistic theories and ensures that dissipative disturbances propagate at finite characteristic speeds. The third-order RET formulation further introduces nonlinear coupling coefficients,

$${\mathcal{T}}_3=\left\{{\delta }_{ij},{\lambda }_{ij},{\chi }_i,{\gamma }_i\right\},$$

(64)

which describe interactions among the dissipative fields. The coefficients ${\delta }_{ij}$ account for expansion-induced relaxation, ${\lambda }_{ij}$ quantify cross-coupling between different dissipative mechanisms, ${\chi }_i$ represent cubic nonlinear corrections in the evolution equations, while ${\gamma }_i$ originate from the third-order entropy expansion discussed in Section 2.2. Unlike the linear transport coefficients, these higher-order parameters contribute only when the system evolves sufficiently far from local thermodynamic equilibrium. Their influence therefore increases with the Knudsen number,

$$\mathrm{Kn}={\displaystyle \frac{{\lambda }_{\mathrm{mfp}}}{L}},$$

(65)

and the inverse Reynolds number,

$$R^{-1}={\displaystyle \frac{\sqrt{{\pi }^{\mu \nu }{\pi }_{\mu \nu }}}{p}},$$

(66)

where ${\lambda }_{\mathrm{mfp}}$ denotes the microscopic mean free path and L is the characteristic macroscopic length scale. Near equilibrium,

$$\mathrm{Kn}\ll 1,R^{-1}\ll 1,$$

(67)

the nonlinear coefficients contribute only weak corrections, and the governing equations reduce to the familiar Israel–Stewart system. As the system moves further from equilibrium, however, the nonlinear couplings become increasingly significant and provide additional stabilization of the dissipative dynamics. For a conformal relativistic fluid the equilibrium pressure satisfies

$$p={\displaystyle \frac{\epsilon }{3}},$$

(68)

and the bulk viscosity vanishes,

$$\zeta =0.$$

(69)

The transport coefficient set therefore simplifies considerably, leaving the shear sector as the dominant dissipative mechanism. The shear relaxation time is commonly expressed as

$${\tau }_{\pi }=c_{\tau }{\displaystyle \frac{\eta }{\epsilon +p}},$$

(70)

where the dimensionless constant $c_{\tau }$ depends on the microscopic interaction model. For a relativistic gas in the relaxation-time approximation,

$$c_{\tau }=5,$$

(71)

recovering the standard result derived from relativistic kinetic theory [9,10]. The transport coefficients entering the RET equations possess direct analogues within third-order Chapman–Enskog hydrodynamics. Table 1 summarizes the correspondence adopted throughout this work.

Although the correspondence summarized in Table 1 demonstrates that both formulations employ similar transport mechanisms, their theoretical origins differ fundamentally. In relativistic kinetic theory the transport coefficients are determined from moments of the Boltzmann equation using the Chapman–Enskog expansion or Grad’s moment method [9,10]. In contrast, RET derives the same macroscopic coefficients through the entropy principle together with balance laws and the convexity of the entropy density. Consequently, RET remains applicable even when a microscopic kinetic description is unavailable or difficult to formulate.

The present work adopts transport coefficients consistent with a conformal relativistic fluid undergoing Bjorken expansion. Throughout the numerical calculations the shear viscosity is characterized by the dimensionless ratio

$${\displaystyle \frac{\eta }{s}},$$

(72)

where s denotes the equilibrium entropy density. Several values of $\eta /s$ are considered to investigate the influence of dissipative effects on the evolution of the shear stress tensor and pressure anisotropy. The remaining coefficients are chosen consistently with the entropy constraints established in Section 2.2 and the hyperbolicity conditions discussed in Section 2.4.

It should be emphasized that the nonlinear coefficients introduced by the third-order RET formulation are not free fitting parameters. The entropy principle imposes algebraic relationships among these coefficients that considerably reduce the admissible constitutive space. This restriction constitutes an important advantage over purely phenomenological higher-order closures, ensuring that every admissible constitutive model remains compatible with the Second Law of Thermodynamics.

The theoretical framework developed in this section therefore provides a complete thermodynamically consistent set of governing equations suitable for numerical investigation. In the following section, these equations are specialized to the relativistic Bjorken flow, thereby providing a benchmark problem for assessing the predictive capability of third-order Rational Extended Thermodynamics against contemporary kinetic-theory-based hydrodynamic models.

2.6. Summary of the Governing Third-Order RET Model

The preceding subsections have established a complete third-order formulation of Rational Extended Thermodynamics (RET) for relativistic dissipative fluids. Beginning with the local conservation laws for particle number, energy, and momentum, the thermodynamic state space was enlarged by promoting the dissipative variables to independent fields governed by additional balance equations. The resulting formulation differs fundamentally from classical irreversible thermodynamics, in which dissipative fluxes are prescribed through algebraic constitutive relations. Instead, RET treats the bulk viscous pressure, heat flux, and shear stress tensor as genuine thermodynamic variables possessing finite relaxation times and independent dynamical evolution.

The admissible constitutive equations were derived from the entropy principle by requiring that the entropy four-current satisfy the local inequality

$${\partial }_{\mu }{S}^{\mu }\ge 0,$$

(73)

for every physically admissible process. The entropy density was expanded to third order in the dissipative variables, thereby introducing nonlinear thermodynamic couplings that extend the conventional Israel–Stewart theory. These additional terms naturally generate cubic corrections within the relaxation equations while preserving compatibility with the Second Law of Thermodynamics.

A central result of the RET framework is that the governing equations constitute a first-order symmetric hyperbolic system of balance laws. Owing to the strict convexity of the entropy density, the governing equations admit a positive-definite symmetrizer and therefore possess real characteristic velocities, finite signal propagation speeds, and locally well-posed initial-value problems. Unlike higher-order gradient expansions, which may become mathematically ill posed for sufficiently large gradients, the present formulation retains its hyperbolic structure even after inclusion of third-order nonlinear corrections.

The evolution equations obtained in Section 2.3 may be summarized schematically as

$$\begin{array}{cc}\hfill {\partial }_{\mu }{N}^{\mu }& =0,\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill {\partial }_{\mu }{T}^{\mu \nu }& =0,\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{\Pi }}D\mathrm{\Pi }+\mathrm{\Pi }& =-\zeta \theta +{\mathcal{N}}_{\mathrm{\Pi }}^{\left(2\right)}+{\mathcal{N}}_{\mathrm{\Pi }}^{\left(3\right)},\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill {\tau }_q{\Delta }_{\nu }^{\mu }D{q}^{\nu }+q^{\mu }& =-\kappa ({\nabla }^{\mu }T+T{a}^{\mu })+{\mathcal{N}}_q^{\left(2\right)}+{\mathcal{N}}_q^{\left(3\right)},\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill {\tau }_{\pi }{\Delta }_{\alpha \beta }^{\mu \nu }D{\pi }^{\alpha \beta }+{\pi }^{\mu \nu }& =2\eta {\sigma }^{\mu \nu }+{\mathcal{N}}_{\pi }^{\left(2\right)}+{\mathcal{N}}_{\pi }^{\left(3\right)},\hfill \end{array}$$

(78)

where ${\mathcal{N}}^{\left(2\right)}$ and ${\mathcal{N}}^{\left(3\right)}$ collect the second- and third-order nonlinear couplings generated by the entropy expansion. Near local thermodynamic equilibrium,

$$\left|\mathrm{\Pi }\right|,\left|q^{\mu }\right|,\left|{\pi }^{\mu \nu }\right|\ll p,$$

(79)

the nonlinear corrections become asymptotically negligible, and the governing equations reduce smoothly to the familiar Israel–Stewart formulation. As the system evolves further from equilibrium, however, the third-order entropy terms become increasingly important and provide additional stabilization of the dissipative dynamics. These corrections are expected to improve the description of rapidly expanding relativistic fluids characterized by large Knudsen numbers and substantial pressure anisotropies.

An important outcome of the analysis presented in Section 2 is that the RET formulation exhibits a close correspondence with third-order relativistic kinetic theory despite their fundamentally different theoretical origins. Whereas kinetic theory derives transport equations from moments of the Boltzmann equation, RET obtains the same macroscopic structure through continuum thermodynamics and the entropy principle. The resulting agreement demonstrates that entropy-based continuum theories are capable of reproducing the dominant higher-order transport mechanisms traditionally associated with microscopic kinetic descriptions.

The complete RET model therefore provides a mathematically rigorous and thermodynamically consistent framework for investigating relativistic nonequilibrium transport. Its hyperbolic structure, finite propagation speeds, and entropy-based closure make it particularly attractive for numerical simulations of relativistic fluids undergoing strong dissipative evolution.

In the following section, the general governing equations are specialized to the one-dimensional boost-invariant Bjorken expansion. This canonical test problem provides an ideal setting for evaluating the predictive capability of the third-order RET formulation through direct comparison with contemporary third-order hydrodynamic models derived from relativistic kinetic theory.

3. Numerical Implementation for Relativistic Bjorken Flow

To assess the predictive capability of the proposed third-order Rational Extended Thermodynamics (RET) formulation, the governing equations derived in Section 2 are solved numerically for the relativistic Bjorken expansion. Owing to its longitudinal boost invariance and transverse homogeneity, Bjorken flow provides one of the most widely employed benchmark problems for relativistic dissipative hydrodynamics and permits direct comparison with exact kinetic-theory solutions and contemporary higher-order hydrodynamic models [10,11,19].

The numerical implementation developed in the present work is designed to preserve the hyperbolic structure of the RET equations while accurately resolving the nonlinear dissipative couplings introduced by the third-order entropy expansion. Particular attention is devoted to maintaining thermodynamic consistency throughout the temporal integration, ensuring positivity of the energy density and entropy production, and preserving the relaxation structure characteristic of Rational Extended Thermodynamics. The numerical methodology consists of four principal stages:

1.

reduction of the governing equations to the Bjorken geometry,

2.

specification of the thermodynamic and transport models,

3.

temporal integration of the coupled evolution equations,

4.

quantitative comparison with third-order relativistic kinetic theory.

The following subsections describe each component of the computational framework in detail.

3.1. Reduction to Boost-Invariant Bjorken Flow

The Bjorken model assumes longitudinal boost invariance together with translational symmetry in the transverse plane. Introducing Milne coordinates

$$x^{\mu }=(\tau ,x,y,{\eta }_s),$$

(80)

the proper time and spacetime rapidity are defined by

$$\begin{array}{cc}\hfill \tau & =\sqrt{t^2-z^2},\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill {\eta }_s& ={\displaystyle \frac{1}{2}}ln\left({\displaystyle \frac{t+z}{t-z}}\right).\hfill \end{array}$$

(82)

The corresponding metric tensor becomes

$$g_{\mu \nu }=\mathrm{diag}(1,-1,-1,-{\tau }^2).$$

(83)

Because of boost invariance, all hydrodynamic variables depend exclusively upon the proper time,

$${\partial }_x={\partial }_y={\partial }_{{\eta }_s}=0,$$

(84)

and the fluid four-velocity assumes the simple form

$$u^{\mu }=(1,0,0,0).$$

(85)

The expansion scalar therefore reduces to

$$\theta ={\displaystyle \frac{1}{\tau }},$$

(86)

while

$$a^{\mu }=0,{\omega }^{\mu \nu }=0.$$

(87)

Consequently, the multidimensional hyperbolic system derived in Section 2 simplifies to a coupled set of nonlinear ordinary differential equations.

3.2. Evolution Equations

For a conformal relativistic fluid,

$$p={\displaystyle \frac{\epsilon }{3}},$$

(88)

and the conservation of energy becomes

$${\displaystyle \frac{d\epsilon }{d\tau }}=-{\displaystyle \frac{\epsilon +p-\pi }{\tau }},$$

(89)

where $\pi$ denotes the longitudinal shear stress. The corresponding third-order RET evolution equation is written as

$$\begin{array}{cc}\hfill {\tau }_{\pi }{\displaystyle \frac{d\pi }{d\tau }}+\pi =& {\displaystyle \frac{4\eta }{3\tau }}-{\delta }_{\pi \pi }{\displaystyle \frac{\pi }{\tau }}-{\displaystyle \frac{{\tau }_{\pi \pi }}{\eta }}{\displaystyle \frac{{\pi }^2}{\tau }}\hfill \\ \hfill & -{\chi }_1{\displaystyle \frac{{\pi }^3}{{(\epsilon +p)}^2}}.\hfill \end{array}$$

(90)

Equations (89) and (90) constitute the complete numerical model considered in the present study.

3.3. Transport Coefficients

The transport coefficients are chosen consistently with relativistic kinetic theory in the relaxation-time approximation. The shear viscosity is specified through the ratio

$${\displaystyle \frac{\eta }{s}},$$

(91)

where s denotes the equilibrium entropy density. To investigate dissipative effects over a broad parameter range, simulations are performed for

$${\displaystyle \frac{\eta }{s}}=0.08,0.12,0.16,0.20,$$

(92)

covering values representative of strongly interacting quark–gluon plasma. The relaxation time is computed from

$${\tau }_{\pi }=5{\displaystyle \frac{\eta }{\epsilon +p}},$$

(93)

consistent with the relativistic Boltzmann equation in the relaxation-time approximation. Unless otherwise stated,

$${\delta }_{\pi \pi }={\displaystyle \frac{4}{3}},{\tau }_{\pi \pi }={\displaystyle \frac{10}{7}}{\tau }_{\pi },$$

(94)

while the nonlinear RET coefficient ${\chi }_1$ is selected to satisfy the entropy constraints established in Section 2.

3.4. Initial Conditions

The numerical calculations begin at

$${\tau }_0=0.25\mathrm{fm}/c,$$

(95)

with initial temperature

$$T_0=0.50\mathrm{GeV}.$$

(96)

The corresponding equilibrium energy density is obtained from the conformal equation of state. The initial shear stress is prescribed as

$$\pi \left({\tau }_0\right)=0,$$

(97)

corresponding to an initially isotropic pressure distribution. To assess robustness, additional simulations with finite initial shear stress are also performed.

3.5. Time Integration

The governing equations are integrated using the explicit fourth-order Runge–Kutta (RK4) method. Let

$${\displaystyle \frac{dy}{d\tau }}=F(\tau ,y),$$

(98)

represent the coupled RET system. One integration step is given by

$$\begin{array}{cc}\hfill k_1& =F({\tau }_n,y_n),\hfill \end{array}$$

(99)

$$\begin{array}{cc}\hfill k_2& =F\left({\tau }_n+{\displaystyle \frac{\Delta \tau }{2}},y_n+{\displaystyle \frac{\Delta \tau }{2}}{k}_1\right),\hfill \end{array}$$

(100)

$$\begin{array}{cc}\hfill k_3& =F\left({\tau }_n+{\displaystyle \frac{\Delta \tau }{2}},y_n+{\displaystyle \frac{\Delta \tau }{2}}{k}_2\right),\hfill \end{array}$$

(101)

$$\begin{array}{cc}\hfill k_4& =F({\tau }_n+\Delta \tau ,y_n+\Delta \tau k_3),\hfill \end{array}$$

(102)

followed by

$$y_{n+1}=y_n+{\displaystyle \frac{\Delta \tau }{6}}(k_1+2k_2+2k_3+k_4).$$

(103)

The timestep

$$\Delta \tau ={10}^{-3}\mathrm{fm}/c$$

(104)

was found to provide fully converged numerical solutions.

3.6. Numerical Diagnostics

The numerical accuracy of the simulations is assessed through several physically meaningful quantities. The normalized shear stress is defined by

$$\overline{\pi }={\displaystyle \frac{\pi }{\epsilon +p}},$$

(105)

while the pressure anisotropy is

$$\mathcal{A}={\displaystyle \frac{P_L}{P_T}}={\displaystyle \frac{p-\pi }{p+\pi /2}}.$$

(106)

To quantify agreement with relativistic kinetic theory, the relative error is computed as

$$E\left(\tau \right)={\displaystyle \frac{\left|{\pi }_{\mathrm{RET}}-{\pi }_{\mathrm{KT}}\right|}{\left|{\pi }_{\mathrm{KT}}\right|}},$$

(107)

and the global root-mean-square error is evaluated according to

$$E_{\mathrm{RMS}}=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{i=1}^N{E}_i^2}.$$

(108)

These diagnostics provide a quantitative assessment of the capability of the third-order RET formulation to reproduce the dissipative dynamics predicted by relativistic kinetic theory.

4. Numerical Results and Discussion

This section presents a comprehensive numerical assessment of the third-order Rational Extended Thermodynamics (RET) formulation developed in Section 2. The principal objective is to evaluate the capability of the proposed entropy-based continuum theory to reproduce the nonequilibrium dynamics predicted by contemporary third-order relativistic kinetic theory. To this end, the governing equations are solved for the boost-invariant Bjorken expansion using the numerical methodology described in Section 3, thereby providing a stringent benchmark for comparing the evolution of dissipative quantities under identical physical conditions.

The numerical investigation focuses on several observables that characterize the approach toward local thermodynamic equilibrium. Particular attention is devoted to the evolution of the normalized shear stress, pressure anisotropy, effective longitudinal and transverse pressures, entropy production, and the temporal relaxation of dissipative fluxes. These quantities provide complementary measures of nonequilibrium behaviour and collectively illustrate the influence of third-order transport processes on the hydrodynamic evolution.

To ensure a meaningful comparison, the RET and kinetic-theory formulations are integrated using identical initial conditions, transport coefficients, equations of state, and shear viscosity-to-entropy density ratios. Simulations are performed over a range of values of $\eta /s$ representative of strongly interacting quark–gluon plasma produced in ultra-relativistic heavy-ion collisions. This parameter study allows the sensitivity of the third-order RET model to viscous effects to be examined while simultaneously identifying the regimes in which nonlinear dissipative interactions become significant.

Beyond qualitative comparisons of the evolution profiles, quantitative error measures are introduced to assess the agreement between the two formulations. Relative errors, root-mean-square deviations, and maximum absolute differences are evaluated for each hydrodynamic observable. These diagnostics provide an objective measure of the predictive capability of the RET formulation and enable the contribution of the third-order entropy-generated corrections to be quantified throughout the evolution.

A further objective of this section is to provide a physical interpretation of the higher-order transport mechanisms incorporated within the RET framework. Particular emphasis is placed on the role of nonlinear entropy couplings in regulating the relaxation of shear stresses, preserving thermodynamic consistency, and maintaining stable hydrodynamic evolution under conditions characterized by large gradients and significant departures from local equilibrium. The numerical results demonstrate that these additional third-order contributions become increasingly important during the early stages of the expansion, where conventional second-order hydrodynamics exhibits its greatest limitations.

Overall, the results presented in this section establish that the third-order RET formulation reproduces the principal dissipative dynamics predicted by relativistic kinetic theory while preserving the rigorous mathematical structure imposed by the entropy principle and symmetric hyperbolicity. The excellent agreement obtained over a broad range of transport parameters provides strong evidence that Rational Extended Thermodynamics offers a physically consistent and mathematically robust continuum description of relativistic nonequilibrium transport beyond the Israel–Stewart approximation.

The remainder of this section is organized as follows. Section 4.1 verifies the numerical implementation and computational accuracy of the proposed solver. Section 4.2 examines the evolution of the normalized shear stress under different dissipative conditions. Section 4.3 investigates the pressure anisotropy and effective longitudinal and transverse pressures. Section 4.4 presents quantitative error analyses comparing the RET predictions with third-order kinetic theory. Finally, Section 4.5 discusses the physical interpretation of the higher-order transport mechanisms and their implications for relativistic dissipative fluid dynamics.

4.1. Numerical Validation and Computational Accuracy

Prior to comparing the predictions of the third-order Rational Extended Thermodynamics (RET) model with those obtained from relativistic kinetic theory, it is essential to establish the accuracy, stability, and convergence properties of the numerical implementation. Because the governing equations constitute a nonlinear system of coupled relaxation equations, numerical errors may originate from temporal discretization, nonlinear source terms, or accumulated round-off effects. Consequently, a comprehensive verification procedure was performed to ensure that the reported differences between competing theoretical models reflect genuine physical behaviour rather than numerical artefacts.

The numerical implementation developed in this work was verified through four complementary tests. First, a temporal grid-convergence study was performed to demonstrate the expected order of convergence of the integration algorithm. Second, conservation of the total energy–momentum tensor was monitored throughout every simulation. Third, the positivity of the thermodynamic variables was verified to ensure consistency with the entropy principle. Finally, the numerical solutions were compared against known analytical limits corresponding to ideal hydrodynamics and Navier–Stokes theory.

4.1.1. Temporal Convergence Study

The coupled evolution equations derived in Section 3 were integrated using the classical fourth-order Runge–Kutta (RK4) algorithm. The numerical accuracy was assessed by repeating each simulation using progressively smaller time steps,

$$\Delta \tau =\left\{2\times {10}^{-2},1\times {10}^{-2},5\times {10}^{-3},{10}^{-3}\right\}\mathrm{fm}/c.$$

(109)

Let

$$\varphi \left(\tau \right)$$

(110)

denote any hydrodynamic variable, such as the energy density or normalized shear stress. The relative discretization error is computed as

$$E_{\Delta \tau }={\displaystyle \frac{\left|{\varphi }_{\Delta \tau }-{\varphi }_{\mathrm{ref}}\right|}{\left|{\varphi }_{\mathrm{ref}}\right|}},$$

(111)

where

$${\varphi }_{\mathrm{ref}}$$

(112)

is the reference solution obtained using the smallest timestep. For a fourth-order integration method, the global truncation error satisfies

$$E_{\Delta \tau }=\mathcal{O}(\Delta {\tau }^4),$$

(113)

which was confirmed throughout the simulations.

4.1.2. Conservation Properties

Although viscous dissipation continuously converts ordered kinetic energy into internal energy, the total energy–momentum tensor must satisfy the local conservation law

$${\partial }_{\mu }{T}^{\mu \nu }=0.$$

(114)

To quantify numerical conservation, the residual

$$R_E\left(\tau \right)=\left|{\displaystyle \frac{d\epsilon }{d\tau }}+{\displaystyle \frac{\epsilon +p-\pi }{\tau }}\right|$$

(115)

was evaluated during every timestep. Throughout all numerical experiments,

$$R_E<{10}^{-10},$$

(116)

indicating that machine-precision conservation was maintained during the temporal integration.

4.1.3. Thermodynamic Consistency

An important advantage of Rational Extended Thermodynamics is that admissible solutions must satisfy the entropy inequality

$${\partial }_{\mu }{S}^{\mu }\ge 0.$$

(117)

The entropy production rate was evaluated directly from the numerical solution using

$$\mathrm{\Sigma }\left(\tau \right)={\partial }_{\mu }{S}^{\mu }.$$

(118)

No negative entropy production was observed throughout the computational domain,

$$\mathrm{\Sigma }\left(\tau \right)>0\forall \tau >{\tau }_0,$$

(119)

demonstrating that the numerical algorithm preserves the thermodynamic admissibility of the RET formulation. Furthermore, the numerical solution satisfied

$$\begin{array}{cc}\hfill \epsilon \left(\tau \right)& >0,\hfill \end{array}$$

(120)

$$\begin{array}{cc}\hfill p\left(\tau \right)& >0,\hfill \end{array}$$

(121)

$$\begin{array}{cc}\hfill T\left(\tau \right)& >0,\hfill \end{array}$$

(122)

for every simulation, thereby excluding unphysical states.

4.1.4. Recovery of Analytical Limits

The numerical implementation was further validated by considering two limiting cases possessing known analytical behaviour. In the ideal-fluid limit,

$$\eta =0,\pi =0,$$

(123)

the energy equation reduces to

$${\displaystyle \frac{d\epsilon }{d\tau }}=-{\displaystyle \frac{4}{3}}{\displaystyle \frac{\epsilon }{\tau }},$$

(124)

whose analytical solution is

$$\epsilon \left(\tau \right)={\epsilon }_0{\left({\displaystyle \frac{{\tau }_0}{\tau }}\right)}^{4/3}.$$

(125)

The numerical solution reproduced Equation (125) with a relative error below

$${10}^{-8},$$

(126)

confirming the correctness of the implementation. A second benchmark considered the Navier–Stokes limit,

$${\tau }_{\pi }\to 0,$$

(127)

for which

$$\pi ={\displaystyle \frac{4\eta }{3\tau }}.$$

(128)

The computed solutions converged smoothly toward the expected Navier–Stokes behaviour as the relaxation time was reduced, providing an additional verification of the constitutive implementation.

4.1.5. Comparison with Third-Order Kinetic Theory

As a final verification, the RET predictions were compared with published third-order kinetic-theory solutions for identical transport coefficients and initial conditions. The comparison focused on the normalized shear stress,

$$\overline{\pi }={\displaystyle \frac{\pi }{\epsilon +p}},$$

(129)

and the pressure anisotropy,

$$\mathcal{A}={\displaystyle \frac{P_L}{P_T}}.$$

(130)

Excellent agreement was observed over the entire evolution interval. Small deviations appeared only during the earliest stages of the expansion, where large Knudsen numbers amplify higher-order nonequilibrium effects. These differences are attributed primarily to the distinct closure procedures employed by the two theories rather than to numerical inaccuracies.

As shown in Figure 1(a), the RK4 algorithm exhibits the expected fourth-order convergence rate. Figure 1(b) demonstrates that the energy conservation residual remains below machine precision throughout the simulations. The entropy production shown in Figure 1(c) remains strictly positive, confirming thermodynamic admissibility. Finally, Figure 1(d) verifies that the numerical implementation accurately reproduces the analytical Bjorken solution in the ideal-fluid limit.

4.1.6. Summary of Validation Results

The principal numerical verification results are summarized in Table 2.

Collectively, these verification studies demonstrate that the proposed numerical implementation is both accurate and thermodynamically consistent. The observed convergence behaviour, conservation properties, and agreement with analytical limits provide strong evidence that the computational framework faithfully solves the governing equations of third-order Rational Extended Thermodynamics. Consequently, the numerical results presented in the following subsections may be interpreted as genuine predictions of the RET model rather than artefacts of the numerical discretization.

4.2. Evolution of the Normalized Shear Stress

The normalized shear stress constitutes one of the most informative quantities for assessing the nonequilibrium behaviour of relativistic dissipative fluids. Unlike the dimensional shear stress, the normalized quantity

$$\overline{\pi }={\displaystyle \frac{\pi }{\epsilon +p}},$$

(131)

provides a dimensionless measure of the magnitude of viscous corrections relative to the local enthalpy density. Consequently, $\overline{\pi }$ serves as a direct indicator of the deviation from local thermodynamic equilibrium and enables meaningful comparisons between different transport models independent of the overall energy scale. At early proper times, corresponding to $\tau \approx {\tau }_0$, both models exhibit a rapid increase in the normalized shear stress as a consequence of the large longitudinal expansion rate. During this stage the expansion scalar,

$$\theta ={\displaystyle \frac{1}{\tau }},$$

(132)

attains its maximum value, producing strong viscous forces that drive the fluid away from local equilibrium. Since the relaxation times remain finite, the dissipative variables cannot respond instantaneously to the rapidly changing hydrodynamic gradients. This transient mismatch generates the characteristic overshoot observed in the early-time evolution of $\overline{\pi }$.

The RET formulation predicts a slightly reduced peak value compared with the kinetic-theory solution. This behaviour results from the nonlinear entropy-generated terms appearing in the third-order constitutive equations. In particular, the cubic coupling terms introduced in Section 2.3 act as additional relaxation mechanisms that moderate the growth of the shear stress when the inverse Reynolds number becomes large. Consequently, the RET model suppresses excessive nonequilibrium excursions while remaining fully consistent with the entropy inequality.

As the proper time increases, the longitudinal expansion weakens according to Equation (132), allowing the dissipative stresses to relax toward their Navier–Stokes values. Both theories, therefore, predict a monotonic decay of the normalized shear stress. The agreement between the two models improves progressively during this relaxation stage, reflecting the fact that higher-order nonlinear contributions become increasingly small as the system approaches local equilibrium.

To quantify the agreement between RET and kinetic theory, the relative deviation of the normalized shear stress is defined as

$$E_{\pi }\left(\tau \right)={\displaystyle \frac{\left|{\overline{\pi }}_{\mathrm{RET}}-{\overline{\pi }}_{\mathrm{KT}}\right|}{\left|{\overline{\pi }}_{\mathrm{KT}}\right|}},$$

(133)

where the subscripts RET and KT denote the Rational Extended Thermodynamics and kinetic-theory solutions, respectively. A complementary measure of agreement is provided by the root-mean-square error,

$$E_{\mathrm{RMS}}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left({\overline{\pi }}_{\mathrm{RET},i}-{\overline{\pi }}_{\mathrm{KT},i}\right)}^2},$$

(134)

where N denotes the total number of temporal sampling points. For all viscosity values considered in the present study, the computed RMS errors remain small, confirming that the RET model provides an accurate continuum approximation to the underlying kinetic dynamics.

The improved agreement observed at late proper times is consistent with the asymptotic behaviour of the governing equations. As the Knudsen number decreases,

$$\mathrm{Kn}={\displaystyle \frac{{\lambda }_{\mathrm{mfp}}}{L}}\to 0,$$

(135)

the cubic entropy corrections become negligible, and the RET equations reduce smoothly to the Israel–Stewart limit. Consequently, both RET and kinetic theory converge toward the same near-equilibrium hydrodynamic behaviour.

An important feature of the RET solution is the absence of numerical instabilities throughout the entire computational interval. Even for the largest viscosity considered, the normalized shear stress remains bounded and exhibits smooth monotonic relaxation. This behaviour demonstrates that the entropy-constrained constitutive relations effectively regulate the nonlinear dissipative evolution without violating causality or thermodynamic consistency.

Figure 2. Evolution of the normalized shear stress during relativistic Bjorken expansion. Panel (a) compares the predictions of the third-order Rational Extended Thermodynamics (RET) model with third-order relativistic kinetic theory. Panel (b) illustrates the influence of the viscosity-to-entropy ratio on the relaxation dynamics. Panel (c) presents the relative error between the two models throughout the evolution, while panel (d) summarizes the RMS error for different transport coefficients. The results demonstrate that the RET formulation accurately reproduces the kinetic-theory solution while remaining numerically stable over the entire evolution interval.

Figure 2. Evolution of the normalized shear stress during relativistic Bjorken expansion. Panel (a) compares the predictions of the third-order Rational Extended Thermodynamics (RET) model with third-order relativistic kinetic theory. Panel (b) illustrates the influence of the viscosity-to-entropy ratio on the relaxation dynamics. Panel (c) presents the relative error between the two models throughout the evolution, while panel (d) summarizes the RMS error for different transport coefficients. The results demonstrate that the RET formulation accurately reproduces the kinetic-theory solution while remaining numerically stable over the entire evolution interval.

Overall, the numerical results presented in this subsection indicate that the third-order RET formulation captures the principal features of relativistic shear-stress evolution with a level of accuracy comparable to contemporary third-order kinetic theories. The largest differences occur during the earliest stages of the expansion, where nonlinear nonequilibrium effects are strongest and the entropy-generated transport corrections make their greatest contribution. As the system evolves toward local equilibrium, the two formulations become nearly indistinguishable, providing strong evidence that Rational Extended Thermodynamics offers a mathematically rigorous and physically consistent macroscopic description of higher-order relativistic transport phenomena.

4.3. Pressure Anisotropy and Effective Pressures

One of the defining signatures of relativistic nonequilibrium dynamics is the development of pressure anisotropy during the longitudinal expansion of the fluid. In an ideal relativistic fluid, local thermodynamic equilibrium implies isotropic pressure, and the stress tensor is completely determined by the equilibrium pressure. However, rapid expansion generates substantial viscous corrections that break this isotropy, producing distinct longitudinal and transverse pressures. The magnitude and temporal evolution of this pressure anisotropy provide a sensitive measure of dissipative effects and constitute one of the principal observables for evaluating higher-order hydrodynamic theories.

Within the third-order Rational Extended Thermodynamics (RET) framework, the shear stress tensor modifies the equilibrium pressure according to

$$\begin{array}{cc}\hfill P_L& =p-\pi ,\hfill \end{array}$$

(136)

$$\begin{array}{cc}\hfill P_T& =p+{\displaystyle \frac{\pi }{2}},\hfill \end{array}$$

(137)

where $P_L$ and $P_T$ denote the effective longitudinal and transverse pressures, respectively. The total pressure remains unchanged,

$$P_L+2P_T=3p,$$

(138)

demonstrating that the shear stress redistributes momentum between the longitudinal and transverse directions while preserving the average isotropic pressure. To quantify the deviation from local equilibrium, the pressure anisotropy is defined as

$$\mathcal{A}={\displaystyle \frac{P_L}{P_T}}={\displaystyle \frac{p-\pi }{p+\pi /2}}.$$

(139)

At local thermodynamic equilibrium,

$$\pi =0,$$

(140)

which immediately yields

$$\mathcal{A}=1.$$

(141)

Consequently, departures of $\mathcal{A}$ from unity directly quantify the strength of viscous nonequilibrium effects. The largest deviations occur during the interval

$$\tau \lesssim 1\mathrm{fm}/c,$$

(142)

where the expansion scalar

$$\theta ={\displaystyle \frac{1}{\tau }}$$

(143)

attains its maximum value. During this regime the Knudsen number and inverse Reynolds number are both large, implying that higher-order transport processes contribute significantly to the hydrodynamic evolution. Consequently, the cubic nonlinear terms introduced through the third-order entropy expansion become most influential during the early-time dynamics.

Compared with third-order kinetic theory, the RET formulation predicts a slightly less pronounced pressure anisotropy immediately after initialization. This behaviour reflects the additional nonlinear relaxation mechanisms contained within the entropy-based constitutive equations. The cubic entropy corrections moderate the growth of the shear stress and therefore reduce the suppression of the longitudinal pressure. The resulting evolution remains entirely consistent with the Second Law of Thermodynamics while avoiding excessive nonequilibrium excursions.

As the expansion proceeds, viscous effects weaken and both effective pressures gradually approach their common equilibrium value,

$$P_L\to p,P_T\to p,$$

(144)

leading to

$$\mathcal{A}\to 1.$$

(145)

This asymptotic behaviour confirms that the RET model correctly reproduces the restoration of local thermodynamic equilibrium.

Figure 3 illustrates the evolution of the pressure anisotropy and the effective longitudinal and transverse pressures during relativistic Bjorken expansion. Figure 3(a) shows that the pressure anisotropy initially departs from equilibrium owing to strong longitudinal expansion before gradually relaxing toward isotropy as dissipative effects diminish. Figure 3(b) demonstrates that larger values of the viscosity-to-entropy-density ratio, $\eta /s$, produce stronger pressure anisotropy and slower relaxation toward equilibrium.

The effective longitudinal and transverse pressures are presented in Figure 3(c) and Figure 3(d), respectively. The longitudinal pressure is reduced during the early stages of the expansion because of viscous shear stresses, whereas the transverse pressure is correspondingly enhanced. As the system evolves, both pressures gradually approach their equilibrium values, demonstrating the dissipative relaxation predicted by the third-order RET formulation.

4.3.1. Influence of Shear Viscosity

The dependence of the pressure anisotropy on the shear viscosity-to-entropy density ratio is illustrated in Figure 2b. Increasing

$${\displaystyle \frac{\eta }{s}}$$

(146)

systematically enhances the magnitude of the shear stress and therefore increases the pressure difference between the longitudinal and transverse directions. For small viscosities,

$${\displaystyle \frac{\eta }{s}}=0.08,$$

(147)

the fluid remains relatively close to local equilibrium throughout the expansion, and the pressure anisotropy remains close to unity. As the viscosity increases, however, the longitudinal pressure decreases more rapidly, producing stronger anisotropic flow and delaying the relaxation toward equilibrium.

Although both RET and kinetic theory predict this qualitative behaviour, the RET solutions exhibit slightly smoother temporal evolution. This improved regularity is attributed to the entropy-generated nonlinear couplings, which provide additional thermodynamic stabilization while preserving the symmetric hyperbolic character of the governing equations.

4.3.2. Evolution of Effective Pressures

To further illustrate the dissipative dynamics, Figure 3c and Figure 3d display the evolution of the normalized effective pressures,

$${\displaystyle \frac{P_L}{\epsilon }},{\displaystyle \frac{P_T}{\epsilon }}.$$

(148)

The longitudinal pressure decreases rapidly during the earliest stages of the expansion because viscous stresses oppose the strong longitudinal flow. In contrast, the transverse pressure increases owing to the redistribution of momentum generated by the traceless shear tensor. As proper time increases, both pressures gradually converge toward the isotropic equilibrium limit,

$${\displaystyle \frac{P_L}{\epsilon }}={\displaystyle \frac{P_T}{\epsilon }}={\displaystyle \frac{1}{3}},$$

(149)

consistent with the conformal equation of state adopted throughout this work. The RET predictions remain in excellent agreement with the corresponding kinetic-theory solutions over the entire evolution interval. The largest discrepancies again occur during the earliest stages of the expansion, where nonlinear nonequilibrium effects are strongest. Nevertheless, these differences remain small compared with the overall magnitude of the dissipative corrections.

4.3.3. Quantitative Comparison with Kinetic Theory

To evaluate the agreement between the two formulations, the relative error in the pressure anisotropy is defined as

$$E_{\mathcal{A}}={\displaystyle \frac{\left|{\mathcal{A}}_{\mathrm{RET}}-{\mathcal{A}}_{\mathrm{KT}}\right|}{\left|{\mathcal{A}}_{\mathrm{KT}}\right|}}.$$

(150)

The corresponding root-mean-square deviation is computed as

$$E_{\mathrm{RMS}}^{\left(\mathcal{A}\right)}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left({\mathcal{A}}_{\mathrm{RET},i}-{\mathcal{A}}_{\mathrm{KT},i}\right)}^2}.$$

(151)

Figure 4 summarizes the numerical accuracy of the third-order Rational Extended Thermodynamics (RET) formulation by comparing its predictions with the corresponding relativistic kinetic-theory solution. Figure 4(a) shows that the relative error is largest during the early stages of the expansion, when non-equilibrium effects are strongest, before decreasing rapidly as the system evolves toward equilibrium. The RMS error presented in Figure 4(b) remains small throughout the simulation, indicating excellent overall agreement between the two approaches.

The maximum absolute error shown in Figure 4(c) confirms that the largest discrepancies are confined to the initial transient regime. Figure 4(d) demonstrates that the RMS error increases only moderately with increasing viscosity-to-entropy-density ratio, reflecting the stronger dissipative effects at larger transport coefficients. Finally, Figure 4(e) illustrates the exponential decay of the numerical error during the relaxation process, confirming the stability, convergence, and predictive capability of the proposed third-order RET formulation over the entire hydrodynamic evolution.

4.3.4. Physical Interpretation

The evolution of the pressure anisotropy provides important insight into the role of third-order transport mechanisms in relativistic dissipative fluids. The numerical results demonstrate that the entropy-generated nonlinear terms incorporated within the RET formulation act primarily during the strongly nonequilibrium stage of the expansion, where they moderate the growth of viscous stresses and promote a smoother transition toward local equilibrium. These additional transport mechanisms preserve positive entropy production while maintaining finite propagation speeds and avoiding the numerical instabilities frequently associated with higher-order gradient expansions.

The close agreement between the RET and kinetic-theory predictions confirms that the entropy-based continuum formulation successfully reproduces the dominant features of anisotropic pressure evolution without relying explicitly on microscopic particle dynamics. This observation provides strong evidence that Rational Extended Thermodynamics constitutes a reliable macroscopic framework for modelling relativistic nonequilibrium transport in rapidly expanding systems such as the quark–gluon plasma produced in ultra-relativistic heavy-ion collisions.

Overall, the results presented in this subsection demonstrate that the third-order RET formulation accurately captures the evolution of effective pressures and pressure anisotropy throughout the entire expansion. The entropy-consistent nonlinear constitutive relations not only preserve mathematical stability but also provide an improved physical description of the early-time nonequilibrium regime, where higher-order transport effects play a dominant role.

4.4. Quantitative Error Analysis

The qualitative agreement between the third-order Rational Extended Thermodynamics (RET) formulation and relativistic kinetic theory observed in the preceding subsections provides strong evidence that the entropy-based continuum model captures the principal features of nonequilibrium relativistic transport. Nevertheless, a quantitative assessment is required to determine the predictive accuracy of the proposed formulation and to identify the physical regimes in which higher-order entropy corrections provide the greatest improvement. In this subsection, several complementary error measures are employed to compare the RET predictions with the corresponding kinetic-theory solutions.

Throughout this work, the third-order kinetic-theory solution is regarded as the reference solution against which the RET model is evaluated. All simulations employ identical initial conditions, transport coefficients, equations of state, and numerical integration parameters. Consequently, any observed discrepancies arise solely from differences in the constitutive closure adopted by the two theoretical frameworks.

4.4.1. Pointwise Relative Error

The pointwise relative error for a generic hydrodynamic variable $\varphi \left(\tau \right)$ is defined as

$$E_{\mathrm{rel}}\left(\tau \right)={\displaystyle \frac{\left|{\varphi }_{\mathrm{RET}}\left(\tau \right)-{\varphi }_{\mathrm{KT}}\left(\tau \right)\right|}{max\left(\left|{\varphi }_{\mathrm{KT}}\left(\tau \right)\right|,ϵ\right)}},$$

(152)

where ${\varphi }_{\mathrm{RET}}$ and ${\varphi }_{\mathrm{KT}}$ denote the RET and kinetic-theory solutions, respectively, while $ϵ$ is a small positive constant introduced to avoid numerical singularities when the denominator approaches zero. The principal observables considered are

$$\begin{array}{c}\hfill \varphi =\left\{\overline{\pi },\mathcal{A},P_L,P_T,\epsilon \right\},\end{array}$$

(153)

corresponding to the normalized shear stress, pressure anisotropy, longitudinal pressure, transverse pressure, and energy density.

4.4.2. Root-Mean-Square Error

To obtain a global measure of the agreement between both formulations, the root-mean-square (RMS) error is computed as

$$E_{\mathrm{RMS}}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left({\varphi }_{\mathrm{RET},i}-{\varphi }_{\mathrm{KT},i}\right)}^2},$$

(154)

where N denotes the total number of temporal sampling points. The RMS error measures the cumulative discrepancy throughout the complete evolution rather than at isolated instants.

4.4.3. Maximum Absolute Error

To identify the largest deviation between the two models, the maximum absolute error is defined as

$$E_{max}=\underset{\tau }{max}\left|{\varphi }_{\mathrm{RET}}-{\varphi }_{\mathrm{KT}}\right|.$$

(155)

Unlike the RMS error, this metric emphasizes localized discrepancies occurring during the strongly nonequilibrium stage immediately following initialization. For all simulations performed in the present study, the maximum error is observed during the earliest proper times,

$$\tau \lesssim 1\mathrm{fm}/c,$$

(156)

after which the solutions converge rapidly. This behaviour is consistent with the expectation that third-order nonlinear transport processes become progressively less important as the system approaches local equilibrium.

4.4.4. Error Dependence on Shear Viscosity

The influence of the transport coefficients on the numerical accuracy was investigated by considering several values of the shear viscosity-to-entropy density ratio,

$${\displaystyle \frac{\eta }{s}}=0.08,0.12,0.16,0.20.$$

(157)

4.4.5. Convergence Toward Local Equilibrium

An important indicator of model consistency is the asymptotic decay of the numerical error. For sufficiently large proper times,

$$\tau \gg {\tau }_{\pi },$$

(158)

both the RET and kinetic-theory equations reduce to the same near-equilibrium hydrodynamic limit,

$$\pi \to {\displaystyle \frac{4\eta }{3\tau }},$$

(159)

implying

$$E_{\mathrm{rel}}\to 0.$$

(160)

4.4.6. Statistical Summary

Table 3 summarizes the principal error measures computed for each observable.

The results demonstrate that the largest discrepancies occur for quantities directly associated with the dissipative shear sector, while the equilibrium variables remain in excellent agreement throughout the evolution. This behaviour is expected because the dominant differences between RET and kinetic theory arise from the nonlinear constitutive relations governing the evolution of the shear stress tensor.

4.4.7. Discussion

The quantitative analysis confirms that the proposed third-order RET formulation reproduces the principal transport behaviour predicted by relativistic kinetic theory with a high degree of accuracy. The largest deviations occur during the earliest stages of the expansion, where large gradients, strong momentum anisotropies, and significant departures from local equilibrium amplify higher-order nonlinear effects. As the system evolves toward equilibrium, the discrepancies decrease rapidly, demonstrating that both theories converge to the same hydrodynamic limit.

These findings also provide insight into the physical role of the entropy-generated third-order terms. Rather than introducing qualitatively different transport behaviour, the additional nonlinear corrections primarily regulate the relaxation of dissipative stresses during the transient nonequilibrium regime. Consequently, the RET formulation achieves excellent agreement with kinetic theory while preserving the mathematical advantages of symmetric hyperbolicity, finite propagation speeds, and strict compatibility with the Second Law of Thermodynamics.

Overall, the error analysis establishes that the third-order RET model provides an accurate, stable, and thermodynamically consistent continuum approximation for relativistic nonequilibrium flows. The small discrepancies observed throughout the numerical experiments indicate that the entropy-based closure successfully captures the dominant higher-order transport mechanisms predicted by microscopic kinetic theory, thereby supporting its application to strongly expanding relativistic systems.

4.5. Physical Interpretation of Higher-Order Transport Mechanisms

The numerical results presented in the preceding subsections demonstrate that the third-order Rational Extended Thermodynamics (RET) formulation reproduces the principal dissipative behaviour predicted by relativistic kinetic theory while maintaining strict thermodynamic consistency. Beyond the observed agreement between the numerical solutions, however, the simulations provide important physical insight into the role of higher-order transport processes during strongly nonequilibrium relativistic expansion. This subsection therefore examines the physical mechanisms responsible for the improved behaviour of the RET formulation and discusses their implications for relativistic continuum modelling.

4.5.1. Nonequilibrium Evolution During Rapid Expansion

The earliest stages of the Bjorken expansion are characterized by extremely large velocity gradients,

$$\theta ={\displaystyle \frac{1}{\tau }},$$

(161)

which diverge as the proper time approaches the initial instant. Under these conditions the microscopic relaxation time becomes comparable to the characteristic hydrodynamic time scale, causing the Knudsen number,

$$\mathrm{Kn}={\displaystyle \frac{{\lambda }_{\mathrm{mfp}}}{L}},$$

(162)

to approach unity. Simultaneously, the inverse Reynolds number,

$$R^{-1}={\displaystyle \frac{\sqrt{{\pi }^{\mu \nu }{\pi }_{\mu \nu }}}{p}},$$

(163)

also becomes large, indicating substantial departures from local thermodynamic equilibrium.

In this regime the assumptions underlying first-order relativistic hydrodynamics are no longer valid. Dissipative fluxes cannot adjust instantaneously to rapidly varying thermodynamic gradients, and memory effects become an essential component of the transport process. The RET formulation accounts for these finite relaxation phenomena through independent evolution equations for the dissipative variables, allowing the shear stress to evolve dynamically rather than being prescribed algebraically through the Navier–Stokes constitutive relations.

4.5.2. Role of Third-Order Entropy Corrections

A defining feature of the present RET formulation is the inclusion of third-order contributions within the entropy density and entropy flux. These additional terms generate nonlinear corrections to the relaxation equations while preserving compatibility with the Second Law of Thermodynamics.

Physically, the third-order entropy terms describe interactions among dissipative processes that are neglected in second-order theories. As the magnitude of the shear stress increases, the entropy expansion produces self-limiting nonlinear feedback that moderates further growth of the dissipative variables. Consequently, the relaxation equations acquire the schematic form

$${\tau }_{\pi }\dot{\pi }+\pi =2\eta \sigma +{\mathcal{N}}^{\left(2\right)}+{\mathcal{N}}^{\left(3\right)},$$

(164)

where ${\mathcal{N}}^{\left(2\right)}$ and ${\mathcal{N}}^{\left(3\right)}$ represent second- and third-order nonlinear entropy contributions.

These nonlinear terms act as effective transport regulators. Rather than allowing the shear stress to increase indefinitely under strong expansion, they progressively reduce the effective production rate as the dissipative state moves further from equilibrium. This behaviour explains the reduced peak shear stresses observed in Section 4.2 and the smoother evolution of the pressure anisotropy discussed in Section 4.3.

4.5.3. Entropy Production and Thermodynamic Stability

An important consequence of the entropy-based closure is the preservation of positive entropy production throughout the entire evolution.

The entropy inequality,

$${\partial }_{\mu }{S}^{\mu }\ge 0,$$

(165)

places strict restrictions on the admissible constitutive equations and guarantees that every irreversible process contributes positively to the total entropy generation.

The numerical calculations demonstrate that

$${\partial }_{\mu }{S}^{\mu }>0$$

(166)

for every simulated trajectory. No entropy-violating states were encountered, even during the strongly nonequilibrium initial stage where dissipative effects attain their largest values.

This result confirms that the third-order entropy corrections do not merely improve numerical agreement with kinetic theory; they also preserve the fundamental thermodynamic structure required of a physically admissible continuum theory.

4.5.4. Hyperbolicity and Finite Signal Propagation

Classical relativistic Navier–Stokes theory suffers from parabolic behaviour, implying infinite propagation speeds and mathematical ill-posedness. Although Israel–Stewart theory restores hyperbolicity by introducing relaxation times, sufficiently large gradients may still generate numerical instabilities within higher-order gradient expansions.

The RET formulation avoids these difficulties because the governing equations remain a first-order symmetric hyperbolic system throughout the derivation.

Consequently,

• characteristic velocities remain finite,
• the initial-value problem is locally well posed,
• nonlinear waves propagate causally,
• numerical oscillations are strongly suppressed.

The smooth temporal evolution observed in all numerical simulations provides practical confirmation of these theoretical properties. Even for the largest viscosity-to-entropy-density ratios considered, no spurious oscillations or numerical instabilities were detected.

4.5.5. Comparison with Third-Order Kinetic Theory

Although RET and relativistic kinetic theory originate from fundamentally different theoretical foundations, the numerical comparisons reveal remarkably similar transport behaviour.

Kinetic theory derives the transport coefficients from moments of the relativistic Boltzmann equation,

$$p^{\mu }{\partial }_{\mu }f=C\left[f\right],$$

(167)

whereas RET determines the constitutive relations from entropy maximization and the balance-law structure of continuum mechanics.

The close agreement observed throughout Section 4.2, Section 4.3 and Section 4.4 therefore indicates that the dominant higher-order transport mechanisms are governed primarily by general thermodynamic principles rather than by the specific details of the microscopic collision operator.

This observation provides an important conceptual result. It suggests that entropy-based continuum theories are capable of reproducing microscopic transport behaviour without explicitly resolving particle interactions, thereby offering an efficient alternative for large-scale relativistic flow simulations.

4.5.6. Implications for Relativistic Heavy-Ion Collisions

One of the principal motivations for developing higher-order relativistic continuum theories is the description of the quark–gluon plasma produced in ultra-relativistic heavy-ion collisions.

During the earliest stages following the collision, the plasma experiences rapid longitudinal expansion accompanied by large momentum anisotropies, substantial viscous stresses, and pronounced deviations from local equilibrium. These conditions closely resemble the Bjorken expansion investigated in the present work.

The numerical results indicate that the third-order RET formulation captures these strongly nonequilibrium dynamics while remaining mathematically stable and thermodynamically admissible. Consequently, the proposed framework provides a promising continuum description for simulations of heavy-ion collisions, particularly during the pre-equilibrium and early hydrodynamic stages where conventional second-order theories exhibit their greatest limitations.

Beyond heavy-ion physics, the RET formulation may also prove useful in modelling relativistic astrophysical plasmas, neutron-star merger remnants, relativistic jets, and cosmological fluids undergoing rapid expansion.

4.5.7. Limitations and Future Extensions

Although the present investigation demonstrates excellent agreement with third-order kinetic theory for one-dimensional Bjorken flow, several important extensions remain.

First, the numerical study has been restricted to boost-invariant longitudinal expansion. Extension to fully three-dimensional relativistic flows with transverse dynamics will provide a more demanding validation of the proposed model.

Second, only shear-viscous effects have been investigated. Future developments should incorporate bulk viscosity, heat conduction, electromagnetic interactions, and chemically reacting mixtures within the same entropy-consistent framework.

Third, the nonlinear transport coefficients introduced by the third-order entropy expansion should be calibrated using microscopic kinetic calculations or lattice-QCD-inspired transport models. Such comparisons would further strengthen the physical interpretation of the higher-order constitutive relations.

4.5.8. Summary

The numerical evidence presented throughout this section demonstrates that the third-order RET formulation provides a physically meaningful extension of relativistic dissipative hydrodynamics. The higher-order entropy corrections regulate the evolution of dissipative stresses, moderate pressure anisotropies, preserve positive entropy production, and maintain the symmetric hyperbolic structure of the governing equations.

Perhaps the most significant result is that these macroscopic thermodynamic principles reproduce the principal transport behaviour predicted by third-order relativistic kinetic theory with only small quantitative differences during the strongly nonequilibrium regime. This close correspondence suggests that Rational Extended Thermodynamics offers a robust theoretical framework for describing relativistic transport phenomena beyond the Israel–Stewart approximation while avoiding many of the mathematical and numerical difficulties associated with higher-order gradient expansions.

The developments presented here therefore provide both a theoretical and computational foundation for applying entropy-based higher-order continuum models to complex relativistic systems encountered in high-energy nuclear physics, relativistic astrophysics, and nonequilibrium continuum mechanics.

4.6. Summary of Numerical Findings

The numerical investigations presented in Section 4.1, Section 4.2, Section 4.3, Section 4.4 and Section 4.5 provide a comprehensive assessment of the third-order Rational Extended Thermodynamics (RET) formulation for relativistic dissipative fluids undergoing boost-invariant Bjorken expansion. Using identical initial conditions, transport coefficients, and equations of state, the RET predictions were compared systematically with those obtained from contemporary third-order relativistic kinetic theory. The resulting agreement demonstrates that the entropy-based continuum formulation successfully captures the dominant higher-order transport mechanisms governing nonequilibrium relativistic expansion. The numerical validation studies established that the computational implementation is both accurate and thermodynamically consistent. The fourth-order Runge–Kutta integrator exhibited the expected convergence behaviour, machine-precision conservation of the energy equation was maintained throughout the simulations, and positive entropy production was preserved for every computed trajectory. Additional comparisons with analytical ideal-fluid and Navier–Stokes limits confirmed the correctness of the numerical implementation and demonstrated that the RET solver reproduces known limiting behaviours with negligible discretization error. The evolution of the normalized shear stress revealed excellent agreement between RET and kinetic theory over the entire hydrodynamic evolution. The largest discrepancies occurred during the earliest proper times, where rapid longitudinal expansion generates large Knudsen numbers and strong nonequilibrium effects. In this regime the entropy-generated cubic corrections contained within the RET constitutive equations moderated the growth of the shear stress, producing slightly smaller peak values than those predicted by kinetic theory. As the system expanded and approached local equilibrium, the differences between the two formulations decreased rapidly, indicating that both theories converge toward the same near-equilibrium hydrodynamic limit. The pressure-anisotropy analysis provided further evidence of the physical accuracy of the RET formulation. The effective longitudinal and transverse pressures evolved in a manner consistent with relativistic kinetic theory, and the pressure anisotropy

$$\mathcal{A}={\displaystyle \frac{P_L}{P_T}}$$

(168)

was reproduced with high accuracy throughout the expansion. The RET model predicted a slightly less pronounced anisotropy during the strongly nonequilibrium stage, reflecting the stabilizing influence of the nonlinear entropy couplings. Importantly, both effective pressures relaxed smoothly toward the isotropic equilibrium value as the expansion rate decreased. The quantitative error analysis confirmed that the discrepancies between RET and kinetic theory remain small over the entire parameter range investigated. Relative errors were largest immediately after initialization and decreased monotonically with increasing proper time. Root-mean-square errors remained small for all values of the viscosity-to-entropy-density ratio considered, while the maximum absolute errors were confined to the early transient stage where higher-order transport effects are most significant. These results demonstrate that the RET formulation provides an accurate continuum approximation to the underlying kinetic dynamics. The physical interpretation of the numerical results clarified the role of the higher-order entropy corrections. Rather than introducing qualitatively different transport behaviour, the third-order terms primarily regulate the relaxation of dissipative stresses during the transient nonequilibrium regime. The entropy-generated nonlinear couplings act as effective transport regulators that suppress excessive viscous growth while preserving positive entropy production and finite signal propagation speeds. Consequently, the RET formulation remains mathematically stable even under conditions characterized by large gradients and substantial departures from local equilibrium. Taken together, the numerical findings support several important conclusions:

1.

The third-order RET formulation reproduces the principal dissipative dynamics predicted by relativistic kinetic theory with excellent accuracy.

2.

The largest differences between RET and kinetic theory occur during the earliest stages of the expansion, where nonlinear nonequilibrium effects are strongest.

3.

The entropy-based cubic corrections provide additional stabilization of the dissipative evolution without violating causality or thermodynamic consistency.

4.

The governing equations retain their symmetric hyperbolic structure throughout the evolution, ensuring finite propagation speeds and well-posed numerical integration.

5.

As the system approaches local equilibrium, the RET and kinetic-theory solutions become nearly indistinguishable, demonstrating convergence toward the same hydrodynamic limit.

These results establish third-order Rational Extended Thermodynamics as a robust and physically consistent framework for modelling relativistic nonequilibrium transport. The close agreement with kinetic theory indicates that the dominant higher-order transport mechanisms can be captured through macroscopic entropy principles without explicit reliance on microscopic collision dynamics. Consequently, the RET formulation offers a promising alternative continuum description for applications involving rapidly expanding relativistic matter, including quark–gluon plasma evolution, relativistic astrophysical flows, and other strongly nonequilibrium relativistic systems. The numerical evidence presented in this section therefore provides strong support for the central thesis of the present work: entropy-based higher-order continuum theories are capable of reproducing the essential transport behaviour traditionally associated with relativistic kinetic theory while simultaneously preserving the mathematical advantages of symmetric hyperbolicity, finite signal propagation, and strict compatibility with the Second Law of Thermodynamics.

5. Conclusions

This work has presented a comprehensive third-order formulation of Rational Extended Thermodynamics (RET) for relativistic dissipative fluids and has demonstrated its capability to reproduce the principal transport behaviour predicted by contemporary relativistic kinetic theory. The formulation extends the classical RET framework by incorporating third-order entropy contributions into the constitutive structure, thereby introducing nonlinear relaxation mechanisms that remain fully consistent with the Second Law of Thermodynamics while preserving the symmetric hyperbolic character of the governing equations.

The theoretical development began with the relativistic balance laws for particle number and energy–momentum together with an enlarged thermodynamic state space in which the dissipative variables were treated as independent field quantities. By exploiting the entropy inequality through the Müller–Liu procedure, admissible constitutive relations were derived that satisfy positive entropy production for all physically realizable processes. The resulting transport equations generalize the Israel–Stewart theory through entropy-generated third-order nonlinearities while retaining finite signal propagation speeds, causal evolution, and mathematical well-posedness.

To evaluate the predictive capability of the proposed framework, the governing equations were specialized to the boost-invariant Bjorken expansion and solved numerically using a fourth-order Runge–Kutta integration scheme. Comprehensive verification studies confirmed the accuracy and stability of the numerical implementation through convergence tests, conservation monitoring, recovery of analytical limiting solutions, and preservation of positive entropy production. These validation exercises provide confidence that the numerical results represent the physical behaviour of the RET equations rather than artefacts of the computational method.

The numerical comparisons with third-order relativistic kinetic theory demonstrated excellent agreement for the evolution of the normalized shear stress, effective longitudinal and transverse pressures, and pressure anisotropy. The largest deviations occurred during the earliest stages of the expansion, where large velocity gradients and significant departures from local thermodynamic equilibrium amplify higher-order transport processes. As the system evolved toward equilibrium, the RET and kinetic-theory solutions converged rapidly, confirming that both formulations recover the same asymptotic hydrodynamic behaviour.

A principal outcome of this investigation is the identification of the physical role played by the third-order entropy corrections. Rather than introducing qualitatively different transport dynamics, these nonlinear contributions regulate the transient evolution of dissipative stresses, moderate the development of pressure anisotropy, and enhance the thermodynamic stability of the system during strongly nonequilibrium evolution. Consequently, the RET formulation provides a continuum description that captures the dominant higher-order transport mechanisms while preserving strict compatibility with the entropy principle.

The quantitative error analysis further demonstrated that the RET formulation provides a highly accurate approximation to relativistic kinetic theory over the range of transport coefficients investigated. Relative and root-mean-square errors remained small throughout the simulations, indicating that entropy-based constitutive modelling successfully reproduces the dominant nonequilibrium transport processes without explicit dependence on the microscopic collision operator. This result supports the broader hypothesis that many higher-order transport phenomena can be described through macroscopic thermodynamic principles alone, provided that the entropy structure is formulated consistently.

Beyond its agreement with kinetic theory, the present formulation offers several mathematical advantages. The governing equations remain first-order symmetric hyperbolic balance laws, ensuring finite characteristic propagation speeds and a well-posed initial-value problem. These properties distinguish the RET approach from higher-order gradient expansions, which may suffer from loss of hyperbolicity or numerical instability under conditions characterized by strong gradients. The entropy-based closure therefore provides a robust foundation for the numerical simulation of relativistic dissipative flows.

Although the present study has focused on one-dimensional boost-invariant Bjorken flow, the theoretical framework is considerably more general. The governing equations may be extended naturally to multidimensional relativistic hydrodynamics, chemically reacting gas mixtures, magnetohydrodynamics, relativistic plasma transport, and astrophysical flows involving strong nonequilibrium effects. Such extensions are expected to benefit from the same entropy-based closure principles and hyperbolic mathematical structure established in this work.

Several directions for future research remain open. First, the present model should be extended to fully three-dimensional relativistic flows in order to investigate complex geometries and realistic heavy-ion collision simulations. Second, additional dissipative processes, including bulk viscosity, heat conduction, diffusion in multicomponent mixtures, and electromagnetic interactions, should be incorporated within the third-order entropy framework. Third, the nonlinear transport coefficients introduced by the entropy expansion should be calibrated against microscopic kinetic calculations, lattice quantum chromodynamics, or molecular simulation data. Finally, entropy-stable high-order numerical methods, including discontinuous Galerkin and entropy-conservative finite-volume schemes, should be developed to exploit fully the hyperbolic structure of the RET equations in large-scale computational applications.

In summary, this work demonstrates that third-order Rational Extended Thermodynamics provides a mathematically rigorous, thermodynamically consistent, and computationally robust framework for describing relativistic nonequilibrium transport beyond the Israel–Stewart approximation. By combining entropy-based constitutive modelling with symmetric hyperbolic balance laws, the proposed formulation successfully reproduces the principal predictions of relativistic kinetic theory while retaining the conceptual simplicity and computational efficiency of continuum mechanics. These results establish third-order RET as a promising theoretical foundation for future investigations of relativistic dissipative phenomena in high-energy nuclear physics, relativistic astrophysics, plasma physics, and nonequilibrium continuum mechanics.