Open mathematical aspects of continuum thermodynamics: hyperbolicity, boundaries and nonlinearities

Thermodynamics is continuously spreading in the engineering practice, which is especially true for non-equilibrium models in continuum problems. Although there are concepts and approaches beyond the classical knowledge, which are known for decades, their mathematical properties and consequences of the generalizations are less-known and are still of high interest in current researches. Therefore, we found it essential to collect the most important and still open mathematical questions related to different continuum thermodynamic approaches. First, we start with the example of Classical Irreversible Thermodynamics (CIT) in order to provide the basis for the more general and complex frameworks, such as the Non-Equilibrium Thermodynamics with Internal Variables (NET-IV) and Rational Extended Thermodynamics (RET). Here, we aim to present that each approach has its specific problems, such as how the initial and boundary conditions can be formulated, how the coefficients in the partial differential equations are connected to each other, and how it affects the appearance of nonlinearities. We present these properties and comparing the approach of NET-IV and RET to each other from these points of view.


Introduction
The famous finding of Fourier, widely known as Fourier's law in which the heat flux q i is proportional with the temperature gradient ∂ i T , serves the basis for almost all the thermal problems in the engineering practice, where λ is known to be the thermal conductivity. Here, and throughout the manuscript, we use index notation with Einstein's summation convention. However, as many experimental results prove, eq. (1) must be generalized in order to provide a reliable model both for describing heat waves and room temperature phenomena [1,2,3,4], occurring on a wide range of spatial scales. This is a task that must be carried out carefully, and in which thermodynamics has an inevitable role with increasing importance. Independently of how one generalizes eq. (1), its consequences must be understood before starting to choose and utilize a model in practical problems, beyond simply trying to model an experimental setting. This stands as a motivation for the present paper: discussing how the thermodynamic modeling process restricts the mathematical properties of the resulting partial differential equations, and how it aids (or fades) our understanding of a physical phenomenon. For instance, generalizing eq. (1) modifies the possibilities for initial and boundary conditions, their analytical and numerical solutions, and their compatibility for coupled problems, i.e., what is a physically feasible, reasonable and solvable model. Therefore, one needs a framework that is constructive, easy to handle for practical problems but still provides insight that is enough for the problem.
In order to fulfill these requirements, many approaches are developed in the last decades. The framework of Classical Irreversible Thermodynamics (CIT) is well-understood for classical phenomena such as diffusive heat conduction or coupled problems (diffusive heat and moisture transport, Peltier and Seebeck effects, etc.) [6,5]. Even this classical approach of CIT contains some interesting and less-known possibilities for modeling. Section 2 is devoted to present the usual procedure of CIT, with discussing its mathematical consequences. Section 3 aims to continue and broaden the classical approach by introducing the so-called internal variables. This is the approach of Non-Equilibrium Thermodynamics with internal variables (NET-IV) [13,18,4]. Here, we focus mainly on nonlinearities, initial and boundary conditions, and its connection with the framework of Rational Extended Thermodynamics (RET) [19,20], b y discussing briefly an intermediate approach, called Extended Irreversible Thermodynamics (EIT) [21]. Both RET and EIT follow the results and procedure of kinetic theory, but on a different level, this is the reason why EIT is called to be intermediate between NET-IV and RET [22]. Therefore, Section 4 presents the basic structure of RET together with its mathematical consequences. More closely, we concentrated on the closure problem particulary, which affects the boundary conditions differently than in the case of NET-IV. In the sense of methodology and problems which are analyzed, this study complements recent review [23] and can be regarded as the other side of the same coin. Table 1. Summarizing the balance equation

