Relativistic Heat Conduction in the Large-Flux Regime

Abstract

We propose a general procedure for evaluating, directly from microphysics, the constitutive relations of heat-conducting fluids in regimes of large fluxes of heat. Our choice of hydrodynamic formalism is Carter’s two-fluid theory, which happens to coincide with Öttinger’s GENERIC theory for relativistic heat conduction. This is a natural framework, as it should correctly describe the relativistic “inertia of heat” as well as the subtle interplay between reversible and irreversible couplings. We provide two concrete applications of our procedure, where the constitutive relations are evaluated, respectively, from maximum entropy hydrodynamics and Chapman–Enskog theory.

1. Introduction

The most widespread theory for relativistic dissipation is the Israel–Stewart theory [1,2], which has proven to be very effective in modeling viscosity and heat conduction in relativistic gases [3] and liquids [4]. The rationale of the Israel–Stewart theory is rooted in Extended Irreversible Thermodynamics [5], which posits that the dissipative fluxes, like the heat flux q, should be treated as non-equilibrium thermodynamic variables. This allows one to define a non-equilibrium entropy density $s\left(q^2\right)$, which can be expanded to the second order in the flux:

$$s\left(q^2\right)=s_{\mathrm{eq}}-{\displaystyle \frac{{\beta }_1}{2T}}{q}^2+𝒪\left(q^4\right).$$

(1)

Using the second law of thermodynamics as a guiding principle, one can then derive some dissipative equations of motion for the fluxes, which resemble Catteno’s model: $\tau \dot{q}+q\propto$ “local gradients” [6]. This model has indeed been shown to be consistent with the kinetic theory of gases [7,8] and the rheological theory of liquids [9,10,11].

It is natural to ask whether we can extend the Israel–Stewart framework beyond the second order in dissipative fluxes. (In this article, we interpret the Israel–Stewart framework as a formulation of transient hydrodynamics, and not as a gradient expansion, see [8], interpretation (iii). This means that q is a dynamical effective field, which parameterizes the displacement of the fluid from local equilibrium. In the case of heat conduction, Israel–Stewart hydrodynamics reduces to the M1 closure scheme [12,13,14]). This question was addressed in the literature for certain dissipative processes, and the resulting formalism seems to depend on the flux under consideration. For the bulk stress, the extension of Israel–Stewart beyond quadratic order can be identified with Hydro+ [15,16,17]. For the shear stress, the extension is called “anisotropic hydrodynamics” [18,19]. In this work, we focus on heat conduction, which is probably the least understood case.

The two most promising extensions of the Israel–Stewart theory for heat conduction beyond quadratic order are Carter’s multifluid theory [20,21,22] (which treats heat as a carrier of inertia) and the GENERIC theory [23]. These were recently proven to be the same mathematical system of equations [11] (just written in different variables). Both approaches define a heat-flux-dependent equation of state, similarly to (1), which can be in principle extrapolated to large values of heat flux. However, to date, no practical procedure has been proposed to compute such an equation of state from microscopic models. This article aims to propose such a procedure.

Throughout the article, we adopt the metric signature $(-,+,+,+)$ and work in natural units: $c=k_B=\hslash =1$.

2. Mathematical Structure of the GENERIC-Multifluid Theory for Heat Conduction

First, let us analyze the theory of Öttinger [23] for relativistic heat conduction. Such theory arises from direct application of the GENERIC framework [24] to relativistic conductive fluids. Since this same theory can also be derived within Carter’s multifluid framework [11], we will refer to it as the GENERIC-Multifluid (GM) theory.

2.1. Non-Equilibrium Thermodynamics

The fields of the GM theory are, by assumption, $\mathsf{\Psi }=\{n,u^{\mu },w_{\mu }\}$. The first two may be interpreted as the rest-frame baryon density and the (Eckart frame [25]) flow velocity. The covector $w_{\mu }$ is an effective non-equilibrium field, usually called “thermal momentum” [22]. The non-equilibrium temperature is defined to be $T=-u^{\mu }{w}_{\mu }$. If the fluid is in the local thermodynamic equilibrium, $w^{\mu }$ must be parallel to $u^{\mu }$ (by isotropy). It follows that the non-negative definite scalar $w^{\mu }{w}_{\mu }+T^2$ can be interpreted as a measure of how far from local equilibrium the fluid is. This motivates (in agreement with Extended Irreversible Thermodynamics [5]) the introduction of a non-equilibrium free-energy density $\mathfrak{F}(T,n,w^{\mu }{w}_{\mu })$, which has an absolute minimum at $w^{\mu }{w}_{\mu }=-T^2$, for fixed values of n and T [26,27]. We define the non-equilibrium entropy density s and chemical potential $\mu$ from the following differential:

