Generalized ballistic-conductive heat transport laws in three-dimensional isotropic materials

General constitutive equations of heat transport with second sound and ballistic propagation in isotropic materials are given using non-equilibrium thermodynamics with internal variables. The consequences of Onsager reciprocity relations between thermodynamic fluxes and forces and positive definiteness of the entropy production are considered. The relation to theories of Extended Thermodynamics is discussed in detail. We provide an explicit expression for all the components of the matrices of the transport coefficients. The expressions are cumbersome but are expected to be useful for computer programming for simulations of the corresponding physical effects.


Basic equations of heat transport coupled with a tensorial internal variable
We consider the balance equations of a rigid heat conductor, i.e. the balance of internal energy and the balance of entropy ρė + q i,i = 0, (1) ρṡ + J i,i = σ (s) . ( Here ρ is the density, e the specific internal energy, q i the current density of the internal energy, the heat flux, s the specific entropy, and J i denotes the entropy flux. The σ s entropy production rate plays a central and constructive role in the theory. i, j, k are spatial indices related to Descartes coordinates, but they can also be considered as abstract spatial indices of vectors and tensors in the sense that they do not refer to particular coordinates [46]; however, it is convenient in case of higher than second-order tensors. A comma in lower indices is for spatial derivation, and upper dot denotes the substantial time derivative (e.g.ė = ∂ t e + v i e ,i , where ∂ t is the partial time derivative). In case of rigid conductors at rest, the relative velocity of the continuum is zero; therefore, the substantial time derivative is equal to the partial time derivative. Regarding the general usage of abstract indices in classical nonrelativistic continuum theories see, e.g. in [47,48]. We introduce an additional internal variable Q i j (a second-order tensor) which will incorporate higherorder effects in heat transport. Its physical meaning is not necessary a priori. However, in order that the reader may set some intuitive feeling of it, it is worth saying that Q i j may be interpreted as the flux of the heat flux (see Refs. [20,25]) in solids, as the pressure tensor in fluids (see Refs. [20,30]), or as the gradient of the heat flux, but here we leave open its meaning since it could also have a structural information about the particular material. We assume that Q i j contributes to the entropy and the entropy flux. The entropy flux must be zero if q i and Q i j are zero, that is in local thermodynamic equilibrium in the absence of heat flux. Therefore, its most general form can be given as where the b i j and B i jk constitutive functions are the Nyíri multipliers, that conveniently represent the deviation from the local equilibrium form of the entropy flux , like their quadratic form in the entropy density [49]. This can be expressed also in an additive form, as the K vector of Müller, [50], if K i = (b i j − δ i j /T )q j + B i jk Q jk .
Expanding the entropy function s(e, q i , Q i j ) up to second-order approximation around a local equilibrium state, we obtain s(e, q i , Q i j ) = s (eq) (e) − 1 2ρ m i j q i q j − 1 2ρ The coefficients m i j and M i jkl have the following symmetries m i j = m ji , M i jkl = M kli j .
Note that (3) and (4) are valid for anisotropic systems too. For isotropic systems, m i j and M i jkl in (4) would reduce to a scalar and the three scalar components conjugate to the three scalar invariants of tensor Q i j , respectively. Thermodynamic stability requires that the inductivity tensors, m i j , M i jkl (see in [19,51]), are positive definite and we assume that they are constant. The entropy production σ (s) , formed by combining (2), (3), and (4), is Inequality (5) expresses the second law of thermodynamics. Following the procedures of non-equilibrium thermodynamics, we obtain the following general three-dimensional anisotropic linear relations between the thermodynamic fluxes b i j − 1 T δ i j , b ji, j − m i jq j , B i jk , B ki j,k − M i jklQkl and forces q i , q j,i , Q i j , Q jk,i b ji, j − m i jq j = L (1) i j q j + L (1,2) i jk q j,k + L (1,3) i jk Q jk + L (1,4) i jkl Q jk,l (6) b i j − 1 T δ i j = L (2,1) i jk q k + L (2) i jkl q k,l + L (2,3) i jkl Q kl + L (2,4) i jklm Q kl,m (7) B ki j,k − M i jklQkl = L (3,1) i jk q k + L (3,2) i jkl q k,l + L (3) i jkl Q kl + L (3,4) i jklm Q kl,m (8) B i jk = L (4,1) i jkl q l + L (4,2) i jklm q l,m + L (4,3) i jklm Q lm + L (4) i jklmn Q lm,n .
Here the conductivity tensors, L (α,β) , and L (γ ) , are restricted by material symmetries and by the second law. Furthermore, reciprocity relations are also to be considered, as we do in Sect. 3.

Onsager reciprocity relations
There are two different justifications of Onsager reciprocity. These are the assumptions regarding microscopic and macroscopic reversibility [52]. The concept of microscopic reversibility goes back to Onsager, [53,54], and assumes a known microstructure, based on the reversal of microscopic velocities. The principle of macroscopic reversibility assumes a particular parity of the physical quantities regarding time reversal, which is originated in the consistency of the balances, and constitutive equations with a time reversal operation [33,55]. Then the physical quantities with even parity are called αand with odd parity as β-type variables. For example, density, entropy, energy, and all thermostatic state variables are of α-type, the velocity, heat flux, and entropy flux are β-type, as one can see from the balances because time derivative changes the parity of the fields (e.g. the time derivative of an α-type variable becomes β-type), but the gradient does not. It is generally assumed that, if the thermodynamic forces are of the same type, then the conductivity tensor is symmetric and when they are of the opposite, the conductivity tensor becomes antisymmetric.
Several theoretical and experimental results support that internal variable-related thermodynamic fluxes and forces do not have definite parities, and both symmetric and antisymmetric parts of the conductivity tensors can be observed [56][57][58]. This is understandable because nothing is assumed about the microscopic structure of the material nor on the physical meaning of Q i j in NET-IV [59]. Therefore, the microscopic reversibility conditions of Onsager cannot be applied, and concept of macroscopic reversibility is not violated, if we assume that the internal variable, Q i j , does not have parity. In the following, we start with the general case, without Onsagerian reciprocity and without any assumption on the parity of the Q i j . Then we investigate the parities separately with symmetric and antisymmetric conductivity tensors. Let us remark that comparison with Extended Thermodynamics identifies Q i j as a pressure tensor or as a flux of the heat flux [43]. In this case, it must have an even character, also because the entropy flux J i and the heat flux q i are β-type, odd quantities.
With the positive sign if Q i j is β-type and with a negative one if it is α-type quantity, the sign changes only if the parity of the respective thermodynamic forces changes, too. Therefore, the L (γ ) tensors, that is the diagonal hypertensors in (6)-(9), do not change the sign.

General isotropic case without assumption on the parity of Q
In the general isotropic case, in which the symmetry properties of the body under consideration are invariant with respect to all rotations and to inversion of the frame of axes, but in which Onsager reciprocity relations are not yet imposed, we have [60] L (1,4) i jkl ≡ L (1,4) L (3,2) i jkl ≡ L (3,2) i jkl = L L (4) i jklmn ≡ L (4) i jklmn = L (4) The coefficients appearing in the entropy (4) are (15) and (16). Furthermore, in the isotropic case (where the symmetry properties of the considered body are invariant only with respect to all rotations of the frame of axes) the third-and fifth-order tensors keep the form L i jk = L ∈ i jk and L i jklm = A 1 ∈ i jk δ lm + A 2 ∈ i jl δ km + A 3 ∈ i jm δ kl + A 4 ∈ ikl δ jm + A 5 ∈ ikm δ l j + A 6 ∈ ilm δ jk , respectively, where ∈ i jk denotes the Levi Civita tensor and the quantities L and A i , i = 1, . . . , 6, are the independent components of the tensors L i jk and L i jklmn , that vanish when there is also the invariance of the properties with respect to the inversion of the axes. Thus, we obtain L (1,2) i jk = L (1,3) i jk = L (2,1) i jk = L (3,1) i jk = 0, L (2,4) i jklm = L (3,4) i jklm = L (4,2) i jklm = L (4,3) i jklm = 0.

Onsager symmetry
Now, we tentatively require Onsager reciprocity relations (10)- (14), as additional restrictions on the coefficients, and explore which further reduction this implies on the number of independent conductivity parameters Then, from (11) 2 L (1,4) i jkl = L (4,1) jkli , and being we obtain L Furthermore, for each isotropic four tensor L i jkl we have the following symmetry relation because of where T 1 , T 2 , and T 3 indicate the independent components of L i jkl . Taking into account the property (33), Onsager relations (12) 1 and (13) 2 are verified in the isotropic case and from (12) 2 we derive L (2,3) i jkl = L (3,2) kli j = L (3,2) i jkl , from which we have L Then, from (24) we obtain L (4) lmni jk = L (4) Adding (24) and (36), using Onsager relation (14) 2 and dividing by 2, we have L (4) i jklmn = C where C (4) 10 , Thus, from relation L (4) i jklmn = L (4) lmni jk the significant components of the isotropic tensor L i jklmn reduce from 15 to 11. Therefore, in case that Onsager reciprocity is imposed, from relations (32), (35), and (37) the number of conductivity parameters is reduced altogether from 34 to 24.

Entropy production
In the general isotropic case, with the aid of relations (27)-(30), (15)-(23), (25) and (26), entropy production (5) can be written as ik q i q k + L (2) i jkl q j,i q k,l + L (3) i jkl Q i j Q kl + L (4) i jklmn Q jk,i Q lm,n + L (1,4) i jkl + L (4,1) l jki q i Q jk,l + L (2,3) i jkl + L (3,2) kl ji In the case where the internal variable Q i j has odd parity, using Onsager relations and (35) and (37), expression (41) takes the form i jklmn Q jk,i Q lm,n + L (1,4) i jkl ± L (1,4) il jk q i Q jk,l + L (2,3) or in extended form From (43), it is seen that the entropy production is a non-negative bilinear form in the components of the heat flux and its gradient, and in the components of the internal variable and its gradient (see in Appendix its matrix representation σ (s) = X α L αβ X β , with X α , X β and L αβ suitable matrices).
The following inequalities can be obtained for the components of the phenomenological tensors, resulting from the fact that all the elements of the main diagonal of the symbolic matrix {L αβ } associated to the bilinear form (43) must be non-negative, representing a condition (only necessary) for the semi-definiteness of the matrix {L αβ } (see Appendix) 3 + 2C 6 + C C (4) 10 ≥ 0, C Relations (46)-(48) come from the non-negativity of the elements of the main diagonal of the sub-matrix L (4) pi jlmn . Moreover, other relations can be obtained from the non-negativity of the major minors P r (r = 1, . . . , 48) of {L αβ }, coming from Sylvester's criterion, that represents a necessary and sufficient condition for the semidefiniteness of the matrix {L αβ }. For instance, the calculation of the major minors up to sixth-order gives the relations (44) 1 , (44) 2 and (45) 1 . The non-negativity of the seventh-order major minor of {L αβ } with 3 , gives the new relation and so on. In Appendix, we give a two-dimensional form of the conductivity matrix {L αβ }, in terms of which the calculation of the conditions of positive definiteness is straightforward.

Rate equations for q i and Q i j in the general case without assumption on the parity of Q i j
Changing indexes i and j in (28), deriving it with respect to x j and substituting it into (27), we deduce where Equation (51) can be written as follows where τ = m L (1) being τ the relaxation time of the heat flux (that, then, has a finite velocity of propagation), λ the heat conductivity and l i have dimension of square length.
In analogous way, if we change i → k, j → i, k → j in equation (30), deriving it with respect to x k and inserting it into (29), we have i.e. where and τ 1 , τ 2 and τ 3 have time dimension. In the rate equations (53) and (57) 24, independent coefficients appear. These equations are the full threedimensional versions of the one-dimensional equations (12)-(13) in [28]. They represent the generalized ballistic-conductive heat transport laws in three-dimensional isotropic materials. Equation (57) can be rewritten by means of three rate equations, splitting the second-order tensor Q i j into its orthogonal components, i.e. where From Eq. (57), we derive the rate equations for Q, Q i j , and The rate equation for Q is (i = j) i.e. where being τ 0 the relaxation time of Q; the rate equation for Q i j is where

The rate equations for q i and Q i j with Onsager reciprocity in the case where Q i j has odd parity
In Sect. 6, we have obtained the rate equations for q i and Q i j and for the scalar part, the deviator of the symmetric part and the skew-symmetric part of Q i j (see (53), (57) or (53) and (68), (71), (74), respectively) without assuming reciprocity relations, but only isotropy. In this section, we derive the heat transport laws in three-dimensional isotropic materials (57), (68), (71), (74) in the form (85), (86), (88) and (90), by using Onsager reciprocity relations (32)- (37), that reduce the number of coefficients in these rate equations from 24 to 21 (when compared to the general isotropic case). The phenomenological equations (27) and (28) remain unchanged (thus also the rate equation (53)), but equations (29) and (30) assume the following form where, with respect to (29) and (30) where The rate Eq. (79) is the same as (57), but with L i replaced by C i (i = 1 . . . 9). We remark that the coefficients l 21 , l 31 , and l 41 in (79) transform according to Onsager relations (32) and (35), so that we have By virtue of (83) and (80) 2 , (81) 1 and (82) 3 , the further conditions for the coefficients are worked out: Thus, using relations (84), the rate equation (79) for Q i j takes the form As in (63), we split the second-order tensor Q i j in its orthogonal components Q, Q i j and Q [i j] , its scalar part, the deviator of its symmetric part, its skew-symmetric part [see (64)- (66)] that, having Q i j odd parity, have also odd parity. In the following, we work out the rate equations for Q, Q i j and Thus, from Eq. (85) we derive: the rate equation for Q (obtained when i = j) where τ 0 is given by (69) 1 and where ∧ τ is given by (72) 1 and where ∨ τ is given by (75) 1 and 6.1 One-dimensional heat transport in the case where Q i j has odd parity In this subsection, we focus on the one-dimensional case, in order to appreciate how the generalization from one dimension to three dimensions analysed in this paper is far from trivial. In the one-dimensional case, we have that the components of B and Q reduce to The system of Eqs. (27)-(30) (in which we use the Onsager relations assuming that Q i j has odd parity) becomes where m > 0, M > 0 (see [28]) and with (·) ,x indicating the derivative of (·) with respect to x. We observe that the system of Eqs. (93)-(96) obtained here is more general of equations (7)- (10) deduced in [28], because of the presence of the phenomenological constant L (1,4) in (93) and (96) and the fact that the coefficients L (1,4) , L (2) , L (2,3) , M, L (3) , C (4) have been obtained from a three-dimensional approach.
In this case, the entropy production (43) assumes the form with C (4) = L 111111 (see matrix (156) of Appendix), or in symbolic matrix notation Because the bilinear form (101) must be non-negative, the matrix A (that is symmetric) associated to this form is non-negative semi-definite, so that the elements of its main diagonal and its major minors must be non-negative Using (94) and (96), Eqs. (93) and (95) become where D = L (1,4) − L (2,3) . In the following, we introduce the relaxation time of the internal variable Q, called τ J : Furthermore, we have supposed the body is at rest, so that material derivative coincides with the partial time derivative (·) ,t . Equations (105) and (106) are analogous to equations (12) and (13) of [28]. For L (2) = C (4) = 0, these equations coincide with those provided in [20] or [25] by assuming Q i j as the flux of the heat flux.
In the following, we will derive heat transport equations analogous but more general of that obtained in [28], where the finite speed of thermal disturbances and the ballistic and diffusive motion of phonons (heat carriers) are taken into account. Instead, in Fourier equation the velocity of heat propagation is infinite. Differentiating Eq. (105) with respect to time, Eq. (106) with respect to the spatial variable x and using Eq. (105) and its second spatial derivative, we can eliminate Q and work out the following generalized ballistic-conductive heat transport law where Equation (108) has been obtained via several differentiations of the linear governing Eqs. (105) and (106). Hence, Eq. (108) is not equivalent to the system of Eqs. (105) and (106). In fact, (108) has a larger set of solutions, coming from the larger number of necessary initial conditions. Thus, we derive where H is defined by (109) and λ is given by (54) 2 . In (110), we see that relaxation time τ q = τ + τ J is given by two contributions: the first comes from the relaxation time of the heat flux [see (54) Thus, we can write where Guyer-Krumhansl equation. In the case where C (4) = M = 0, the heat Eq. (108) becomes then, we work out where l, having the dimension of a length which may be interpreted as an average, mean free path of the heat carriers (phonons), i.e. the average length between successive collision amongst them. We observe that only in Guyer-Krumhansl heat equation the coefficient multiplying the field q ,x x has the physical meaning of l 2 .
Cahn-Hilliard-type equationIn the case where C (4) = M = m = 0, the heat Eq. (108) becomes from which we obtain Jeffreys-type equation(or double-lag model [61]) In the case where C (4) = L (2) = m = D = 0 (then, H = 0), the heat Eq. (108) becomes thus we derive: We note that in the Jeffreys-type heat equation τ J is the relaxation time of q.
What is specially worth in this subsection is not only the ability to obtain many situations studied up to now, but specially the fact that the coefficients appearing in the one-dimensional case are complicated combinations of the independent coefficients appearing in the three-dimensional case. Thus, measurements in one dimension are not sufficient to give information in the general three-dimensional situation, which is the only one able to exhibit the basic meaning of each coefficient. We emphasize that Jeffrey type, Maxwell-Cattaneo-Vernotte and Fourier equations are the same as in [28].

Rate equations for q i and Q i j in the isotropic case where Q i j have even parity
In Sect. 6, we have obtained the rate Eqs. (53) and (57) for the heat flux q i and the internal variable Q i j , respectively, in the general isotropic case without assumptions regarding the parity of the internal variable Q i j (q i is odd and Q i j can be of odd or even type) and then we have not discussed Onsager reciprocity relations. In Sect. 7, we have shown how these rate equations transform supposing the odd parity of Q i j . In this section, we treat the case where Q i j has even parity and it is very easy to see that the rate Eq. (53) remains unchanged (as in the odd parity case). Instead, the rate Eq. (57) (that takes the form (85) when we assume the even parity of Q i j ) transforms in where the quantities l 21 , l 31 and l 41 take the following form in which Onsager symmetry relations (32) and (35) have been applied, with negative sign. The other coefficients continue to have the same definitions given in Sect. 6, but the coefficients of q i, j and q j,i have different signs with respect to those in (85). Furthermore, in this considered case relations (84) 1,2 become Finally, the rate equations for the orthogonal components Q, Q i j and Q [i j] , that are still of even type, remain formally unchanged from Eqs.
7.1 One-dimensional isotropic heat transport in the assumption that Q i j has even parity Taking into account expressions (92), using the Onsager relations (10)- (14) in the case where the internal variable Q i j has an even parity, equations (27)-(30) take the form where only Eqs. (134) and (135) are different from (95) and (96) because of the signs of the first terms in their right-hand sides. As consequence of this difference, we have that the entropy production (101) takes the new reduced form so that the associated matrix A (that is diagonal and thus symmetric) takes the following diagonal form so that only relations (103) are still true. Furthermore, using (133) and (135), Eqs. (132) and (134) become where D = L (1,4) − L (2,3) . Relation (138) is equal to (105), while (139) has opposite sign in its right-hand side with respect to (106) (D continues to have the same value). Finally, deriving Eq. (138) with respect to time, Eq. (139) with respect to the spatial variable x and using Eq.
(138) and its second spatial derivative, we can eliminate Q and work out the same the heat transport Eq. (108) (and than (110)) where the only difference consists in the fact that the quantity H defined by (109) takes the new form Thus, H is always positive in the case of even parity of Q i j .

Special cases of heat transport equation in the assumption that Q i j has even parity
Applying the procedures used in Sect. 7.2 to obtain from (108) (and also (110)) special cases, we derive the following results: a) The ballistic-conductive heat transport Eqs. (113) and (114), being C (4) = L (2) = 0 and H = L (3) L (2) + D 2 , take the following form different from (113) and (114)

Discussion and conclusions
In this paper, ballistic-conductive heat transport in isotropic materials has been treated in the framework of Nonequilibrium Thermodynamics with Internal Variables (NET-IV). Onsager reciprocity has also been considered (in both particular cases in which Q i j have been assumed to be odd or even with respect to macroscopic time reversal) and the consequences were derived. For the sake of fast applicability, the explicit expressions for the components of the conductive matrix are given in Appendix in the two cases. 1 Our approach, NET-IV, is general and universal. It characterizes the deviation from local equilibrium both in the entropy density and in the entropy flux in the simplest possible functional forms. The entropy density depends on the internal variables quadratically, in order to preserve the concavity, that is thermodynamic stability. The entropy flux depends on the internal variables linearly; therefore, it disappears when they are zero. As long as these two physical conditions and the entropy inequality are valid, the derived consequences are also valid. The generality of the assumptions ensures the universality of the final evolution equations. Here we have considered a strictly linear theory, when the m and M tensors and the conductivity tensors, L (α,β) and L (γ ) , are constant.
The conditions of positive definiteness of the corresponding conductivity matrix can be calculated directly with the help of computer algebra programs. Though the expressions are very cumbersome, it should be noted that every coefficient appearing in them corresponds in principle to an observable phenomenon. Instead, the much simpler one-dimensional case may grasp essential qualitative features, but its coefficients are a combination of three-dimensional coefficients giving a deeper and more complete description. We have obtained a complete set of equations for generalized ballistic-conductive heat transport in three-dimensional isotropic rigid conductors for the variables T, q i , Q i j . These are the balance of internal energy (1) with the caloric equation of state s eq (e) = 1/T and the balance-type constitutive Eqs. There are two different aspects of ballistic heat transport in continua. From the point of view of kinetic theory, it is the propagation of phonons without collisions with the lattice. Then heat is reflected only at the boundaries of the medium. This microscopic understanding is the foundation of the so-called ballistic-diffusive integrodifferential model of Chen [64][65][66][67][68]. The kinetic theory and macroscopic considerations are mixed, the distribution function f is split into two parts, one for ballistic phonons and the other referred to diffusive phonons. Also, internal energy and heat flux are decomposed into ballistic and diffusive components. This approach leads to two independent continuum representations. First, it is a particular boundary condition for continuum theories that can also be introduced to second-sound models, like Guyer-Krunhansl equation [69]. On the other hand, for ballistic phonons, the speed of propagation is equal to the speed of 'first' sound, the speed of elastic waves in the medium. The speed of propagation is independent of the boundary conditions in a continuum approach, and this is the meaning of the ballistic terminology in our theory, following Rational Extended Thermodynamics (RET) [22,70]. It is also remarkable that Chen's model is equivalent to an extended continuum heat transport theory, where the coexistence of two kinds of heat carriers (ballistic and diffusive phonons) is assumed as it was shown by Lebon et al. [71,72] and investigated in [73].
Theories of Extended Thermodynamics (ET) assume that the constitutive equations are local, and the rate equations are written in a hierarchical series of balances, where the dissipative fluxes appear as densities in the consecutive balance. These assumptions are consequences of the definition of the macroscopic fields as moments of the single-particle phase space probability density and the Boltzmann equation. In our case, with internal variables, this structure is the consequence of the second law and can be observed on the left-hand side of (6) and (8). Then essential aspects of ET are well represented. On the other hand, NET-IV has many material coefficients that are missing in ET, in particular in Rational Extended Thermodynamics, where only the two relaxation times of the Callaway collision integral represent the material properties. This property of RET is attractive, but the price is not only that the validity of the theory is connected to the particularities of the microscopic model, but also that the speed of the ballistic propagation, the speed of elastic waves, can be obtained exactly only by considering the complete moment series, or practically by using dozens of evolution equations (with consecutively increasing tensorial orders) [22]. The low number of material coefficients leads to many evolution equations in modelling ballistic propagation of heat.
Giving the three-dimensional structure of ET and NET-IV for heat transport in case of isotropic materials opens the field to build and solve realistic models of two-and three-dimensional experimental setups, where the two theories lead to different predictions. To appreciate some of the original aspects of this work, let us eventually comment Eqs. (28) and (30) for b i j and B i jk and their consequences on the entropy flux. It is well known in the literature [25,74,75] that one of the expression of the entropy flux is Note then that the constitutive equation for the entropy flux (3), when b i j and B i jk are given by (28) and (30), lead to a richer expression than (143), namely Thus, in our analysis the extended entropy flux is more general than (143) and plays an important role in the thermodynamic consistency of couplings with the heat flux q i and tensorial internal variables as Q i j .

Appendix
Here, we give a two-dimensional symmetric explicit representation of the conductivity matrix {L αβ }. This form is useful when the conditions of positive definiteness have to be calculated. Though the explicit writing is cumbersome, it is especially useful when an abstract notation is not sufficient, but explicit calculations must be done, or when a computer program for solving equations or carrying out numerical simulations must be implemented.