Integral form
From now on, we use the upper dot notation for the substantial time derivative, i.e.,˙≡ d dt = ∂ t + v i ∂ i , which shows how a quantity changes in time at a fixed material point while it moves with velocity v i , and, ∂ t and ∂ i are the partial time and space derivatives, respectively.
With omitting the details, the following substantial balance equations can be derived for the thermomechanical 1 processes of non-polar 2 continua: (1) Conservation of mass implies the so-called continuity equatioṅ (2) Cauchy's first equation of motion ensures the conservation of linear momentum. Here σ ij denotes the Cauchy stress tensor, which is the opposite of the pressure tensor P ij and g i is the field strength. (3) The conservation of angular momentum implies the symmetry of the Cauchy stress, i.e., which is valid only for non-polar media, and usually called as Cauchy's second equation of motion. (4) Starting with the conservation of total energy, one can derive the balance of internal energy or the first law of thermodynamics, which reads as where e denotes the specific internal energy, q i the heat flux density (more precisely the conductive current of internal energy), q V is the volumetric heat source density, σ ij ∂ j v i is the mechanical power, for which σ ij ∂ j v i = σ ij (∂ j v i + ∂ i v j ) /2 holds in the light of (4).
These balance equations -prescribed for conservative quantities -are governed by the entropy balance, which is in substantial form, and s is the mass specific entropy, (J S ) i is the entropy current density and Σ S is the entropy production, which is non-negative and measures the dissipation along a process. While the first law of thermodynamics is given in (5), the second law is more complex. Although, it is expressed formally by (6), we provide here a more detailed definition, following Ván [14].
(1) There exists independent themodynamical bodies, which can be described by extensive state quantities and intensive state functions. One body is characterized by N extensive state quantities (X a , a ∈ 1, 2, . . . , N ). The thermodynamical state space is spanned by the Descartes product of all X a . We will denote the intensive state functions with Y a , a ∈ 1, 2, . . . , N , which are the functions of the extensive ones. For fluids and gases with one component, the extensive state quantities are the mass M , the volume V and the internal energy E. (2) The entropy is the potential function of the vector space spanned by the intensive thermodynamical state quantites, characterized by the Gibbs relation as In other words, the intensive state functions are the partial derivatives of entropy, i.e., For fluids and gases the Gibbs relation is where the intensive state functions are defined by the partial derivatives in which T denotes the temperature, p is the (static) pressure and µ is the chemical potential. (3) The entropy is an extensive thermodynamical function, called as locality condition. More precisely, (a) the entropy is an Euler homogeneous function of any of its variables, e.g., (b) for all scalar X A , one can be introduce the X A -specific entropy, which is the function of the corresponding specific quantities, e.g., the X 1 -specific entropy s is defined as Usually the mass specific entropy s(e, v) and the entropy density ρ S (ρ E , ρ) is used, where e is the (mass) specific internal energy, v is the specific volume, and ρ E = E V is the internal energy density. (c) The Euler relation for entropy reads as which, for fluids and gases, is In the light of (13), the entropy can be given as where the expression in the bracket is the mass specific entropy. As a consequence of statement 3, the Gibbs relation follows for s(e, v): 3 (4) The X 1 -specific entropy is a concave function of its variables, i.e., This statement leads to the internal or material stability criteria. Particularly for fluids and gases, these criteria are where c v is the isochoric specific heat and χ T is the isothermal compressibility. All these together form the concept of entropy, and the second law is expressed through its balance equation. Moreover, CIT has one more crucial assumption, which is called local equilibrium hypothesis. Usually, the classical literature defines it as the state of a material point is close to the equilibrium, that is, all the spatial points of a continuum is close to equilibrium, and quasi-equilibrium states are following each other from time instant to time instant. Such definition suffers from many ambiguities, e.g., how the notions 'close', and 'quasi-equilibrium' can be defined. Instead, we recommend using the following definition. Local equilibrium is when the Gibbs-relation remains valid out of homogeneous equilibrium. In other words, the same X a extensives are used to characterize both the equilibrium state and the process out of equilibrium.
However, it is clear by now that the local equilibrium hypothesis has a limited region of validity. For instance, heterogeneous materials, fast processes, very high heat fluxes, low or high-temperature conditions, etc., are all capable of resulting in a non-equilibrium phenomenon. Consequently, in the generalizations of CIT, other variables are introduced as well, called non-equilibrium variables.
2.1. Exploiting the second law. Equations (2), (3) and (5) determine the time evolution of mass, momentum and energy fields, and containing five unknown quantities (ρ, v i , σ ij , e, q i , but g i and q V are prescribed functions) overall. From mathematical poing of view, one needs a 'closure', the relation between densities and their fluxes. These equations are the so-called constitutive equations, (roughly, the material laws), which describe the relationships among (J X ) i and x.
Thermodynamically consistent constitutive equations can be derived with the help of entropy and entropy production. The reversible-or the entropy preserving-contribution of the constituve equations is embedded into the entropy function, while the irreversible-or dissipative-contribution, which increases entropy (thus results in thermodynamically stable processes) can be derived via evaluating the entropy production. For mathematically rigorous techniques, we refer to the work of Colemann-Noll and Liu [15,16].
In what follows, we provide a simple example by utilizing a mathematically less precise, but and 'easier to understand and use' method. Example: Heat conduction in rigid bodies -Fourier heat conduction. Here, we consider a rigid body at rest w.r.t. a given reference frame, thus both the velocity field and the volume change are identically zero. Furthermore, we neglect the volumetric heat source density, too. Consequently, the balance of internal energy and the Gibbs relation are Let us evaluate now the l.h.s. of the entropy balance (6), constrained by the balance of internal energy (19). With separating a full divergence term, we have First, by separating the full divergence from ρṡ, we obtained the usual definition of the entropy current density, i.e., it is not necessarily an apriori assumption on contrary to the usual approach, and reads In case of the generalized theories, (21) is different and offers the possibility to include spatially nonlocal effects in the model (roughly speaking, higher-order spatial derivatives). Finally, one must ensure the positive semi-definiteness of the entropy production for which the simplest solution is to assume a linear dependence on reciprocal temperature gradient, where the heat current density appears to be as a linear function of ∂ i 1 T . Moreover, l (T ) is a positive function of temperature, and λ(T ) is called thermal conductivity. There are several -infinite -other opportunities to ensure the positive-semidefiniteness of entropy production. For example, where l 1 (T ), l 3 (T ), l 5 (T ) are positive functions, and the positiveness of Σ S is immediately ensured. Analogously, assuming l(T, ∂ i T ) coefficient results in a similar solution. It is still unknown how these nonlinear solutions could relate to the experiments. Overall, heat conduction in a rigid body can be described by the balance of internal energy, and the constitutive equation, called Fourier's law in this case: As we mentioned above, the temperature is a function of the internal energy. Assuming that T (e) is invertable, and the specific internal energy is linear in the temperature with the specific heat c as e(T ) = cT , then substituting (26) into (25), the heat conduction equation for temperature is which is a parabolic second order partial differential equation. From mathematical point of view, the usual boundary conditions are the following.
(1) First kind or Dirichlet boundary condition: the temperature is prescribed on the boundary: here t and r i are the time and space coordinates, respectively. (2) Second kind or Neumann boundary condition: the normal component of the heat flux is prescribed on the boundary, that is, Applying the Fourier equation (26), (29) can be rewritten to (3) Third kind or Robin boundary condition (convection type): the linear combination of the temperature and the normal component of heat flux is prescribed on the boundary. For instance, the Newtonian cooling boundary is where h is the heat transfer coefficient and T ∞ is the ambient temperature. Here we add some comments on the boundary conditions. Let us assume, that the material parameters c and λ can be considered as constant. In this case one can introduce the thermal diffusivity α := λ ρc , and the heat conduction equation reduces toṪ However, since the differential operator˙− α∂ i,i is linear, a similar equation can be derived for the heat current density as well,q in which q i stands as the only field variable. These can be called as the temperature and heat flux representations of the system (25)- (26). For some particular situations, such as the modeling of the laser flash experiment, qrepresentation fits well. Notably, in this case, the usual second kind boundary becomes first kind, and it is not possible anymore to prescribe the temperature directly on the boundary. However, knowing the q i field, the temperature can be easily recovered by integration, with a further restriction on a reference point. Similarly, the constitutive equations for heat conducting Newtonian fluids, i.e., the Navier-Stokes-Fourier equations is also possible to derive in the same way [4].