$$d\mathfrak{F}=-sdT+\mu dn+{\displaystyle \frac{\sigma }{2}}d(w^{\mu }{w}_{\mu }+T^2).$$

(2)

At equilibrium, the thermodynamic coefficient $\sigma$ is necessarily positive (for $\mathfrak{F}$ to be in a minimum [27]). Additionally, we can define the thermodynamic energy density $\epsilon =\mathfrak{F}+Ts$ and the thermodynamic pressure $P=-\mathfrak{F}+\mu n$, as in standard thermodynamics [28]. The thermodynamic identities below follow directly from the above definitions:

$$\begin{array}{cc}\hfill d\epsilon =& Tds+\mu dn+{\displaystyle \frac{\sigma }{2}}d(w^{\mu }{w}_{\mu }+T^2),\hfill \\ \hfill dP=& sdT+nd\mu -{\displaystyle \frac{\sigma }{2}}d(w^{\mu }{w}_{\mu }+T^2),\hfill \\ \hfill \epsilon +P=& Ts+\mu n.\hfill \end{array}$$

(3)

2.2. Hydrodynamic Constitutive Relations

The effective fields $\mathsf{\Psi }$ are not observable. They are just mathematical degrees of freedom that we use to parameterize the macroscopic state of the system. Within relativistic hydrodynamics, the relevant physical observables are the following fluxes: $T^{\mu \nu }$ (the stress-energy tensor), $s^{\mu }$ (the entropy current), and $n^{\mu }$ (the baryon current). Thus, we need some formulas to express these fluxes in terms of the effective fields $\mathsf{\Psi }$. Such formulas are usually referred to as constitutive relations. For the GM theory, the constitutive relations are postulated to be

$$\begin{array}{cc}\hfill T^{\mu \nu }=& P{g}^{\mu \nu }+(\epsilon +P-\sigma T^2)u^{\mu }{u}^{\nu }+\sigma w^{\mu }{w}^{\nu },\hfill \\ \hfill s^{\mu }=& (s-\sigma T)u^{\mu }+\sigma w^{\mu },\hfill \\ \hfill n^{\mu }=& n{u}^{\mu }.\hfill \end{array}$$

(4)

These are just the most natural constitutive relations that one can write working in the Eckart frame, i.e., assuming that $n^{\mu }\propto u^{\mu }$, $n=-n^{\mu }{u}_{\mu }$, and $\epsilon =T^{\mu \nu }{u}_{\mu }{u}_{\nu }$. Indeed, we defined the heat flux vector

$$q^{\mu }=\sigma T(w^{\mu }-T{u}^{\mu }),$$

(5)

which satisfies the orthogonality condition $u_{\mu }{q}^{\mu }=0$, and we can rewrite the constitutive relations (4) as follows:

$$\begin{array}{cc}\hfill T^{\mu \nu }=& P{g}^{\mu \nu }+(\epsilon +P)u^{\mu }{u}^{\nu }+u^{\mu }{q}^{\nu }+q^{\mu }{u}^{\nu }+{\displaystyle \frac{q^{\mu }{q}^{\nu }}{T^2\sigma }},\hfill \\ \hfill s^{\mu }=& s{u}^{\mu }+{\displaystyle \frac{q^{\mu }}{T}},\hfill \\ \hfill n^{\mu }=& n{u}^{\mu }.\hfill \end{array}$$

(6)

These can be interpreted as non-perturbative generalizations of the Israel–Stewart constitutive relations [2,29,30] (in the Eckart frame). It should be kept in mind that all the thermodynamic variables may depend on the heat flux in a fully non-linear manner. Thus, the present theory is in principle applicable in regimes with large fluxes of heat.

2.3. Consistency with Relativistic Thermodynamics

Let us verify that the above theory is consistent with the principles of relativistic thermodynamics in the van Kampen–Israel formulation [1,31,32,33,34]. Using Equations (3) and (4), one can easily prove the following identities:

$$\begin{array}{cc}\hfill T{s}^{\mu }=& P{u}^{\mu }-\mu n^{\mu }-u_{\nu }{T}^{\nu \mu },\hfill \\ \hfill Td{s}^{\mu }=& -\mu d{n}^{\mu }-u_{\nu }d{T}^{\nu \mu }-2\sigma u^{[\mu }{w}^{\nu ]}d{w}_{\nu }.\hfill \end{array}$$

(7)

The first equation is Israel’s covariant Euler relation [32]. Note that, while in general this is an equilibrium identity, in the GM theory, it happens to hold also in the presence of a heat flux. The second equation coincides with Israel’s covariant Gibbs relation if and only if $w^{\mu }=T{u}^{\mu }$. This implies that the fluid is in local thermodynamic equilibrium if and only if $q^{\mu }=0$ [see Equation (5)], i.e., there is no flow of heat across the fluid. Thus, the theory is indeed consistent with the principles of van Kampen–Israel thermodynamics.

It is also straightforward to verify that the GM theory describes a multifluid of Carter [26,35]. In fact, if we express the free energy as a function $\mathfrak{F}(w_{\alpha },n^{\alpha },g^{\alpha \beta })$, we have the following partial derivatives:

$${\displaystyle \frac{\partial \mathfrak{F}}{\partial w_{\alpha }}}{|}_{n^{\beta },g^{\beta \gamma }}=s^{\alpha },{\displaystyle \frac{\partial \mathfrak{F}}{\partial n^{\alpha }}}{|}_{w_{\beta },g^{\beta \gamma }}=-\mu u_{\alpha }+{\displaystyle \frac{s-\sigma T}{n}}(w_{\alpha }-T{u}_{\alpha }),2{\displaystyle \frac{\partial \mathfrak{F}}{\partial g^{\alpha \beta }}}{|}_{n^{\beta },w_{\gamma }}=T_{\alpha \beta }-P{g}_{\alpha \beta },$$

(8)

which are consistent with Carter’s theory in the generating function formulation [36]. This implies that, as long as the second law of thermodynamics is respected, and the equation of state for $\mathfrak{F}$ is prescribed in accordance with the requirements listed in [37], the GM theory is linearly causal [38] and covariantly stable [39], both dynamically and thermodynamically [40,41].

2.4. Equations of Motion

To complete the theory, we need to prescribe some equations of motion for the fields $\mathsf{\Psi }=\{n,u^{\mu },w_{\mu }\}$. Since the algebraic degrees of freedom are eight, we need eight independent equations of motion. Out of these, five are the conservation laws ${\nabla }_{\mu }{n}^{\mu }=0$ and ${\nabla }_{\mu }{T}^{\mu \nu }=0$. The remaining three are derived to guarantee consistency with the principles of GENERIC [24]. The simplest equation of motion fulfilling all the requirements is [23]

$$u^{\mu }({\nabla }_{\mu }{w}_{\nu }-{\nabla }_{\nu }{w}_{\mu })=-{\displaystyle \frac{1}{\tau }}(w_{\nu }-T{u}_{\nu }),$$

(9)

where $\tau (T,n,w^{\mu }{w}_{\mu })>0$ can be interpreted as the relaxation time. In (9), there are only three independent equations, since contraction of both sides with $u^{\nu }$ returns a trivial identity “$0=0$”. Consistency with GENERIC automatically entails consistency with the Onsager–Casimir principle [42,43] and with the second law of thermodynamics. Indeed, with the aid of the second equation of (7), we can explicitly evaluate the entropy production rate:

$$T{\nabla }_{\mu }{s}^{\mu }={\displaystyle \frac{\sigma }{\tau }}(w^{\mu }{w}_{\mu }+T^2),$$

(10)

which is non-negative definite for arbitrary values of $\mathsf{\Psi }$.

To get a better insight into the physical content of Equations (9) and (10), we can express them in terms of the heat flux vector (5). The result is

$$\begin{array}{cc}\hfill \tau \sigma T{\mathcal{L}}_u{\left({\displaystyle \frac{q}{\sigma T}}\right)}_{\nu }+q_{\nu }=& -\tau \sigma T({g^{\mu }}_{\nu }+u^{\mu }{u}_{\nu })(T{u}^{\lambda }{\nabla }_{\lambda }{u}_{\mu }+{\nabla }_{\mu }T),\hfill \\ \hfill T{\nabla }_{\mu }{s}^{\mu }=& {\displaystyle \frac{q^{\mu }{q}_{\mu }}{\tau \sigma T^2}}.\hfill \end{array}$$

