Calculations for Extended Thermodynamics of dense gases up to whatever order and with all the symmetries

The 14 moments model for dense gases, introduced in the last years by Arima, Taniguchi Ruggeri, Sugiyama, is here considered. They have found the closure of the balance equations up to second order with respect to equilibrium; here the closure is found up to whatever order with respect to equilibrium, but for a more constrained system where more symmetry conditions are imposed and this in agreement with the suggestion of the kinetic theory. The results, when restricted at second order with respect to equilibrium, are the same of the previously cited model but under the further restriction of full symmetries.


Introduction
Starting point of this research is the article [1] which belongs to the framework of Extended Thermodynamics. Some of the original papers on this subject are [2], [3] while more recent papers are [4]- [18] and the theory has the advantage to furnish hyperbolic field equations, with finite speeds of propagation of shock waves and very interesting analytical properties. It starts from a given set of balance equations where some arbitrary functions appear; restrictions on these arbitrariness are obtained by imposing the entropy principle and the relativity principle. However, these restrictions were so strong to allow only particular state functions; for example, the function p = p(ρ, T ) relating the pressure p with the mass density ρ and the absolute temperature T , was determined except for a single variable function so that it was adapt to describe only particular gases or continuum. This drawback has been overcome in [1] and other articles such as [19]- [34] by considering two blocks of balance equations; for example, in the 14 moments case treated in [1], they are where The first 2 components of P N are zero because the first 2 components of equations (1) 1 are the conservation laws of mass and momentum; the first component of Q E is zero because the first component of equations (1) 2 is the conservation laws of energy. The whole block (1) 2 can be considered an "Energy Block". The equations (1) can be written in a more compact form as where F A = (F N , G E ) , F kA = (F kN , G kE ) , P A = (P N , Q E ) .
In the whole set (2), F A are the independent variables, while F kij , G ki , P ij , Q i are constitutive functions. Restrictions on their generalities are obtained by imposing 1. The Entropy Principle which guarantees the existence of an entropy density h and an entropy flux h k such that the equation holds whatever solution of the equations (2). Thanks to Liu' s Theorem [35], [36], this is equivalent to assuming the existence of Lagrange Multipliers µ A such that An idea conceived by Ruggeri is to define the 4-potentials h ′ , h ′k as so that eqs. (4) 1,2 become which are equivalent to if the Lagrange Multipliers are taken as independent variables. A nice consequence of eqs. (6) is that the field equations assume the symmetric form. Other restrictions are given by 2. The symmetry conditions, that is the second component of F N is equal to the first component of F kN , the third component of F N is equal to the second component of F kN , the second component of G E is equal to the first component of G kE . Moreover, F ij , F kij and G ki are symmetric tensors. The symmetry of F kij and G ki is motivated by the kinetic counterpart of this theory (see section 4 of ref. [19]), even if was not imposed in [1] in order to have a more general model. We propose, in a future article, to remove this further constraint. Thanks to eqs. (6) these conditions may be expressed as where we have assumed the decomposition µ A = (µ, µ i , µ ij , λ, λ i ) for the Lagrange Multipliers. Moreover µ ij is a symmetric tensor. The next conditions come from 3. The Galilean Relativity Principle.
We prefer to devote an entire section, the next one, to describe how to impose this principle. The result, combined with the above conditions coming from the Entropy Principle and the Symmetry Conditions will be that a scalar function H exists, such that In section 3 of the present article, we find the general solution up to whatever order with respect to equilibrium, of the conditions (9), (10). In section 4 we will see the implications of this solution to a second order theory, coming back to the moments as independent variables; this will allow us to recover, as first order approximation for the functions F kij and G ki , the same result of [1]; similarly, for the second order approximation for the entropy density and its flux. So the present results generalize that article and also confirm it, because we guarantee that the equations at subsequent orders don' t give other restrictions on the first order theory.