Non-Equilibrium Thermodynamics with Internal Variables
While the classical thermal and mechanical equations -Hooke, Navier-Stokes, Fourier -are enough for the majority of engineering problems, these models have their limitations as any other model as well. For instance, modeling rheological phenomena occurring in solids require the generalization of Hooke's law. One possible way to do so is to introduce dynamic degrees of freedom [17], aimed to catch the strain rate dependence by including a non-equilibrium state variable. A more general approach is to use internal variables, which are not necessarily zero in equilibrium, contrary to dynamic degrees of freedom. This is more apparent in situations related to microcrack modeling, in which the density function of cracks can be non-zero in equilibrium [24]. In the following, we place our focus on internal variables through the example of the Maxwell-Cattaneo-Vernotte (MCV) and Guyer-Krumhansl (GK) equations [25,26,27,28,29].

Maxwell-Cattaneo-Vernotte equation.
In both cases, we follow the usual procedure of non-equilibrium thermodynamics [18,4], in which the classical, local equilibrium part of the specific entropy s is shifted by a positive semidefinit function of the internal variable ξ i in a way to preserve the convexity properties of the entropy: As a special choice, one can decide to identify ξ i with q i , and thus achieve compatbility with EIT. Applying the classical entropy flux (J S ) i = q i /T , one obtains the following expression for the entropy production with considering ξ i ≡ q i : and the balance of internal energy, is exploited, in which the source terms are neglected. Here we encounter multiple possibilities how to solve the entropy inequality (35), similarly to the Fourier equation. First, seeking linear solutions for isotropic materials, we find the following in which one can define the thermal conductivity λ = l/T 2 and the relaxation time τ = ρml, and obtain the final form of the MCV equation as Although this derivation is well-known since Gyarmati [30], there are multiple mathematical aspects, which must be discussed in detail. First, the Eq. (37) shows a sort of 'basic form' of the constitutive equation, where the relationship among the coefficients is visible. It is of high importance when the simplest nonlinearities, e.g., temperaturedependent thermal conductivity is tried to be implemented: it immediately affects the relaxation time as well with far-reaching consequences [31]. Therefore, it is not possible to arbitrarily include such nonlinearities without knowing the 'basic form' of a constitutive equation. It is more apparent for the Guyer-Krumhansl equation. Furthermore, since q i became a state variable, it is thermodynamically possible to assume l = l(T, q i ), which case is entirely open for future investigation both theoretically and experimentally. An analogous but distinct approach is also published recently [32], in which the nonlinearities are treated on a different ground with a similar q i dependence. Second, Eq. (37) is only the simplest linear solution, and the entropy inequality enables to find infinitely many nonlinear ones, too. For instance, the 'simplest' nonlinear solution arises as a polynomial of ∂ i 1 T − ρmq i . However, even the linear case is able to bring up difficulties, still open questions both on physical and mathematical grounds. Such a difficulty emerges when non-uniform initial temperature distribution is accounted. In case of Fourier's law, it is quite straightforward how to handle that situation: the gradient of the initial temperature field immediately leads to the initial state of the heat flux field. In case of the MCV equation, one has to solve the differential equation (38) first, together with a constraint to determine the integration coefficient. It seems to be 'trivial', but let us keep in mind that the MCV equation is the simplest extension of the Fourier equation with merely one time derivative. 6 3.2. Guyer-Krumhansl equation. The situation becomes more difficult for the GK equation. Its derivation requires a generalized entropy flux as well, beyond the specific entropy (34). In the framework of EIT, (J S ) i = q i /T + µ∂ i q j q j , (µ ∈ R + 0 ) is assumed, which can exactly be recovered in the approach of NET-IV using a current multiplier, or Nyíri-multiplier B ij [33]. Here we have two options: either (J S ) j = B ij q i or (J S ) j = q j /T + B ij q i are valid. Calculating the entropy production, for the first option we obtain and for the second one, it reads In both cases, B ij is defined through the entropy inequality, and can recover the EIT approach exactly by assuming the simplest linear solution. We note that both cases allow a more general solution for B ij , e.g., B ij = µ ijkl ∂ l q k in which case µ ijkl is a fourth order tensor.
In order to keep the discussion focused, let us restrict ourselves to the one-dimensional situation and find the following as a solution of the entropy inequality: which forms the GK equation: where Notably, the Nyíri-multiplier is an efficient tool to realize coupling between quantities with different tensorial order, which would not be possible otherwise for isotropic materials due to the Curie principle. Physically, the spatially nonlocal extension of a hyperbolic model, i.e., the appearance of new spatial derivatives in the constitutive equation comparing to the MCV (38), leads to a parabolic model [4]. Consequently, a Nyíri-multiplier ensures a parabolic envelope for hyperbolic equations. Also, in general the occurrence of the nonlocal behaviour implies a coupling between different tensorial order quantities. Here, it can be interpreted as a coupling between q i and ∂ j q i , which ∂ j q i represents the rate of a 'thermal pressure' [34].
Regarding the state variable dependent coefficients, the importance of the 'basic form' is more apparent here. For instance, l 2 = l 2 (T ) : R + → R + , (l 2 ∈ C 1 ) introduces the T-dependence into all other coefficients. It is shown in [31] on the example of the MCV equation that achieving linear T-dependence in τ and λ necessarily implies ρ = ρ(T ), therefore the mechanics must be included. Although λ = λ(T ) looks completely natural in the classical case, for non-Fourier problems this convenience vanishes. Moreover, if l 1 = l 1 (T ) : R + → R + , (l 1 ∈ C 1 ) holds, one encounters extra derivatives in (43), e.g., which seems to be straightforward if the GK equation is given in its 'basic form'; however, this is not the case in the literature, in general. Therefore we suggest that speaking about non-Fourier models is safe only in the case when the complete picture behind the constitutive equation is provided, especially for nonlinear problems. Otherwise, the implementation of any state variable dependence is unwise and implies unkown simplifications as well.
The treatment of non-uniform initial conditions becomes significantly difficult contrary to Fourier's law. Since the constitutive equation (43) consists of both time and spatial dependence, one must include the boundary conditions in order to determine the initial heat flux distribution, too. This is unusual in the classical approach. Furthermore, even the initially uniform (hence ∂ i T = 0 i ) temperature distribution could allow a more general, possibly nonuniform heat flux distribution. The only requirement is to satisfy the following partial differential equation at the initial time instant: The trivial solution is a constant q at t = 0, enabling the non-zero case, too. Now, it reveals the next obstacle to applying the classical approach for boundary conditions. For the Fourier equation, q i and T are connected by the equality (1), thus both field variables are adequate to define the usual boundary conditions. This does not hold for the GK equation, the nonlocal term ∂ xx q does not allow the parallel prescription of boundaries for both variables [35,36]. This is the primary disadvantage of finite element methods in which all field variables are placed onto the same node, hence producing false solutions [37]. This attribute of the GK equation can be handled in two ways. The first one is more suitable for analytical solutions; choose a primary field variable, which is used to prescribe the boundary conditions. Then, starting from the 'basic form', it is possible to eliminate the other field variable, at least in the linear regime. This can be either T or q-representation of the GK equation. It could also be advantageous for the Fourier heat equation. The second option suggests keeping the evolution equations 'untouched', i.e., in the form of a system of partial differential equations. When the evolution equations are thermodynamically compatible, it is possible to use a numerical method with staggered fields [35], in which the primary field variable is kept on the boundary, while the other is shifted. Consequently, the shifted variable does not require to define boundary conditions, avoiding the problematic parallel description, and the structure of the evolution equations ensures the applicability of staggered fields. As a consequence of the proper treatment of boundary conditions, the solutions will satisfy the maximum principle, therefore the unphysical ones are excluded. This is not the case in [38], where the temperature evolution presents negative values as well.