(11)

The consistency with the Israel–Stewart theory [1,2] in the limit of small heat fluxes is evident. The Lie derivative ${\mathcal{L}}_u$ in the first equation automatically accounts for the coupling with the vorticity predicted by kinetic theory [3], and it guarantees that $q^{\nu }$ remains orthogonal to $u^{\nu }$ at all times.

3. Evaluation of the Constitutive Relations from Microphysics

In the previous section, we outlined a general hydrodynamic framework for describing relativistic heat conduction non-perturbatively. Now, we need a procedure for computing the non-equilibrium equation of state $\mathfrak{F}(n,u^{\mu },w_{\mu })$ from microphysics. The main difficulty is that the field $w_{\mu }$ does not have a straightforward physical interpretation. Indeed, even T itself is not clearly defined (out of equilibrium [25]). This may open the doors to all sorts of ambiguities when trying to connect hydrodynamics with other levels of description, like kinetic theory. Here, we present a simple (and rigorous) procedure that allows one to circumvent all interpretative difficulties and to evaluate $\mathfrak{F}$ unambiguously.

3.1. General Strategy

Pick a spacetime event $\mathcal{P}$ and move to the local rest frame of the fluid. Align the $x^1$ axis with the heat flux vector $q^{\mu }\left(\mathcal{P}\right)$. Then, the constitutive relations (4) and (6) can be expressed in components as follows:

$$T^{\mu \nu }\left(\mathcal{P}\right)=\left[\begin{array}{cccc}\epsilon & \sigma Tw& 0& 0\\ \sigma Tw& P+\sigma w^2& 0& 0\\ 0& 0& P& 0\\ 0& 0& 0& P\end{array}\right]=\left[\begin{array}{cccc}\epsilon & q& 0& 0\\ q& P_L& 0& 0\\ 0& 0& P_T& 0\\ 0& 0& 0& P_T\end{array}\right],$$

(12)

where $w=\sqrt{w^{\mu }{w}_{\mu }+T^2}$ is the non-equilibrium excursion, $q=\sqrt{q^{\mu }{q}_{\mu }}$ is the heat flux magnitude, $P_L=P+q^2/\left(T^2\sigma \right)$ is the longitudinal pressure, and $P_T=P$ is the transversal pressure. Comparing the two matrices above, and recalling the second law (10), we obtain

$${\displaystyle \frac{w}{T}}={\displaystyle \frac{P_L-P_T}{q}},T^2\sigma ={\displaystyle \frac{q^2}{P_L-P_T}},T={\displaystyle \frac{P_L-P_T}{\tau {\nabla }_{\mu }{s}^{\mu }}}.$$

(13)

If one has a microscopic model for the heat flux (e.g., from kinetic theory), they can evaluate q, $P_L$, $P_T$, $\tau$, and ${\nabla }_{\mu }{s}^{\mu }$ explicitly. Note that there is no ambiguity in the kinetic definition of each of these quantities. Thus, there is no ambiguity over the exact values of w, $\sigma$, and T for a given kinetic state. Varying the state, we can reconstruct the function $\sigma =\sigma (T,n,w^2)$. Additionally, the equilibrium free-energy density ${\mathfrak{F}}_{\mathrm{eq}}(T,n)=\mathfrak{F}(T,n,w^2=0)$ is known from statistical mechanics. Therefore, from (2), one can finally compute the non-equilibrium free energy

$$\mathfrak{F}(T,n,w^{\mu }{w}_{\mu })={\mathfrak{F}}_{\mathrm{eq}}(T,n)+{\displaystyle \frac{1}{2}}{\int }_0^{w^{\mu }{w}_{\mu }+T^2}\sigma (T,n,w^2)d\left(w^2\right).$$

(14)

Once $\mathfrak{F}$ is known, all the constitutive relations can be computed through partial differentiation.

3.2. Two Simple Examples