The Galilean Relativity Principle
There are two ways to impose this principle. One of these is to decompose the variables F A , F kA , P A , µ A in their corresponding non convective partsF A ,F kA ,P A ,μ A and in velocity dependent parts, where the velocity is defined by This decomposition can be written as After that, all the conditions are expressed in terms of the non convective parts of the variables. This procedure is described in [2], [36] for the case considering only the block (1) 1 and is followed in [1] for the whole set (1). Another way to impose this principle leads to easier calculations; it is described in [37] for the case considering only the the block (1) 1 and here we show how it is adapt also for the whole set (1). First of all, we need to know the transformation law of the variables between two reference frames moving one with respect to the other with a translational motion with constant translational velocity v τ . To know it, we may rewrite (12) in both frames, that is where the index a denotes quantities in the absolute reference frame and index r denotes quantities in the relative one;F B ,F kB ,ĥ ′ ,ĥ ′k ,μ B haven't the index a, nor the index r because they are independent from the reference frame. Now we can use a property of the matrix X A B ( v) which is a consequence of its definition (13) and reads So we may contract (14) 6,7 with X C A (− v r ) so obtaininĝ which can be substituted in (14) 1,2 . The result is Now we use another property of the matrix X A B ( v) which is a consequence of its definition (13) and reads Moreover, we use the well known property In this way the equations (16) become Finally, we deduceĥ ′ ,ĥ ′k andμ A from (14) 8,9,10 and substitute them in (14) 3,4,5 so obtaining where for the last one we have also used a contraction with (20) and (21) give the requested transformation law between the two reference frames and it is very interesting that it looks like eqs. (12). Now, if the Lagrange Multipliers are taken as independent variables, eqs. (21) 3 are only a change of independent variables from µ a B to µ r C , while (19) 1 , (20), (21) 1,2 are conditions because they involve constitutive functions where the form of the functions F A , F kA , h ′ , h ′k don' t depend on the reference frame for the Galilean Relativity Principle. If we substitute µ a B from eq. (21) 3 in (22) 1−4 and then substitute the result in (19) 1 , (20), (21) 1,2 , we obtain Well, these expressions calculated in v i τ = 0 are nothing more than eqs. (22) 5−8 , as we expected. But, for the Galilean Relativity Principle they must be coincident for whatever value of v i τ ; this amounts to say that the derivatives of (23) with respect to v i τ must hold. This constraint can be written explicitly more easily if we take into account that µ r which can be written explicitly by use of (13) and reads Consequently, the derivatives of (23) 3,4 with respect to v i τ become where we have omitted the index a denoting variables in the absolute reference frame because they remain unchanged if we change v i τ with −v i τ , that is, if we exchange the absolute and the relative reference frames. It is not necessary to impose the derivatives of (23) 1,2 with respect to v i τ because they are consequences of (26) and (6). Consequently, the Galilean Relativity Principle amounts simply in the 2 equations (26). So we have to find the most general functions satisfying (7) and (26). After that, we have to use eqs. (6) 1 to obtain the Lagrange Multipliers in terms of the variables F A . By substituting them in (6) 2 and in h ′ = h ′ (µ A ), h ′k = h ′k (µ A ) we obtain the constitutive functions in terms of the variables F A . If we want the non convective parts of our expressions, it suffices to calculate the left hand side of eqs. (6) 1 in v = 0 so that they becomê From this equation we obtain the Lagrange Multipliers in terms ofF A (Obviously, they will beμ A ) and after that substitute them in h ′ = h ′ (µ A ), h ′k = h ′k (µ A ) (the last of which will in effect beĥ ′k ) and intoF kA = ∂h ′k ∂µ A , that is eq. (6) 2 calculated in v = 0. It has to be noted that from (11) it followsF i = 0, so that one of the equations (27) is 0 = ∂h ′ ∂µ i ; this doesn' t mean that h ′ doesn' t depend on µ i , but this is simply an implicit function defining jointly with the other equations (27) the quantitiesμ A in terms ofF A . We note also here the ground to settle µ i = 0 at equilibrium: in fact, in this state we have µ ij = 0, λ i = 0 so that, for the Representation Theorems, ∂h ′ ∂µ i is proportional to µ i and ∂h ′ ∂µ i = 0 implies µ i = 0. By using a procedure similar to that of the paper [37], we can prove that we obtain the same results of the firstly described approach. Now, from (7) 2 it follows ∂h ′[i ∂µ j] = 0; this equation, together with (7) 1 are equivalent to assuming the existence of a scalar function H such that the above mentioned (8) holds. In fact, the integrability conditions for (8) are exactly (7) 1 and ∂h ′[i ∂µ j] = 0. Thanks to (8), we can rewrite (7) and (26) as the above mentioned eq. (9) and We note now that the derivative of (28) 1 with respect to µ k is equal to the derivative of (28) 2 with respect to µ; similarly, the derivative of (28) 1 with respect to λ k is equal to the derivative of (28) 2 with respect to λ, as it can be seen by using also eqs. (9). Consequently, the left hand side of eq. (28) 1 is a vectorial function depending only on two scalars µ, λ and on a symmetric tensor µ ij . For the Representation Theorems [38]- [46], it can be only zero.