3.3.
Further remarks about boundary conditions. Finally, we mention a few general remarks about the formulation of boundary conditions. Let us recall the Robin-type boundary, In fact, the Fourier's law is its differential analogue, and that type of boundary represents the interaction between a point-like body and its environment through a thermal resistance. However, in principal, exchanging the constitutive equation would affect that definition, too. For example, considering the MCV equation, analogously the convective type boundary would be One could argue for and against, but our aim here is to reveal that it is still an open, unsolved problem for non-Fourier models. This idea is already revealed in [39] for heat conduction models based on kinetic theory. However, the experimental verification is still missing, and based on [39], the differences could be negligible in certain cases, but overall, there are no decisive results. For a more detailed discussion about the boundary condition possibilities in internal variable theories, we refer to the work of Cimmelli [40].

3.4.
Coupled heat and mass transport. Previously, we focused our discussion on the structure of heat conduction equations. The same obstacles emerge for coupled systems, too; however, that parts serves as a bridge between NET-IV and RET. As the simplest demonstration for generalized coupled equations, we take the example of the extended Navier-Stokes-Fourier (NSF) equations. Their detailed thermodynamic background can be found in [22], presenting their basis for approaches of NET-IV and RET. Briefly, the classical, local equilibrium part of the specific entropy s(e, ρ) (with ρ being the mass density) is extended in the same way as we have seen on the example of heat conduction, but now the mass density also plays an essential role. Thus the specific entropy for the generalized case is in which Π ij is the viscous pressure, denotes its deviatoric part. Since the spherical and deviatoric parts are linearly independent of each other, that separation must be realized from the beginning, and utilized for the Nyíri-multiplier in the entropy current as follows: Omitting the detailed calculation of entropy production σ s , we arrive to 8 with n > 0, l 11 l 22 − (l 12 + l 21 ) 2 /4 > 0, k 11 k 22 − (k 12 + k 21 ) 2 /4 > 0, (56) as the simplest linear solution for isotropic materials. Rearringing the coefficients with considering the strictly linear case, it reads where the coefficients are formed as τ 1,2,3 are the corresponding relaxation times, ν and η are the shear and bulk viscosities, respectively; the α and β are the coupling coefficients. We observe again that the appearance of the nonlocal terms necessarily implies a connection between different tensorial order quantities. While the utilization of Nyíri-multipliers allows to achieve such a coupling both for the deviatoric and spherical parts, it also results in a parabolic model (57). This parabolic structure includes further possibilities than its hyperbolic counterpart from RET (i.e., when k 22 = 0, l 22 = 0). Consequently, the usual procedure of NET-IV leads to a more general model in which the coefficients are connected by means of (56) and (58) without any further constraint, on the contrary to RET. In parallel, it could stand as a disadvantage as well due to the higher number of parameters to fit. As a side note, we may draw attention to a recent work [41] which retains the local equilibrium assumption, and thus Gibbs relation, but allows for higher order coupling for entropy current similar to the one obtained by Nyíri-multipliers. As a result, coupled constitutive relations are obtained, but not of the rate type. Staying in the linear and quasi-linear regime, the compatibility between RET and NET-IV (consequently, with EIT, too) is proved in [22]. However, the situation could be different in the nonlinear case, which becomes apparent regarding the rarefied gas experiments [42]. The importance of these experiments lies in the mass density dependence of the coefficients (58). It is revealed in these experiments that the speed of sound considerably changes with respect to the mass density, keeping the temperature and the sound frequency ω constant. According to the original measurements, the outcome is the function of the frequency -pressure ratio of ω/p only. However, when the ω and T are constants, the density remains a mere variable to change in order to reveal the change in the speed of sound. That is, decreasing the density of gas from its normal state down to the low-pressure region (e.g., 10 Pa [42]), continuousl y enhance the role of the coupling between the fluid and thermal equations together with their hyperbolic character (i.e., the relaxation times become higher and higher). Therefore, these coefficients must be mass density-dependent. Indeed, it turned out to be true and crucial.
On one hand, in order to recover the ω/p scaling in the dispersion relations, it demands a particular mapping, i.e., z(ρ) : R + → R + where z could be any coefficient from (58). Moreover, this scaling also requires constant viscosities, thermal conductivity, and ideal gas state equations. On the other hand, these requirements are contradictory with viscosity vs. mass density measurements, in which the viscosity changes considerably for both the dense and rarefied regions. While models based on kinetic theory restrict themselves to such a scaling, the freedom in NET-IV offers further possibilities; for instance, it is easily possible to implement any measured mass density function with neither reaching further obstacles nor contradicting the basic principles of NET-IV.
Overall, we emphasize again that the implementation of nonlinearities requires special attention as the behavior of the generalized constitutive equations differs from the classical ones.