Suppose that heat is transported by a single branch of quasi-particle excitations, which have a long mean free path and carry zero net baryon number (so their motion does not modify the value of $u^{\mu }$). For simplicity, we assume that the variable T, defined in (13), fully characterizes the energy distribution of such excitations. It follows that the non-equilibrium excursion w only affects the angular distribution of the excitation momenta, but not the magnitude of the momenta. Thus, we can express the heat flux and the pressure anisotropy in the following form:

$$\begin{array}{cc}\hfill q=& v_{\mathrm{ch}}(T,n)\mathfrak{R}(T,n){\int }_{-1}^{+1}f(cos\theta )cos\theta d(cos\theta ),\hfill \\ \hfill P_L-P_T=& \mathfrak{R}(T,n){\int }_{-1}^{+1}f(cos\theta ){\displaystyle \frac{3{(cos\theta )}^2-1}{2}}d(cos\theta ),\hfill \end{array}$$

(15)

where $\mathfrak{R}(T,n)$ is the average stress content of the excitation branch, $v_{\mathrm{ch}}(T,n)$ is the characteristic speed of the branch, and $f(cos\theta )$ is the (normalized) angular distribution of the momenta in the branch. (Note that the fluid possesses other excitation branches, with shorter mean free path (some of which carry net baryon number). All of these other excitations manage to thermalize (being short-lived). Hence, $\mathfrak{R}$ does not describe the totality of the stress trace $2P_T+P_L$. Instead, the stress trace can be decomposed as $\mathfrak{R}+3\mathfrak{C}$, where $\mathfrak{C}$ is an isotropic piece. At this stage, we do not need to model $\mathfrak{C}$ explicitly, because it does not contribute to the difference $P_L-P_T$, being isotropic). To determine f, we need a kinetic model for flux-limited diffusion. There are two popular proposals in the literature [44]. The first, due to Minerbo [45], postulates that the angular distribution f should maximize the entropy for the given value of heat flux q. The second, due to Levermore [12,46], is an approximate solution of the Boltzmann equation for the long-lived excitations, with a derivation that goes back to Chapman and Enskog [47]. Both approaches lead to a formula for f that depends on a free parameter $Z\in [0,+\infty ]$. The exact expressions for f are provided below:

$$\begin{array}{ccc}\hfill \mathrm{Minerbo}:& \hfill f\left(x\right)& ={\displaystyle \frac{Z{e}^{Zx}}{2sinhZ}},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \mathrm{Levermore}:& f(x)={\displaystyle \frac{1}{2Z(cothZ-x)}}.\hfill \end{array}$$

(17)

When $Z=0$, the distribution is isotropic. When $Z=+\infty$, all of the excitations travel in the direction of the heat flux. In Appendix A, we sketch the derivation of (15)–(17). Both models lead to the same prescription for the heat flux as a function of Z, namely $q=v_{\mathrm{ch}}\mathfrak{R}(cothZ-1/Z)$, but predict different pressure anisotropies (see Figure 1, left panel):

$$\begin{array}{cc}\hfill \mathrm{Minerbo}:& {\displaystyle \frac{P_L-P_T}{\mathfrak{R}}}={\displaystyle \frac{Z^2-3ZcothZ+3}{Z^2}},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \mathrm{Levermore}:& {\displaystyle \frac{P_L-P_T}{\mathfrak{R}}}={\displaystyle \frac{3cothZ(ZcothZ-1)}{2Z}}-{\displaystyle \frac{1}{2}}.\hfill \end{array}$$

(19)

Plugging these formulas into (13), we obtain two alternative constitutive relations $\sigma \left(w\right)$, which are plotted in Figure 1, right panel. There is no analytical expression because, in both cases, the relation $\sigma \left(w\right)$ is given in a parametric form, $\left\{w\right(Z),\sigma (Z\left)\right\}$, and the dependence of w on Z does not admit an analytic inverse. However, we can fit the relations using a polynomial approximation. Below, we report a good compromise between analytical simplicity and accuracy (see dashed lines in Figure 1):