The solution up to whatever order
Let us firstly show a particular solution of (9) and (10). Let ψ n (µ, λ) be a family of functions constrained by Let us define the function In Appendix 1 we will show that eqs. (9) and (10) are satisfied if H is replaced by H 1 ; in other words, H = H 1 is a particular solution of our conditions. Moreover, (H 1 ) eq. = ψ 0 (µ, λ) which is an arbitrary two-variables function, such as H eq. . So we can identify and define In this way the conditions (9) and (10) become and we have also (∆H) eq. = 0 .
We now note that from (33) 3,4 we deduce that the function ∂∆H ∂µ k has all the derivatives with respect to µ i , µ ij , λ i which are symmetric tensors, so that its Taylor ' s expansion around equilibrium is where H is a symmetric tensor depending only on the scalars µ and λ, so that it has the form where I p denotes the set of all non negative integers r such that r + p is even. But also ∂∆H ∂µ has all the derivatives which are symmetric tensors; in fact, the derivatives of (33) 3,4 with respect to µ are The derivatives of (33) 1,2 with respect to µ ab are whose skew-symmetric parts with respect to b and i are thanks to eq. (33) 3,4 . By using (36) we obtain that also ∂H ∂µ has all the derivatives which are symmetric tensors, so that its expansion around equilibrium is where the derivative ofH q,r has been introduced for later convenience and without loss of generality. By integrating (39) we obtain But the function ∆H is present in (33), (8) only through its derivatives with respect to µ and µ k so that the functionH doesn't effect the results and, consequently, it is not restrictive to assume that it is zero. By substituting (40) in (36) we see that the function ∆H has derivatives which are symmetric tensors and we have also its expansion, that is where we have defined H 0,q,r =H q,r ; in this way, the term of (41) with p = 0 is (40), while for the other terms we can change index according to the law p = 1 + P and obtain the remaining part of (36). Now an interesting consequence of (33) and (34) is that PROPERTY 1: " The expansion of ∆H up to order n ≥ 1 with respect to equilibrium is a polynomial of degree n − 1 in the variable µ." Thanks to this property, it is not restrictive to assume for ∆H a polynomial expansion of infty degree in the variable µ; we simply expect that the equations will stop by itself the terms with higher degree. Consequently, it is not restrictive to assume that the functions H p,q,r in (41) can be written as So, even if µ is not zero at equilibrium, for what concerns ∆H, we can do an expansion also around µ = 0; obviously, the situation is different for the particular solution H 1 reported in eq. (30). We report the proof of Property 1, in Appendix 2.
If we substitute (41) and (42) in (33) 3,4 , we obtain identities. If we substitute (41) and (42) After that, we see that (43) 1 is satisfied as a consequence of (44). Let us focus now our attention to eq. (43) 2 ; for p = 0, 1 it becomes where, for (45) 2 we have used (44) with p = 2. After that, eq. (43) 2 with use of eq. (44) • in the case with p even, gives (45) 1 with (q + p 2 , s + p 2 ) instead of (q, s), • in the case with p odd, gives (45) 2 with (q + p−1 2 , s + p−1 2 ) instead of (q, s). But we have now to impose (33) 5 . To this end, let us take its derivatives with respect to µ i 1 , · · · , µ i P , µ h 1 k 1 , · · · , µ h Q k Q , λ j 1 , · · · , λ j R and let us calculate the result at equilibrium; in this way we obtain where overlined indexes denote symmetrization over those indexes, after that the other one (round brackets around indexes) has been taken. (Note that, in the fourth term the index R − 1 appears; despite this fact, the equations holds also for R = 0 but in this case, this fourth term is not present as it is remembered also by the factor R). Now, the first, second, fifth and sixth term can be put together so that the above expression becomes This relation, for P = 0, 1 reads with the agreement that the last terms are not present in the case R = 0. (For eq. (47) 2 we have used (44) with p = 2). For the other values of P , eq. (46) with use of eq. (44) • in the case with p even, gives (47) 1 with (q + p 2 , s + p 2 ) instead of (q, s), • in the case with p odd, gives (47) 2 with (q + p−1 2 , s + p−1 2 ) instead of (q, s). Summarizing the results, we have that (44) gives H P,Q,R,s in terms of H 0,Q,R,s and H 1,Q,R,s , while eqs. (45) and (47) give restrictions on H 0,Q,R,s and H 1,Q,R,s . Let us now see how to solve these restrictions. First of all, eq. (34) expressed in terms of (41) and (42) becomes We have also already imposed that ∆H becomes zero when calculated in µ = 0 and µ k = 0. This can be expressed in terms of (41) and (42) as From (47) 1 with R = 0, by using (48) and (49), we find Let us now apply (45) 2 with (r − 1, s − 1) instead of (r, s) and, subsequently, (45) 1 , so obtaining where the values of k are restricted by 2k ≤ r and k ≤ s. So, if s ≤ r−1 2 ( we recall that r is odd because the sum of the first and third index must be an even number), we can take k = s so obtaining H 1,q,r,s in terms of H 1,q,r,0 ; if s ≥ r+1 2 we take k = r−1 2 so obtaining H 1,q,r,s in terms of H 1,q+ r−1 2 ,1,s− r−1 2 which is zero for (45) 2 and (50). So we have obtained that By applying this result, (45) 1 becomes where, for the value with r = 0 we have used (50). So every thing is determined in terms of H 1,q,r,0 which remains arbitrary up to now. After that (45) 1,2 become consequences of (52) and (53). Regarding (47) 1 , we have already imposed it for R = 0; if R ≥ 2, s ≥ R 2 it is an identity, as a consequence of (52) and (53). But if R ≥ 2, s ≤ R−2 2 by using (52) and (53) it becomes for s = 0, while for s ≥ 1 it becomes and this is ∂ 2s ∂λ 2s of the relation obtained from (54) with (Q + s, R − 2s) instead of (Q, R). Regarding (47) 2 , for R ≥ 3, s ≥ R−1 2 it is an identity, as a consequence of (52) and (53). Similarly, if s ≤ R−3 2 , R ≥ 3 by using (52)  and this is ∂ 2s+1 ∂λ 2s+1 of the relation obtained from (54) with Q = q + s + 1, R = r − 2s − 1. Finally, in the case R = 1 by using (52), (53) and (50) it becomes an identity. We note that (54) can be written as so that it allows to determine, with an integration, the function H 1,Q+1,R−1,0 in terms of H 1,Q,R−1,0 . So H 1,0,R−1,0 remains arbitrary and also the family of constants arising from the integrations.