Rational Extended Thermodynamics
Extended thermodynamics emerged in the sixties of the last century [43] as a method to overcome the paradox of infinite speed of pulse propagation, which appears in the Fourier model of heat conduction and propagation of shear diffusion in Newtonian fluids. Although these models are widespread in scientific literature, as well as engineering applications, their physical inconsistency with fundamental physical constraints seemed disturbing.
The solution to this problem was first found within the framework of thermodynamics of irreversible processes (TIP), and therefore called Extended TIP. Further study related it more closely, in physical aspects, to the kinetic theory of gases, especially Grad's method of moments [44], and to the hyperbolic systems of balance laws, in mathematical aspects. This direction is called Rational Extended Thermodynamics (RET), capturing at the same time methods and spirits of rational and extended thermodynamics [19].
It is clear from the outset that RET relies both on sound physical principles and desirable formal mathematical structure with the aim to establish physically/thermodynamically consistent models that predict a finite speed of pulse propagation. To that end, we shall discuss both aspects of this theory in the sequel.
The section will be organized as follows. First, basic principles and formal structure of the governing equations will be presented, along with important mathematical consequences like symmetrization. Second, an alternative approach called Molecular Extended Thermodynamics (MET), based upon the application of the Maximum Entropy Principle (MEP) will be presented. Since the major issue in non-equilibrium thermodynamics is the closure problem, we shall demonstrate how these two approaches in conjunction lead to proper and complete closure of the governing equations. Finally, we shall tackle the important problem of boundary conditions. 4.1. Formalism of RET. Rational extended thermodynamics, as its name tells, has a rational foundation with twofold meaning of the adjective: on one hand, it is based upon physical principles; on the other, it constructs appropriate physical theory by deduction. Physical principles which serve as its basis are: (1) governing equations are of balance type; they are also invariant with respect to Galilean transformations; (2) constitutive functions depend on local (in point) values of field variables; (3) governing equations are adjoined with the entropy balance law with convex entropy density and non-negative entropy production. The balance laws structure provides natural weak formulation and possibility for shock waves. Galilean invariance [45] determines the velocity dependence of constitutive quantities. The local character of the constitutive functions restricts the balance laws to the quasi-linear system of first-order partial differential equations. Finally, entropy inequality secures physical consistency, but also provides a tool for the transformation of governing equations into symmetric hyperbolic form.
If we denote by u(t, x) ∈ R n the vector of field variables, then generic form of governing equations in RET is: where F(u) ∈ R n is vector of densities, F k (u) ∈ R n is vector of fluxes, Φ k (u) ∈ R n is vector of non-convective fluxes, and f (u) ∈ R n is vector of productions. Every solution u(t, x) of the system (59) is called thermodynamic process.
If the vector of field variables is split into velocity v ∈ R 3 and vector of objective quantities w ∈ R n−3 , u = (v, w), then invariance with respect to Galilean transformations restricts the velocity dependence of densities, fluxes and productions by means of matrix function X(v): where hat denotes velocity independent functions determined by: where X(0) = I. Entropy inequality is a supplementary balance law that reads: where h is entropy density, h k is entropy flux, and Σ is entropy production. Since (61), like balance laws (59), must be invariant with respect to Galilean transformations, we have: It restricts the structure of balance laws in more fundamental way: it must be satisfied for any thermodynamic process, i.e. for any solution of (59). Since (59) and (61) are both quasi-linear partial differential equations, they must be compatible. There comes entropy principle into play. In procedure established by Liu [16], entropy balance law (61) is regarded as fundamental equation, and balance laws (59) are treated as constraints. Following the method of multipliers, problem with constraints is transformed into a problem without them by introduction of the co-vector u ′ (u), so that following compatibility condition is established: It is easily transformed into: dh = u ′ · dF, dh k = u ′ · dF k , Σ = u ′ · f ≥ 0. (62) We should make a caesura at this stage, before we proceed with further analysis of the entropy principle. Namely, compatibility conditions (62) bring two important novelties. First, once the co-vector of multipliers u ′ is determined from (62) 1 , (62) 2 formally determines the structure of entropy flux. In other words, entropy flux is no longer a priori defined, but becomes a constitutive quantity [46]. Second, residual inequality (62) 3 determines the structure of source terms so that entropy production is semi-definite. This clearly implies that dissipation comes from production terms.
If we rearrange the compatibility condition (62) 1 , the following result emerges: where h ′ is a new scalar potential, conjugate to h. Applying the same procedure to (62) 2 , we obtain new potentials by means of Legendre transform of entropy density and entropy flux with respect to densities and fluxes, respectively: Convexity of the entropy density h with respect to F allows the transformation F → u ′ to be invertible, not only locally, but globally. Then, from (63) one obtains: (59) can be put into symmetric form: Because of the privileged character of multipliers u ′ , that now act as new field variables, they deserved special name the main field [47], or entropy variables [48]. A few remarks may be given at this stage. Liu's procedure [16] brought a new concept in the study of compatibility between governing equations and entropy balance law. It introduced the multipliers as mediators, but their role went beyond compatibility issues. It appeared that symmetric form of balance laws coincided with Godunov's [49] analysis of gas dynamics equations, and symmetrization results of Friedrichs and Lax [50]. This structure was further exploited in [51], and the hierarchy of hyperbolic principal subsystems was recognized. Furthermore, it was shown in the same study that sub-characteristic condition for the characteristic speeds appears as a consequence of sub-entropy law. Finally, hyperbolic systems of balance laws with convex entropy were subject of mathematical analysis of smoothness and global existence of solutions [48,52].
Although the formalism of RET looks rather abstract at first sight, it was applied to numerous particular systems describing different physical processes. From the very beginning, its principal application was related to non-equilibrium thermodynamics, and much of its popularity was gained by remarkable equivalence with Grad's moment method in the kinetic theory of gases [44], especially 13 moments model. It was recognized that in balance laws (59) densities and fluxes may be endowed with hierarchical structure, with increasing tensorial order of components: In kinetic theory of gases F i1i2...ip are moments of the distribution function of order p. In principle, this hierarchy could be infinite, and in that case it is equivalent to Boltzmann equation. However, when it is truncated at certain order, one must close the system because productions and the last flux are undetermined. Grad expanded the velocity distribution function in terms of generalized Hermite polynomials, and resolved the closure problem by construction of its approximate form. RET, instead, relies on entropy principle and exploits it for the closure problem.
Relation to the kinetic theory of gases gave RET additional boost, but also put the limits: for a long time, it was restricted to the analysis of monatomic gases, since there was no systematic procedure for derivation of a hyperbolic system of balance laws for polyatomic gases. This obstacle was overcome recently. Starting from the first result [53], a vast amount of results were accumulated and is still growing. An up-to-date collection of main results may be found in [20]. The crucial distinction with respect to classical monatomic theory is that governing equations are in two hierarchies instead of one. In particular, the model of viscous, heat conducting, polyatomic gas is described by 14 moments, whose governing equations are structured as follows:

11
F −hierarchy is momentum-like, while G−hierarchy is energy-like, and in physical variables their densities read: where t ij = −(p + Π)δ ij + σ ij is the stress tensor, p is the hydrostatic pressure, Π is the dynamic pressure, and σ ij is the traceless viscous stress, while q i is the heat flux. In such a way, the last equation of F −hierarchy in (65) becomes balance law for stress tensor, while the last equation of G−hierarchy becomes balance law for heat flux. It must be emphasized that, except in the simplest cases of 5 or 6 moments equations, closure obtained by means of Liu's method of multipliers can be achieved only in the neighborhood of local equilibrium. In that sense it is not exact, but approximate. As a consequence, productions P ij and Q i , and entropy density are determined only approximately. In particular, we have the productions: and the entropy density and entropy flux: where s is specific entropy in equilibrium. In (67), τ S , τ Π and τ q are relaxation times, that are functions of objective quantities, typically mass density ρ and temperature T .
Here we face important limitations in the practical application of described closure procedure. First, higher-order theories admit only approximate closure, thus restricting the application to processes in a certain neighborhood of the local equilibrium. This is also manifested through the loss of hyperbolicity in certain regions of state space. That was the subject of [54]. Second, production terms are determined up to scalar positive functions, relaxation times, that are related to transport coefficients. This is a typical limitation of macroscopic theory, and cannot be overcome within its framework. To do that, one must step outside, either towards some more refined theoretical approach, or towards experimental data.