$$\begin{array}{cc}\hfill \mathrm{Minerbo}:& {\displaystyle \frac{T^2\sigma }{\mathfrak{R}{v}_{\mathrm{ch}}^2}}={\displaystyle \frac{5}{3}}-{\displaystyle \frac{3}{2}}{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^2+{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^4-{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^6,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \mathrm{Levermore}:& {\displaystyle \frac{T^2\sigma }{\mathfrak{R}{v}_{\mathrm{ch}}^2}}={\displaystyle \frac{5}{6}}+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^4.\hfill \end{array}$$

(21)

These approximations are designed to be very accurate up to $v_{\mathrm{ch}}w/T\sim 0.5$. At larger heat fluxes, the accuracy is slightly lower (∼5 % error). However, at maximum heat flux, namely for $v_{\mathrm{ch}}w/T=1$, the polynomial approximation becomes exact. As a consistency check, we note that, in the limit of small q, the conductivity coefficient $\kappa =\tau \sigma T$, as predicted by Equation (11), has the correct scaling [48] for both models:

$$\kappa \propto {\displaystyle \frac{\tau \mathfrak{R}{v}_{\mathrm{ch}}^2}{T}}.$$

(22)

Using Equation (14), we finally obtain the non-equilibrium free-energy density for both models (see Figure 2):

$$\begin{array}{cc}\hfill \mathrm{Minerbo}:& \mathfrak{F}={\mathfrak{F}}_{\mathrm{eq}}+\left[{\displaystyle \frac{5}{6}}{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^2-{\displaystyle \frac{3}{8}}{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^4+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^6-{\displaystyle \frac{1}{48}}{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^8\right]\mathfrak{R},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \mathrm{Levermore}:& \mathfrak{F}={\mathfrak{F}}_{\mathrm{eq}}+\left[{\displaystyle \frac{5}{12}}{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^2+{\displaystyle \frac{1}{36}}{\left({\displaystyle \frac{v_{\mathrm{ch}}w}{T}}\right)}^6\right]\mathfrak{R}.\hfill \end{array}$$

(24)

These are the equations of state we were looking for.

4. Conclusions

At present, we don’t know whether the GM theory is applicable outside of the Israel–Stewart regime. For that to happen, three conditions need to be met. First, the dynamics of heat must be fully characterized by a single non-equilibrium structural variable $w_{\mu }$. Secondly, it should still be possible to define an extended thermodynamic theory involving the non-equilibrium excursion parameter $w^{\mu }{w}_{\mu }+T^2$. Finally, the equation of motion for $w_{\mu }$ should be governed by GENERIC dynamics. And all of this must be true for large values of heat flux. Admittedly, these are quite strong assumptions to digest. However, given the success of GENERIC in describing complex fluids [49], it may happen that certain relativistic liquids indeed fulfill the requirements.

The main danger when dealing with “far-from-equilibrium” theories of this kind is the risk of non-falsifiability. There is so much freedom in the construction of the non-equilibrium equation of state $\mathfrak{F}(T,n,w^{\mu }{w}_{\mu })$ that it is virtually possible to fit any given data a posteriori, by simply adjusting the equation of state to the needs. This would likely result in overfitted fluid models. To avoid this problem, one should know the non-equilibrium equation of state before fitting the data with hydrodynamics. Ideally, the (theoretical) error bars of the non-equilibrium equation of state should be smaller than the (experimental) error bars of the data points.

Here, we have proposed a simple procedure for evaluating $\mathfrak{F}(T,n,w^{\mu }{w}_{\mu })$ for any given microscopic model. This procedure has the advantage of being free of intrinsic uncertainties. Rather than coming up with a statistical interpretation of $\mathfrak{F}$, which would suffer from ambiguities related to the unclear microscopic definition of $w_{\mu }$, we adopted a more rigorous approach: We showed that T, $w^{\mu }{w}_{\mu }$, and $\partial \mathfrak{F}/\partial \left(w^{\mu }{w}_{\mu }\right)$ can all be expressed in terms of physical observables that are unambiguously defined in any kinetic theory. Therefore, if the GM theory holds for large values of q (which is admittedly a big “if”), there is one and only one free energy for each given microscopic model. This makes the GM theory at least falsifiable.

We have tested the method with two simple kinetic models of flux-limited diffusion. The first, due to Minerbo [45], based on the maximum entropy principle, and the second, due to Levermore [12,46], based on the Chapman–Enskog procedure [47]. The resulting non-equilibrium free energies are reported in Figure 2. Their qualitative behavior is reasonable. For example, the non-equilibrium deviations of the free energy are of the order of the pressure anisotropy, which is rather natural, considering that $\mathfrak{F}=\mu n-P$. Indeed, all the scaling laws agree with the expectations.