A second order theory
For a second order theory, we need the expressions of h ′ and h ′k up to third order. By using the results of the previous sections, we can see that they are determined in terms of ψ 1 , H 1,0,1,0 , H 1,1,1,0 , H 1,2,1,0 , H 1,0,3,0 and they are h ′ = h 0,0,0 + h 0,1,0 µ ll + Here ψ 1 (µ, λ), H 1,0,1,0 (λ), H 1,0,3,0 (λ) are arbitrary functions while H Now we want express our closure by taking the moments as independent variables. To this end we will indicate with () eq. the quantity () calculated at equilibrium when this is defined as the state with µ i = 0, µ ij = 0, λ i = 0; we will indicate with () Eq. the quantity () calculated at equilibrium when this is defined as the state with ,F ij = 0,Ĝ i = 0. Both states are coincident at the order zero.

The closure at equilibrium
Let us substitute (56) and (57) in (6) 1 and calculate the right hand sides at equilibrium and the left hand side at Equilibrium and for v i = 0. So we obtain where p is the pressure and ε the internal energy density. The same procedure applied to (6) 2 and (5) 1 givesF From the second of these expressions and from the definition h Eq. = ρs, we obtain d s − 2λ[d ε + pd (1/ρ)] = 0, where we have used h 0,1,0 = − 1 2λ h 0,0,0 which comes out from (58) 1,2 . The result, thanks to the Gibbs Relation, allows to identify as in [1] that with T the absolute temperature. We will see in the sequel that, if we take as variables describing equilibrium the absolute temperature T and the chemical potential µ Eq. , then all our closure will be determined in terms of the arbitrary functions ψ 1 (µ, λ), H 1,0,1,0 (λ), H 1,0,3,0 (λ) and of the two constants arising from integration of (60). Also p, ρ and ε will be functions depending on them; in particular, for a third order theory we will have Instead of this, if we take as variables describing equilibrium ρ and T we will see that more difficult calculations are needed. To reach this end, let us now in the remaining part of this subsection to prepare some results which will be used later. Let us use (61) 1 to find µ Eq. as a function of (ρ, T ). By taking its derivatives with respect to ρ and T , we obtain Similarly, by taking the derivatives with respect to ρ and T of we find By taking the derivatives with respect to ρ and T of ε = 1 2 we find Well, from (67) and (69) 2 we deduce By substituting in (69) 1 we obtain which is well known because from the Gibbs Relation it follows whose integrability condition is just (71). Eqs. (65), thanks to (70), now become Regarding the closure for the fluxes at first order, we have found (62) 1,2 , where the coefficient φ 0,0,1 appears. This is given by (59) 2 ; to find a relation between it and h 0,0,0 given by (58) 1 , let us take the derivative By using also (65) 1 , (68), (66), (61) 1 and (67), we obtain By comparing eqs. (62) of the present article with eqs. (47) of [1] calculated at equilibrium, we find the same result with Moreover, in [1] there was only the condition (44) 1 on β 1 (because (44) 4 can be considered as the definition of h 4 ) and this condition is exactly the present condition (74). Consequently, we also can define h 4 from Making the point of situation, at equilibrium we have 3 state functions describing the particular material under consideration: p(ρ, T ), ε(ρ, T ), β 1 (ρ, T ). From the mathematical viewpoint they are arbitrary except that the first 2 of them are linked by (71) and the third one is determined by (74), except for an arbitrary function of the single variable T . In the sequel, we will always use the identification (75).
Other expressions which will be useful in the sequel are the following ones These relations are consequences of (58), (61), (59) 2 and (75); moreover, in the passage denoted with (1) = we have used eq. (66) and (68), while in the passage denoted with (2) = we have used eq. (66). We can now proceed and note that from (70) it follows that and that the determinant D of the matrix is equal to Moreover, we have: where we have used (58) in the first passage and (77) 2 in the second passage, where we have used (58) in the first passage and (61) 1 in the second passage, where we have used (59) 3 and (58) 1 in the first passage and (66) in the second passage, where we have used (59) 1 , (58) 1 , (66), where we have used (59) 1 and (64) 1 , where we have used (59) 1 and (58) 1 in the first passage and a part of (77) 2 in the second passage. Now, from (58) 6 and (59) 2 it follows h 0,0,2 = ∂ ∂λ φ 0,0,1 . By using this result, we may find where in the passage denoted with = we have used eq. (74). So, up to now, for the first order we have found every thing in terms of p(ρ, T ), ε(ρ, T ), φ 0,0,1 (ρ, T ) and φ 0,1,1 . On the last two of these functions we have the condition (74) and the following ones where in the passage denoted with (1) = we have used (59) 4 , and (73) 1 ; moreover, in the passage denoted with (2) = we have used eq. (66) and (61) 3 .
where in the passage denoted with = we have used eq. (66) and (61) 3 . Now we note that (59) 2,4 , (89) and (90) are 4 equations in the 3 unknowns ∂ψ 1 ∂µ , ∂ 2 ψ 1 ∂λ∂µ , ∂ 3 ψ 1 ∂λ 2 ∂µ ; it follows that the determinant of the complete matrix must be zero, that is Now this determinant can be simplified by adding to its first column the second one multiplied by 2 λ and the third one multiplied by 2 λ 2 ; moreover, in line 3, column 4 we can substitute ∂ ∂λ H 1,0,1,0 from (60) 1 , so that it becomes (It is interesting that now H 1,1,1,0 disappears and the result depends on φ 0,0,1 only through ∂φ 0,0,1 ∂T ). The resulting expression can be used to find This result can be rewritten by inserting ∂φ 0,0,1 ∂T from (76), so that it becomes After that, from (87) with (76) we obtain Moreover, it will be useful for the sequel to know the determinant where in the second passage we have used (82), (81) and (93).
To conclude this preparatory work, let us note that where in the passage denoted with from which we deduce (96)