4.2.
Molecular extended thermodynamics-application of the maximum entropy principle. While RET successfully combined the methodology of continuum mechanics, i.e., application of the entropy principle, with the hierarchical structure of balance laws resembling the moment method of Grad, molecular extended thermodynamics (also termed extended thermodynamics of moments), relied on the application of Maximum Entropy Principle (MEP). MEP is a generalization of the well-known result of classical thermodynamics that entropy reaches a maximum in equilibrium.
In the kinetic theory of gases, state of the system is described by velocity distribution function f (t, x, ξ), where ξ ∈ R 3 is the molecular velocity. Its evolution is determined by the Boltzmann equation: where Q(f, f ) is integral collision operator that describes the mechanism of interaction between the particles. In this context MEP acquires variational character. Namely, it assumes that the actual state of the system is monitored through macroscopic quantities determined as moments of f . Therefore, for the kinetic entropy defined as functional h = h[f ], the actual distribution function is the one that maximizes h, subject to constraints which define macroscopic variables as moments of f .
The solution of the variational problem depends on the choice of macroscopic variables. If the constraints are to be ρ, ρv, and ρε, one obtains the Maxwellian equilibrium distribution function as a maximizer. Kogan [55] showed that extending the list of constraints by the stress tensor and heat flux leads to Grad's distribution function as a maximizer. This gave a stimulus to further application of MEP and led to a striking conclusion that exact solution of the variational problem formulated above is equivalent to the results of extended thermodynamics obtained by the procedure described in the previous subsection. Here, an explanation should be given for the word exact. Namely, the solution to the variational problem is not particularly challenging-it is an exponential function of the polynomial of molecular velocities multiplied by appropriate Lagrange multipliers. However, computation of Lagrange multipliers in terms of macroscopic variables could be problematic if the largest degree of the moment is odd-in that case, the moments are diverging. That was observed already by Kogan, and it was further analyzed by Dreyer [56], leading to the approximation of the exact solution in the neighborhood of Maxwellian distribution. That was the way in which the closure problem was resolved in the case of odd highest moments. More mathematically intoned formulation of the MEP was given by Levermore [57].
Although MEP may be regarded as just another route to the hyperbolic system of balance laws already given by RET, there is an interest in its development. It may fill the gap between the endpoint reached by RET and the fully closed system of balance laws. An illustrative example will be the hyperbolic system of balance laws for polyatomic gases.
Starting from the statistical mechanical considerations [58], it was shown [59,60] that one may use a single additional continuous variable I to describe the internal degrees of freedom of molecules in polyatomic gas, and to recover the internal energy density of polytropic gases. Within this kinetic framework, it was formulated a MEP for polyatomic gases [61] that recovered the results (65)-(68) of RET. More precisely, for the kinetic entropy defined as: one searches for the maximizer f (t, x, ξ, I) subject to constraints: where C = ξ − v is the peculiar velocity. A non-negative measure ϕ(I)dI is property of the model. By taking ϕ(I) = I α , for α > −1, one obtains the following 14 moments velocity distribution function: where f E and q(T ): are equilibrium distribution and auxiliary function, respectively. Densities (66) are recovered as moments of f 14 through (70), as well as entropy and entropy flux (68) using the relation α = (D − 5)/3. However, the main advantage comes from computation of the productions (67) as moments, i.e. weak form of the collision operator: In [62] it was used collision cross section that describes variables hard spheres, with an additional parameter s > −3/2, so that relaxation times were explicitly computed: This result was used to adapt the parameter s so that model yields physically correct value of Prandtl number. In such a way, closure of 14 moments model is completed by means of simultaneous application of RET and MEP.
We would like to draw attention to the question of non-linearity in the context of RET, especially related to production terms. To that end, we mention two examples. The first one is the multi-temperature mixture of Euler fluids. The model was established within the framework of RET in [63]. Closure obtained in this case was exact. Nevertheless, there remained a gap in the sense of phenomenological positive semi-definite matrices. More precise characterization was impossible within the macroscopic approach. Only in the special case of binary mixture it was determined the explicit form of the source terms [64,65,66]. Recently [67], using the kinetic theory approach, it was determined the non-linear structure of source terms in the multi-temperature mixture. Since they were structurally different from those obtained by RET, compatibility was achieved in a linear approximation of both. In such a way, phenomenological coefficients of RET were again determined in conjunction with the kinetic theory approach. The second example is concerned with non-linear closure using MEP by means of higher degree of multipliers used for expansion in the neighborhood of Maxwellian distribution. This leads to non-linear closure for non-convective fluxes, as shown in [68,69].
4.3. Note on boundary conditions. Derivation of appropriate boundary conditions is a tremendous open problem in RET. There is no widely accepted approach, and this note will provide some of the standpoints and procedures that appeared in the literature.
One possible approach is developed and exploited in the works of Struchtrup and Torrilhon [70,71]. The main features that distinguish this approach will be briefly explained. First, analysis is built up around so-called R13 equations-parabolic regularization of Grad's 13 moments equations obtained by adding the terms of Super-Burnett order. Construction of the model may be justified as regularization in the neighborhood of pseudo-equilibrium manifold, determined by Grad's velocity distribution, just as NSF model presents regularization in the vicinity of equilibrium manifold determined by the Maxwellian [70]. To determine the unknown non-convective fluxes for higher-order equations, the authors proceed in a classical way and exploit the entropy production rate. The outcome is non-local constitutive relations that bring regularizing terms in the model.
Second, the key for the derivation of boundary conditions is the entropy production rate due to collisions of the atoms with a solid wall. It is determined from the linear form of the entropy inequality, and the entropy flux at the wall. According to [71], since the entropy production rate must be non-negative, entropy flux out of the gas wall interface must be larger than the entropy flux into the interface. This implies the bilinear structure of the entropy production rate and determines the fluxes at the boundary in terms of their driving forces. Such boundary conditions generalize the Maxwell-Smoluchowski velocity slip-temperature jump conditions. Another approach is proposed by Barbera et al. [72] in the context of heat transfer problems in rarefied gases. It was motivated by the observation that rarefied gases, and the moments which describe its state, are subject to thermal fluctuations. They are so rapid that gas cannot adjust its fields to the actual values of non-controllable data neither in the bulk region, nor at the boundary. Therefore, it is assumed that the gas adjusts its field variables to the mean values of fluctuating data. Thermal fluctuations are determined by the entropy evaluated for Grad's approximation of the velocity distribution function.
These two examples demonstrate that the derivation of boundary data for higher-order moments in RET is a demanding task, which cannot be resolved in a universal manner. In both approaches mentioned above, boundary conditions are consequences of some general physical principle/insight. However, we still need a comprehensive comparative analysis of different approaches.

Study closure
Since the seminal paper of Coleman and Noll [15], the second law of thermodynamics took the role of constraint on the admissible forms of constitutive relations. This introduced the notion of thermodynamic consistency for the continuum models. Nevertheless, certain questions still remain open, and our aim was to enlighten some mathematical aspects of continuum theories related to the models of non-equilibrium/irreversible processes.
Throughout the study, one may follow the evolution of the idea proposed by Coleman and Noll. Within the framework of CIT (Section 2), after the discussion of balance laws and the notion of entropy, the second law was exploited in the analysis of heat conduction in rigid bodies. Apart from classical result-linear relation between the heat flux and temperature gradient-a non-linear closure was proposed along with an inherent problem of determination of the phenomenological coefficients. This line of analysis was pursued in NET-IV (Section 3) in which internal variables were introduced, and generalizations of the heat conduction problem were discussed. This led to the Maxwell-Cattaneo-Vernotte, and, after the generalization of entropy flux, to the Guyer-Krumhansl equation. The study in this framework was then expanded to the problems with coupled mass and heat transport, which naturally motivated comparison with EIT and RET. Finally, within RET (Section 4), it was presented the formalism that yields the governing equations, which predict finite speeds of propagation of disturbances. The list of field variables is extended, the structure of equations is strict, and the entropy balance law yields a closure for non-convective fluxes and source terms. Although the entropy balance law retained the central role, this résumé demonstrates that the introduction of new assumptions could enrich the picture of non-equilibrium/irreversible processes. The way in which it is achieved depends on the aims and restrictions of a certain approach. For example, NET-IV relies on internal variables and (extended) entropy balance law to produce appropriate constitutive relations, which eventually get the rate type form. On the other hand, RET imposes governing equations in the form of a hyperbolic system of balance laws for an extended set of variables. To determine the phenomenological coefficients, NET-IV relies on experimental data, whereas RET may use the possibility of looking at first principles, i.e., kinetic theory of gases, to close the system.
Another issue that may be followed is the problem of initial and boundary conditions. This is a problem that has two equally important sides-mathematical and physical. Boundary conditions should be formulated such that the well-posedness of the problem is preserved. At the same time, the data for initial and boundary conditions are supposed to be measured or controlled, and this may cause the problem in generalized approaches. Our review shows that within CIT, it is a classical and well-understood problem, both mathematically and experimentally. NET-IV proposes a generalization of convection-type boundary data, based upon constitutive relation, in a manner that is analogous to CIT. However, RET, which is primarily concerned with the behavior of rarefied gases, grabs for physical arguments that may yield appropriate information about the conditions on the solid boundary. From this perspective, it becomes apparent that there is no unifying approach to the problem of boundary conditions. It is still an open and challenging field for study.
Lastly, there is one more question which will remain unanswered, perhaps because it may be related to taste or philosophical standpoint. As mentioned in analysis of NET-IV, equations obtained in this way may be regarded as parabolic envelope to corresponding (subordinate) hyperbolic system. On the other hand, it is a standard procedure in RET to recover the parabolic (thermodynamic) limit by asymptotic methods in the small relaxation time limit [73] (even R13 equations may be obtained in this way starting from higher order hyperbolic system). It is really a matter of one's preferences which system will be regarded as "fundamental", and which one will be treated as its approximation.
As an outlook, we also mention the framework of GENERIC (General Equation for the Non-Equilibrium Reversible-Irreversible Coupling) [74,75,76,77], which can be interpreted as a thermodynamically consistent generalization of Hamiltonian mechanics, and as such, it raises additional mathematical questions, such as the geometric structure behind thermodynamics. GENERIC derives the reversible and irreversible contributions separately from each other, which looks artificial in some cases, such as Fourier heat conduction as a purely dissipative process, but in a more general and complex situation, it is a productive method for deriving dynamical models. That sort of separation can be useful in the construction of specific numerical methods, see, e.g., [78,79,80,81]. The comparison of GENERIC with other non-equilibrium thermodynamical frameworks is currently under investigation [82,83].