4.3
The entropy density and its flux, up to second order.
From properties of homogeneous functions, we have that the expression of first order of the entropy density is equal to h (1) = π A ∂h (1) ∂π A where π A denotes π,F <ij> , q i . Moreover, from (4) 1 it follows It follows that h (1) = π A ∂h ∂π A (0) = 0, that is the homogeneous part of the entropy density at first order is zero. By proceeding in a similar way for the entropy density at second order, we obtain where in the passage denoted with We appreciate that this result is exactly the same of eq. (55) in [1] with h 3 given in their eq. (43) 2 . Also the expression for h 2 given in eq. (43) 1 of [1] corresponds to the present (80) 2 .
For what concerns the entropy flux h k , we have that (4) 2 implies which can be calculated in v j = 0 so transforming themselves in the corresponding ones denoted with anˆ. By calculating them at equilibrium, the first two will say thatĥ k (1) doesn' t depend on π, nor on F <ab> , while the third one becomes ∂ĥ k (1) ∂q a = 1 T δ ka from which we obtain the well known expressionĥ k We note now that, by substituting φ 0,1,1 from (105) 2 in (104), we obtain By calculating (108) at first order with respect to equilibrium, thanks also to eq. (110), (103), (99), (102), (104), they become ∂ĥ k Kδ k<a q b> , which is exactly the same which can be obtained from the calculations in [1] and in the particular case with all the symmetry conditions, that is L = 5 6 K.

A comparison with the results of [1].
We have already seen that the entropy density found with the present approach has at equilibrium the same expression than in [1]; similarly, for the first order expression h (1) = 0 and for the second order expression (106). Similarly, for the non convective part of the entropy flux in both articles it is zero at equilibrium, has the form (109) at first order with respect to equilibrium, and the same form (111) at second order. Regarding the non convective part of F kij , in both papers it is zero at equilibrium and have the same expression at first order, that is (110) of the present paper and (47) 1,2 of [1]. Regarding G kill , its expression at equilibrium here is given by (62) 2 which is the same of [1], except for identifying the present function φ 0,0,1 with the function β 1 of [1], as it has be done in eq. (75). For what concerns G kill (1) , we obtain the same expression if we link the function φ 0,1,1 of the present paper with the function K of [1] through eq. (105) 2 ; the expressions are (105) 1 in the present paper and (47) 3 of [1], obviously in the case L = 5 6 K of full symmetries. It is worth noting that in (48) 2 of [1] we find also the definition of β 3 which can be also interpreted as the definition of K, that is Moreover, in (48) 1 of [1] also β 2 is defined and under the further assumption of full symmetries (that is L = 5 6 K), this definition is equivalent to which links the function φ 0,0,1 of the present article with the function β 3 of [1]. After that, • The equations (43)  Let us prove that H = H 1 , with H 1 given by (30) and ψ n constrained by (29), is a particular solution of (9) and (10).
B Appendix 2: Proof of the PROPERTY 1.
Let us the property 1 with the iterative procedure and let ∆H n denote the homogeneous part of ∆H of order n with respect to equilibrium. We have, • Case n = 1: The equation (33) 5 at equilibrium, thanks to (34) becomes 2 ∂ 2 ∆H 1 ∂µ∂µ ki λ = 0 from which we have that ∂∆H 1 ∂µ can depend only on µ, µ i , λ, λ c ; but the representation theorems show that no scalar function of order 1 with respect to equilibrium can depend only on these variables. It follows that ∂∆H 1 ∂µ = 0 so that ∆H 1 is of degree zero with respect to µ and the property is verified for this case.
• Case n ≥ 2: Let us suppose, for the iterative hypothesis that ∆H up to order n ≥ 1 with respect to equilibrium is a polynomial of degree n − 1 in the variable µ; we proceed now to prove that this property holds also with n + 1 instead of n.
In fact, eq. (33) 1 up to order n − 1 gives ∂ 2 ∆H n ∂µ∂µ ij = ∂ 2 ∆H n+1 ∂µ i ∂µ j from which we have ∆H n+1 = P n−2 + ∆H n+1 i (µ, µ ab , λ, λ c )µ i + ∆H n+1 0 (µ, µ ab , λ, λ c ) where P n−2 is a polynomial of degree n − 2 in µ and which is at least quadratic in µ j . After that, eq. (33) 2 up to order n gives which, thanks to (115) becomes This relation, calculated in µ j = 0 gives because P n−2 is at least quadratic in µ j . The derivative of (116) with respect to µ j , calculated then in µ j = 0, is from which ∂∆H n+1 j ∂λ i = P ij n−1 with P ij n−1 a polynomial of degree n − 1 in µ. By integrating this relation, we obtain ∆H n+1 i = P i n−1 + f i n−1 (µ, µ ab , λ, ) where P i n−1 is a polynomial of degree n − 1 in µ. But, for the Representation Theorems, a vectorial function such as f i n−1 is zero because it depends only on scalars and on a second order tensor. It follows that By using this result, eq. (117) can be integrated and gives ∂∆H n+1 with P i n a polynomial of degree n in µ. Now we impose eq. (33) 5 at order n and see that its first, second, fifth and sixth terms are of degree n − 2 in µ so that we have 2 ∂ 2 ∆H n+1 ∂µ∂µ ki λ + 2 ∂ 2 ∆H n+1 ∂µ k ∂µ ij λ j = Q n−2 with Q n−2 a polynomial of degree n − 2 in µ. This relation, thanks to (115) becomes 2λ ∂ 2 ∆H n+1 a ∂µ∂µ ki µ a + 2λ ∂ 2 ∆H n+1 0 ∂µ∂µ ki + 2λ j ∂ 2 P n−2 ∂µ k ∂µ ij + 2λ j ∂ ∂∂µ ij ∆H n+1 k = Z n−2 with Z n−2 a polynomial of degree n − 2 in µ. This relation, calculated in µ j = 0, thanks to (118) and to the fact that P n−2 is at least quadratic in µ j , gives But a function depending only on µ and λ cannot be of order n + 1 with respect to equilibrium; it follows that f (µ, λ) = 0. Consequently, (120), (118) and (115) give that ∆H n+1 is a polynomial of degree n in µ and this completes the proof. We note also that the proof of PROPERTY 1 has not required the conditions (33) 3,4 so that it holds also without the symmetry conditions (7) 4,5 ; this will be useful for a future work when we will study this model without these two last symmetry conditions.