Temperature Distribution through a Nanofilm by Means of a Ballistic-Diffusive Approach

As microelectronic devices are important in many applications, their heat management needs to be improved, in order to prolong their lifetime, and to reduce the risk of damage. In nanomaterials, heat transport shows different behaviors than what can be observed at macroscopic sizes. Studying heat transport through nanofilms is a necessary tool for nanodevice thermal management. This work proposes a thermodynamic model incorporating both ballistic, introduced by non-local effects, and diffusive phonon transport. Extended thermodynamics principles are used in order to develop a constitutive equation for the ballistic behavior of heat conduction at small-length scales. Being an irreversible process, the present two-temperature model contains a one-way transition of ballistic to diffusive phonons as time proceeds. The model is compared to the classical Fourier and Cattaneo laws. These laws were not able to present the non-locality that our model shows, which is present in cases when the length scale of the material is of the same order of magnitude or smaller than the phonon mean free path, i.e., when the Knudsen number Kn ≤ O(1). Moreover, for small Kn numbers, our model predicted behaviors close to that of the classical laws, with a weak temperature jump at both sides of the nanofilm. However, as Kn increases, the behavior changes completely, the ballistic component becomes more important, and the temperature jump at both sides of the nanofilms becomes more pronounced. For comparison, a model using Fourier’s and Cattaneo’s laws with an effective thermal conductivity has shown, with reasonable qualitative comparison for small Knudsen numbers and large times.


Introduction
An increasing attention towards nanotechnology has yielded new developments concerning transport phenomena in nanomaterials. It is well known that heat transfer at micro-and nanoscales shows different behaviors from that at macroscales [1,2]. Thermal management is important, due to the miniaturization of electronic devices, which necessitates more attention, to prevent overheating [3]. Local hotspots can reduce remarkably the lifetimes of such electronic devices, and even cause damage to those devices [4,5]. Other studies that investigate thermal management of devices at the nanoscale include the cooling of microfluidic electronic systems [6], enhancing the thermal conductance of 2D materials [7], and improving the efficiency of thermoelectric devices [8], to mention a few. Modelling can provide useful information [9], in order to attain the aforementioned objectives, whether they be sought after or refrained from. Quite some modelling work has been performed in analyzing thermal management in nanofilms or devices at the nanoscale. Works on phonon transport in silicon nanofilms [10,11] make use of numerical methods for predicting thermal conductance.
More specific numerical methods are also employed, such as non-equilibrium molecular dynamics [12], the often-used phonon Boltzmann equation [13], and also phonon hydrodynamics [14]. These studies do not make the distinction between the different types of phonon carriers. Indeed, as the characteristic size of the system becomes of the order of magnitude of the phonon mean free path, ballistic features begin to play a role. To account for the ballistic effects, some studies proposed ballistic-diffusive phonon transport via Monte Carlo simulations [15,16]. Other works regarded the ballistic-type of phonon transport as an introduction of a non-local effect, obtaining via the phonon Boltzmann equation the Guyer-Krumhansl equation [17][18][19]. Whilst these works appeared to provide quite correct predictions (and as will be mentioned later, our model will give comparable results), these works either demand a high computing cost [15,16], or they lack the distinction between the diffusive and ballistic components [17][18][19]. Although they do not make an explicit distinction between the ballistic and diffusive components, this does not necessarily mean that the model is incorrect (it can indeed predict correctly); such a distinction can provide more information on the behavior of the phonons. The model that we propose in this work contains both non-local effects, and makes distinctions between the diffusive and ballistic components. Furthermore, it will be developed from thermodynamic principles (and not from numerical nor kinetic bases), which makes it easily connectable to other physical phenomena, such as mass transport, electric current, and so on [1]. In a greater framework, this work would fit easily without great adaptations nor high computing efforts.
Whether heat transport is of diffusive of a ballistic nature depends mainly on the size of the system. A typical way to make a rough quantitative distinction between these two heat transport mechanisms is the so-called Knudsen number Kn ≡ /L, where denotes the mean free path of the carriers in question (phonons in case of heat transfer) and L a characteristic dimension of the system under study. For large Knudsen numbers (typically larger than one), the heat transport is referred to as ballistic, i.e., dominated by phonon collisions with the walls. If the Knudsen number is much smaller than one, heat transport is simply diffusive, i.e., dominated by phonon-phonon collisions inside the system, and described by Fourier's law. For small-length systems (of the order of magnitude of the mean free path of the carrier or smaller), as well as for high-frequency processes (of which the inverse is of the order of magnitude of the relaxation time of the carrier or smaller), Fourier's law is no longer valid. These observations justify the need to propose a constitutive law that generalizes Fourier's law, emanating from extended well-known established thermodynamic laws. Although there are various ways by which this can be achieved [20][21][22], we prefer to present a solution that is easily connectable with existing balance equations, and recognizable with more simplified constitutive equations. The classical Fourier law: relating the heat flux vector q to the temperature gradient ∇T, with λ denoting the thermal conductivity, is not applicable at short time periods and small spatial scales. In order to account for high frequencies, Fourier's law has been generalized by Cattaneo [23] under the form: with τ designating the relaxation time of the heat flux and ∂ t the partial time derivative. Nonetheless, Cattaneo's relation cannot be used for nanosystems. For this, we propose an extended, generalized constitutional equation for heat conduction, which is valid for nanosystems, systems with high frequencies, and systems with large relaxation times. For this, we use Extended Non-Equilibrium Thermodynamics [1,24], proposing a procedure that can be easily generalized to other transport phenomena, such as electrical conduction and matter diffusion, but we limit ourselves here to heat transport. Applied to heat transport, the theory has as main characteristic to upgrade the heat flux, and higher-order heat fluxes to the rank of independent variables as the internal energy (temperature). The obtained result will then be of use for describing ballistic heat transport, when we propose a two-temperature model, in order to describe heat transport through a thin slab. At the

Methods: Constitutional Equation for Ballistic Heat Transport
For rigid material, we can define an entropy function

Methods: Constitutional Equation for Ballistic Heat Transport
For rigid material, we can define an entropy function . Let us assume that this entropy ( , ) is depending on both the internal energy and the heat flux vector . The corresponding Gibbs equation in a rigid conductor (i.e., ≡ ) at rest writes as: wherein and are measured per unit volume, 1 is a phenomenological coefficient assumed to be independent of , and a dot stands for the scalar product. The coefficient 1 will be identified later on, where it will be shown to be related to the relaxation time and the thermal conductivity . It should be noted that (3) does not take into account any terms that allow for coping with nonlocal effects, which are dominant at small length scales. To take them into account, it is suggested to introduce a hierarchy of fluxes (1) , (2) ,..., ( ) with (1) being identical to the heat flux vector , where (2) (a tensor of rank two) is the flux of , (3) the flux of (2) , and so on. Up to the ℎorder flux, the generalized Gibbs equation can then be written as: wherein the symbol ⨂ denotes the inner product of the corresponding tensors, and = 2,3, … . Moreover, the time evolution of the entropy is governed by a general balance equation, which can be written as: wherein stands for the entropy flux and for the rate of entropy production per unit volume, which is positive definite to fulfil the second law of thermodynamics. In order to deduce an expression for , we use (4) for , and we need a constitutive relation for in terms of the set of variables established earlier. It is natural to expect that is not simply given by the classical expression −1 [25], but that it will depend on all the higher order fluxes up to order , namely: with designating the phenomenological coefficients. We substitute in (5) the expressions of and by (4) and (6), respectively. We then replace derivative ∂ t u by the energy conservation law, = −∇ • , which leads to the expression for the entropy production: The above is a bilinear expression in fluxes and forces (the quantities between parentheses), which suggests the following hierarchy of linear flux-force relations: It is the purpose to replace the set of relations (8)-(9) by one single constitutive equation that takes into account all the ℎ order fluxes. For the sake of clarity, we will make the development up to the fourth order, and generalize afterwards. The fourth order flux equation ( = 4 in Equation (9)) is given by: Let us now take the divergence (∇) of (10), substituting it subsequently in (9), with = 3, and omitting any flux, > 4. This gives: . Let us assume that this entropy Inventions 2019, 4, x a comparison is made with Fourier's and Cattaneo's laws, adapted with a developed ef conductivity.

Methods: Constitutional Equation for Ballistic Heat Transport
For rigid material, we can define an entropy function . Let us assume that this e is depending on both the internal energy and the heat flux vector . The corres equation in a rigid conductor (i.e., ≡ ) at rest writes as: wherein and are measured per unit volume, 1 is a phenomenological coefficie be independent of , and a dot stands for the scalar product. The coefficient 1 wi later on, where it will be shown to be related to the relaxation time and the therm . It should be noted that (3) does not take into account any terms that allow for co local effects, which are dominant at small length scales. To take them into account, it introduce a hierarchy of fluxes (1) , (2) ,..., ( ) with (1) being identical to the hea where (2) (a tensor of rank two) is the flux of , (3) the flux of (2) , and so on. order flux, the generalized Gibbs equation can then be written as: wherein the symbol ⨂ denotes the inner product of the corresponding tensors, an Moreover, the time evolution of the entropy is governed by a general balance equation written as: wherein stands for the entropy flux and for the rate of entropy production pe which is positive definite to fulfil the second law of thermodynamics. In order expression for , we use (4) for , and we need a constitutive relation for in ter variables established earlier. It is natural to expect that is not simply given b expression −1 [25], but that it will depend on all the higher order fluxes up to orde with designating the phenomenological coefficients. We substitute in (5) the exp and by (4) and (6), respectively. We then replace derivative ∂ t u by the energy con = −∇ • , which leads to the expression for the entropy production: The above is a bilinear expression in fluxes and forces (the quantities between which suggests the following hierarchy of linear flux-force relations: It is the purpose to replace the set of relations (8)-(9) by one single constitutive takes into account all the ℎ order fluxes. For the sake of clarity, we will make the d to the fourth order, and generalize afterwards. The fourth order flux equation ( = 4 i is given by: Let us now take the divergence (∇) of (10), substituting it subsequently in (9), w omitting any flux, > 4. This gives: (u, q) is depending on both the internal energy u and the heat flux vector q. The corresponding Gibbs equation in a rigid conductor (i.e., d t ≡ ∂ t ) at rest writes as:

ethods: Constitutional Equation for Ballistic Heat Transport
For rigid material, we can define an entropy function . Let us assume that this entropy ( , ) pending on both the internal energy and the heat flux vector . The corresponding Gibbs tion in a rigid conductor (i.e., ≡ ) at rest writes as: rein and are measured per unit volume, 1 is a phenomenological coefficient assumed to dependent of , and a dot stands for the scalar product. The coefficient 1 will be identified on, where it will be shown to be related to the relaxation time and the thermal conductivity should be noted that (3) does not take into account any terms that allow for coping with noneffects, which are dominant at small length scales. To take them into account, it is suggested to duce a hierarchy of fluxes (1) , (2) ,..., ( ) with (1) being identical to the heat flux vector , re (2) (a tensor of rank two) is the flux of , (3) the flux of (2) , and so on. Up to the ℎ r flux, the generalized Gibbs equation can then be written as: ( , , (1) rein the symbol ⨂ denotes the inner product of the corresponding tensors, and = 2,3, … . eover, the time evolution of the entropy is governed by a general balance equation, which can be ten as: rein stands for the entropy flux and for the rate of entropy production per unit volume, h is positive definite to fulfil the second law of thermodynamics. In order to deduce an ession for , we use (4) for , and we need a constitutive relation for in terms of the set of bles established earlier. It is natural to expect that is not simply given by the classical ession −1 [25], but that it will depend on all the higher order fluxes up to order , namely: designating the phenomenological coefficients. We substitute in (5) the expressions of by (4) and (6), respectively. We then replace derivative ∂ t u by the energy conservation law, −∇ • , which leads to the expression for the entropy production: The above is a bilinear expression in fluxes and forces (the quantities between parentheses), h suggests the following hierarchy of linear flux-force relations: It is the purpose to replace the set of relations (8)-(9) by one single constitutive equation that s into account all the ℎ order fluxes. For the sake of clarity, we will make the development up e fourth order, and generalize afterwards. The fourth order flux equation ( = 4 in Equation (9)) ven by: ≡ ) at rest writes as: per unit volume, 1 is a phenomenological coefficient assumed to stands for the scalar product. The coefficient 1 will be identified to be related to the relaxation time and the thermal conductivity s not take into account any terms that allow for coping with nonat small length scales. To take them into account, it is suggested to (1) , (2) ,..., ( ) with (1) being identical to the heat flux vector , ) is the flux of , (3) the flux of (2) , and so on. Up to the ℎequation can then be written as: the inner product of the corresponding tensors, and = 2,3, … . he entropy is governed by a general balance equation, which can be y flux and for the rate of entropy production per unit volume, lfil the second law of thermodynamics. In order to deduce an , and we need a constitutive relation for in terms of the set of is natural to expect that is not simply given by the classical will depend on all the higher order fluxes up to order , namely: enological coefficients. We substitute in (5) the expressions of y. We then replace derivative ∂ t u by the energy conservation law, xpression for the entropy production: ression in fluxes and forces (the quantities between parentheses), rarchy of linear flux-force relations: the set of relations (8)-(9) by one single constitutive equation that er fluxes. For the sake of clarity, we will make the development up e afterwards. The fourth order flux equation ( = 4 in Equation (9)) nce (∇) of (10), substituting it subsequently in (9), with = 3, and es: and u are measured per unit volume, γ 1 is a phenomenological coefficient assumed to be independent of q, and a dot stands for the scalar product. The coefficient γ 1 will be identified later on, where it will be shown to be related to the relaxation time τ and the thermal conductivity λ. It should be noted that (3) does not take into account any terms that allow for coping with non-local effects, which are dominant at small length scales. To take them into account, it is suggested to introduce a hierarchy of fluxes Q (1) , Q (2) ,..., Q (n) with Q (1) being identical to the heat flux vector q, where Q (2) (a tensor of rank two) is the flux of q, Q (3) the flux of Q (2) , and so on. Up to the N th -order flux, the generalized Gibbs equation can then be written as: unit volume, 1 is a phenomenological coefficient assumed to nds for the scalar product. The coefficient 1 will be identified e related to the relaxation time and the thermal conductivity ot take into account any terms that allow for coping with nonsmall length scales. To take them into account, it is suggested to (2) ,..., ( ) with (1) being identical to the heat flux vector , s the flux of , (3) the flux of (2) , and so on. Up to the ℎuation can then be written as: inner product of the corresponding tensors, and = 2,3, … . entropy is governed by a general balance equation, which can be lux and for the rate of entropy production per unit volume, the second law of thermodynamics. In order to deduce an , and we need a constitutive relation for in terms of the set of atural to expect that is not simply given by the classical l depend on all the higher order fluxes up to order , namely: ological coefficients. We substitute in (5) the expressions of We then replace derivative ∂ t u by the energy conservation law, ression for the entropy production: sion in fluxes and forces (the quantities between parentheses), chy of linear flux-force relations: wherein the symbol denotes the inner product of the corresponding tensors, and n = 2, 3, . . . N. Moreover, the time evolution of the entropy is governed by a general balance equation, which can be written as: a comparison is made with Fourier's and Cattaneo's laws, adapted with a developed effective thermal conductivity.

Methods: Constitutional Equation for Ballistic Heat Transport
For rigid material, we can define an entropy function . Let us assume that this entropy ( , ) is depending on both the internal energy and the heat flux vector . The corresponding Gibbs equation in a rigid conductor (i.e., ≡ ) at rest writes as: wherein and are measured per unit volume, 1 is a phenomenological coefficient assumed to be independent of , and a dot stands for the scalar product. The coefficient 1 will be identified later on, where it will be shown to be related to the relaxation time and the thermal conductivity . It should be noted that (3) does not take into account any terms that allow for coping with nonlocal effects, which are dominant at small length scales. To take them into account, it is suggested to introduce a hierarchy of fluxes (1) , (2) ,..., ( ) with (1) being identical to the heat flux vector , where (2) (a tensor of rank two) is the flux of , (3) the flux of (2) , and so on. Up to the ℎorder flux, the generalized Gibbs equation can then be written as: wherein the symbol ⨂ denotes the inner product of the corresponding tensors, and = 2,3, … . Moreover, the time evolution of the entropy is governed by a general balance equation, which can be written as: wherein stands for the entropy flux and for the rate of entropy production per unit volume, which is positive definite to fulfil the second law of thermodynamics. In order to deduce an expression for , we use (4) for , and we need a constitutive relation for in terms of the set of variables established earlier. It is natural to expect that is not simply given by the classical expression −1 [25], but that it will depend on all the higher order fluxes up to order , namely: with designating the phenomenological coefficients. We substitute in (5) the expressions of and by (4) and (6), respectively. We then replace derivative ∂ t u by the energy conservation law, = −∇ • , which leads to the expression for the entropy production: wherein J s stands for the entropy flux and η s for the rate of entropy production per unit volume, which is positive definite to fulfil the second law of thermodynamics. In order to deduce an expression for η s , we use (4) for ∂ t

of 16
Fourier's and Cattaneo's laws, adapted with a developed effective thermal al Equation for Ballistic Heat Transport can define an entropy function . Let us assume that this entropy ( , ) internal energy and the heat flux vector . The corresponding Gibbs tor (i.e., ≡ ) at rest writes as: asured per unit volume, 1 is a phenomenological coefficient assumed to a dot stands for the scalar product. The coefficient 1 will be identified hown to be related to the relaxation time and the thermal conductivity (3) does not take into account any terms that allow for coping with nonminant at small length scales. To take them into account, it is suggested to uxes (1) , (2) ,..., ( ) with (1) being identical to the heat flux vector , nk two) is the flux of , (3) the flux of (2) , and so on. Up to the ℎ -Gibbs equation can then be written as: enotes the inner product of the corresponding tensors, and = 2,3, … . ion of the entropy is governed by a general balance equation, which can be entropy flux and for the rate of entropy production per unit volume, e to fulfil the second law of thermodynamics. In order to deduce an (4) for , and we need a constitutive relation for in terms of the set of ier. It is natural to expect that is not simply given by the classical that it will depend on all the higher order fluxes up to order , namely: henomenological coefficients. We substitute in (5) the expressions of ectively. We then replace derivative ∂ t u by the energy conservation law, , and we need a constitutive relation for J s in terms of the set of variables established earlier. It is natural to expect that J s is not simply given by the classical expression T −1 q [25], but that it will depend on all the higher order fluxes up to order n, namely: with β n designating the phenomenological coefficients. We substitute in (5)

Methods: Constitutional Equation for Ballistic Heat Transport
For rigid material, we can define an entropy function . Let us assume th is depending on both the internal energy and the heat flux vector . The equation in a rigid conductor (i.e., ≡ ) at rest writes as: wherein and are measured per unit volume, 1 is a phenomenological be independent of , and a dot stands for the scalar product. The coefficien later on, where it will be shown to be related to the relaxation time and the . It should be noted that (3) does not take into account any terms that allow local effects, which are dominant at small length scales. To take them into acc introduce a hierarchy of fluxes (1) , (2) ,..., ( ) with (1) being identical to where (2) (a tensor of rank two) is the flux of , (3) the flux of (2) , and order flux, the generalized Gibbs equation can then be written as: wherein the symbol ⨂ denotes the inner product of the corresponding tens Moreover, the time evolution of the entropy is governed by a general balance e written as: wherein stands for the entropy flux and for the rate of entropy produ which is positive definite to fulfil the second law of thermodynamics. In expression for , we use (4) for , and we need a constitutive relation for and J s by (4) and (6), respectively. We then replace derivative ∂ t u by the energy conservation law, ∂ t u = −∇·q, which leads to the expression for the entropy production: The above is a bilinear expression in fluxes and forces (the quantities between parentheses), which suggests the following hierarchy of linear flux-force relations: It is the purpose to replace the set of relations (8)-(9) by one single constitutive equation that takes into account all the N th order fluxes. For the sake of clarity, we will make the development up to the fourth order, and generalize afterwards. The fourth order flux equation (n = 4 in Equation (9)) is given by: Let us now take the divergence (∇) of (10), substituting it subsequently in (9), with n = 3, and omitting any flux, n > 4. This gives: Applying again operator divergence on (11), substituting it on its turn in (9) but with n = 2 and repeating the same operation until arriving at Equation (9) with n = 1 (the latter of which is equivalent to Equation (8)), the final result is: This relation may be viewed as a generalization of Fourier's and/or Cattaneo's relation, up to the fourth-order flux Q (4) . However, the use of higher order fluxes is not very convenient from a practical point of view, since their physical meaning, being not directly measurable, is not clear. This calls for proposing a more suitable representation, wherein all the high-order fluxes have been expressed solely in the heat flux.
Therefore, the purpose is to express Equation (12) exclusively in terms of the more physical quantities, the temperature T and the classical heat flux q. The higher-order fluxes Q (n) are eliminated by making use of the very definition of the high-order flux, namely: In virtue of (13) and after some rearrangements, Equation (12) can be rewritten as: By extending the exercise to higher-order fluxes with n = N (where theoretically, N → ∞ , but practically often not necessary), we obtain the following general equation governing heat conduction in rigid bodies: Setting γ 1 = β 1 = 0, Equation (8) reduces to Fourier's law, from which it follows that ν 1 = 1 λT 2 . Comparing the unities in (8) indicates that the dimension of the ratio γ 1 /ν 1 is that of time (t), Defining the time unit as the heat flux (q) and relaxation time τ 1 , one has γ 1 = τ 1 ν 1 = τ 1 λT 2 . As the heat flux q is expressed in W/m 2 , the unity of Q (2) is W/(m·s), according to (13): [t] Q (n) , where [L] is a unity in length. As a direct consequence of the property that Q (2) is the flux of q, one has in (8) that For further application, we neglect higher orders of the relaxation times, i.e., τ 2 = τ 3 = . . . = τ n = 0, which means that the phenomenological coefficients β n and γ n , with n ≥ 2, and ν n , and n ≥ 3 (note that although β 2 and γ 2 are related to τ 2 , ν 2 is rather related to τ 1 ), can be omitted. Equation (15) becomes: Inventions 2019, 4, 2

of 16
We can see clearly differences between (16) and the conventional constitutional laws for heat transport (1) and (2), but also some similarities. For systems where its characteristic size is much larger than the mean free path , Equation (16) reduces to Cattaneo's equation. This corresponds also to neglecting any higher order of heat fluxes (reducing the whole deduction in this section to only Equations (3) and (5), with J s = T −1 q). If, furthermore, the changes in the heat flux occur in a time span that is much larger than the relaxation time, we obtain Fourier's law. Other well-known methods that are used, as those mentioned in the Introduction, comprise the phonon-Boltzmann equation, the Monte Carlo simulations, and Molecular Dynamics. The Monte Carlo simulations and Molecular Dynamics methods do not require higher-order variables, but they demand high computing times, and are not flexible for modelling various scales at a time, although efforts have been performed in performing multi-scale molecular dynamics [26]. Solving the phonon-Boltzmann equation also asks for a high computational time. These are issues that are not a problem for the model presented in this work. However, the definition of the ballistic and diffusive temperatures is not well-defined, which limits direct interpretation. This can be circumvented by noticing that the ballistic and diffusive temperatures are actually defined through their respective internal energies. Assuming that the total internal energy is a measure for the real temperature, the evaluation of the total temperature becomes, therefore, correct. Assuming also that the ballistic and diffusive internal energies sum up linearly to obtain the total internal energy, allows for an interpretation of the contribution of the ballistic and diffusive phonon regimes via the internal energies. The present development could, of course be extended by introducing non-linear heat fluxes. However, such an extension, albeit theoretically not an issue, would be numerically cumbersome to perform, also inducing also boundary conditions that may not be well defined. Non-linear contributions of the heat flux (at first and higher orders) could be important when the length scale becomes much smaller than the mean free path of the phonons. Nevertheless, for most materials, a linear heat flux consideration still gives satisfactory results [1,27], making it very reasonable and consistent to use linear heat fluxes for predicting material properties in several applications.

Temperature Distribution in Nanoscaled Materials
We have discussed that, depending on the size, the type of heat transfer can be either ballistic or diffusive. In order to appreciate the contribution of both the ballistic and diffusive parts of the heat flux, we propose a two-temperature model.

Two-Temperature Model
In order to find the temperature distribution, we need to assume the coexistence of two kinds of heat carriers: diffusive phonons which undergo multiple collisions within the core of the system, and ballistic phonons originating at the boundaries and experiencing collisions with the walls. This representation is called the ballistic-diffusion model. As such, the internal energy and the heat flux are decomposed into a ballistic and a diffusive component, which gives u = u b + u d and q = q b + q d , where the indexes b and d refer to ballistic and diffusive contributions, respectively. The space of state variables is constituted by the quantities u d , u b , q d and q b . Let us introduce the diffusive and ballistic temperatures T d and T b , defined respectively by T d = u d /c d and where c d and c b stand for the heat capacity per unit volume, and are positive quantities to guarantee the stability of the equilibrium state. Stating that c d = c b = c, leading to the total temperature being defined by T = u/c, or equivalently, T = T d + T b . Although the quantities T d and T b bear some analogy with the classical definition of the temperature, it should however, be realized that, strictly speaking, these quantities do not represent temperatures in the usual sense, but must be considered as a measure of the internal energies. After having defined the state variables, we specify their behavior in the course of time and space. The evolutions of the internal energies u d and u b , are governed by the classical energy balance laws: while the total internal energy u = u d + u b satisfies the first law of thermodynamics: The quantities r d , r b , and r designate source terms, which may be either positive or negative. In virtue of (19), one has to satisfy r d + r b = r. In the absence of energy sources (r = 0), one simply has that r d = −r b . Based on kinetic theory considerations [25], it is shown that: where the minus sign indicates that ballistic carriers can be converted into diffusive ones, but that the inverse is not possible. Moreover, τ b is the relaxation time of the ballistic energy flux q b . Now, the evolution equation for the fluxes are to be developed. Concerning the diffusive phonons, it is assumed that they satisfy Cattaneo's equation, to cope with their high frequency properties. As a consequence, we are allowed to write: wherein the relaxation time τ d and the thermal conductivity coefficient λ d are positive quantities to meet the requirements of the stability of equilibrium and the positivity of the entropy production, respectively [23,28]. However, expression (21) is not able to describe the ballistic regime, which is mainly influenced by non-local effects, as most of the ballistic carriers cross the system without experiencing collisions, except with the boundaries. As shown before [17,18,28], this situation is satisfactorily described by Equation (16): where it is reminded that b is identified by the mean free path of the phonons, with the subscript indicating that (only in the present section) a distinction is made between ballistic and diffusive contributions for the mean free path to account for the same distinction for the heat flux. Note that the terms involving the space derivatives of the heat flux vector account for the non-local effects, and they are important when the spatial scale of variation of the heat flux is comparable to the mean free path of the heat carriers. Expressions (17), (18), (21), and (22) provide the basic set of eight scalar evolution equations for the eight unknowns u d , u b , q d , and q b . Eliminating the fluxes, using the aforementioned definitions and balance equations, we obtain two equations for the diffusive and ballistic contributions of the phonons in terms of internal energies. Thin films are defined here as films, of which the thickness L may be of the same order of magnitude, or even smaller than the mean free path of the phonons. Heat capacity and thermal conductivity are assumed to be constant, and to take the same values for the diffusive and ballistic phonons (λ = 1 3 cv d,b d,b ), while internal energy sources are absent (r = 0). Initially, the system is at uniform energy u 0 , or using an equivalent terminology, at the "quasi-temperature" T 0 , related to u 0 by u 0 = cT 0 . Often, for thin films in real applications, the other dimensions are much larger, so that we can neglect their effect and consider only a one-dimensional model, with coordinate z. At t = 0, the lower surface z = 0 is suddenly brought to the "quasi-temperature" T l = T 0 + ∆T, while the upper surface z = L is kept at "quasi-temperature" T 0 , equal to the so-defined reference temperature. For further purpose, we introduce the Knudsen numbers Kn i = i /L (i = d, b), which can be given the more general form: It is convenient to use dimensionless quantities for numerical simulations, using the following scalings (the asterisks denoting the dimensionless parameters): with θ d , θ b , and θ (= θ d + θ b ) designating the non-dimensional energies (or temperatures) associated with the diffusive, ballistic, and total energies, respectively. Note that u d (z = L) and u b (z = L) are chosen, such that they are taken as reference values (the sum of which corresponds to the aforementioned u 0 ), which are in the dimensionless form equal to zero, i.e., θ d (z = L) = θ b (z = L) = 0 (as seen later). The new evolution equations for the dimensionless temperature distributions now take the form:

Boundary and Initial Conditions
The purpose is to solve for θ d (z * , t * ) and θ b (z * , t * ). At t = 0, the sample is at uniform temperature T 0 , which implies that the total energy is given by u(z, 0) = u d (z, 0) + u b (z, 0) = cT 0 . However, it is reasonable to suppose that at short times, the ballistic phonons are dominant, so that the initial energy will be essentially of ballistic nature, leading to u b (z, 0) = cT 0 , or in dimensionless notation: Throughout the sample, at time t = 0, the heat flux q is also zero; as a consequence of the energy balance (17), (18), it is checked that initially, ∂θ(z * ,t * ) ∂t * = 0. This result remains, in particular, satisfied under the assumptions: The formulation of the boundary conditions is a more delicate problem. Their importance has to be underlined, because in nanomaterials, their influence is felt throughout the whole system. To satisfy the conditions θ(0, t * ) = 1 and θ(1, t * ) = 0, the simplest tentative would be to suppose that, at z * = 0, θ b (0, t * ) = 1, together with θ d (0, t * ) = 0, while at z * = 1, the temperatures of both the ballistic and diffusive constituents would be zero. However, such expressions are too simple, and do not, in particular, cope with temperature jumps, due to thermal boundary resistance. This is the reason for why we have considered the following boundary conditions for the ballistic carriers: The quantity a, which represents the temperature jump of the ballistic phonons at the face z * = 0 at t * = 0, is taken as being equal to 1 2 . This value may be understood statistically. Since the temperature boundary condition at z * = 0 actually represents an internal energy boundary condition, it can be said that the ballistic phonons that are generated at the heated face are formed by half of the carriers at the initial internal energy θ b = 0, and the other half at the value θ b = 1, corresponding to the energy at the face where the temperature is suddenly increased. A discussion on this result can be found in [25], where it was mentioned that, by solving Boltzmann's equation, θ b (0, t * ) = 1 2 at the heated boundary z * = 0. A posteriori, it is shown later on that this value leads to results, which match satisfactorily well with other different approaches. Concerning the diffusive carriers, we assume that both of the interfaces are black phonon emitters and absorbers, implying that the boundaries are made of incident diffusive carriers only. Combining Cattaneo's equation [23] and Marshak's boundary condition [29,30] for black body thermal radiation, one obtains for z * = 1, 0: where the factor (Kn d /Kn b ) 2 stems from the non-equality of the relaxation times.

Comparison with Fourier and Cattaneo Laws
The solutions are compared with those stemming from the combination of the dimensionless Fourier and Cattaneo laws (1)-(2) with the energy balance (19), with r = 0: with θ f and θ c being the dimensionless Fourier and Cattaneo temperatures. Since no ballistic-diffusion transitions are considered in these latter laws, the boundary conditions (29)-(31) are not used. Instead, the classical (for the example considered in this paper) boundary conditions are used: The same initial conditions as for the two-temperature model are used.

Numerical Method
The model is solved by using a semi-implicit difference method for the spatial discretization and the forward Euler method for the temporal discretization for each time step t * . At the beginning of the simulation, the initial conditions are imposed. At a certain time t * , the equations for the temperatures are written in matrix form, A * V = B, where A is a matrix of dimension (n − 2) × (n − 2) (n being the number of nodes in the discretized spatial domain, and two nodes being dedicated to the boundary conditions; the bulk equations are thus discretized on (n − 2) nodes), containing the coefficients of the temperatures at time t * + 1, V, a vector of dimension (n − 2) containing the unknown temperature values at the spatial nodes at time t * + 1, and B being a vector of dimension (n − 2), containing the known temperature values at time t * . Due to the spatial gradients in the bulk equations, the values at nodes 2 and n − 1 depend on the boundary values at nodes 1 and n, respectively. The vector V is calculated as V = A −1 B. The boundary condition for the ballistic temperature at the hot side is imposed and held at that value at each time step. The boundary condition for the ballistic temperature at the cold side as well as the boundary conditions of the diffusive temperature at both sides at time t * + 1 are then obtained through the bulk values at time t * + 1. This procedure is repeated until the relative errors of the bulk and boundary values of both the ballistic and diffusive temperatures between the previous and present loops become smaller than 10 −4 . The obtained values at time t * + 1 are stocked in the output matrix, and used as the known values for the next time step t * + 2. This procedure is outlined in a schematic form in Figure 1. boundary condition for the ballistic temperature at the cold side as well as the boundary conditions of the diffusive temperature at both sides at time * + 1 are then obtained through the bulk values at time * + 1. This procedure is repeated until the relative errors of the bulk and boundary values of both the ballistic and diffusive temperatures between the previous and present loops become smaller than 10 . The obtained values at time * + 1 are stocked in the output matrix, and used as the known values for the next time step * + 2. This procedure is outlined in a schematic form in Figure 1.

Temperature Distribution of the Two-Temperature Model
The model described in the previous section has been solved using Mathematica. In Figure 2, the non-dimensional temperature profiles are shown for n = 0.1, 1 and 10, where n = n = Kn. These values for n are chosen, since they all represent cases where both the diffusive and ballistic regimes are present, with a dominant diffusive regime for n = 0.1, and a dominant ballistic one for n = 10. The results are presented for three non-dimensional times: * = 1, 10, and 100. In order to emphasize the relative importance of the two contributions in our model, both the ballistic and diffusive components of the temperature are shown next to the non-dimensional temperature (labelled as "Total"). For comparison, the results from Fourier's and Cattaneo's models are presented as well. We notice that generally, the region close to the hot side ( * = 0) is mainly dominated by the ballistic component contribution, which decreases with space, while the diffusive one is increasing

Temperature Distribution of the Two-Temperature Model
The model described in the previous section has been solved using Mathematica. In Figure 2, the non-dimensional temperature profiles are shown for Kn = 0.1, 1 and 10, where Kn d = Kn b = Kn. These values for Kn are chosen, since they all represent cases where both the diffusive and ballistic regimes are present, with a dominant diffusive regime for Kn = 0.1, and a dominant ballistic one for Kn = 10. The results are presented for three non-dimensional times: t * = 1, 10, and 100. In order to emphasize the relative importance of the two contributions in our model, both the ballistic and diffusive components of the temperature are shown next to the non-dimensional temperature (labelled as "Total"). For comparison, the results from Fourier's and Cattaneo's models are presented as well. We notice that generally, the region close to the hot side (z * = 0) is mainly dominated by the ballistic component contribution, which decreases with space, while the diffusive one is increasing up to a maximum (not so discernible for Kn = 10), after which one observes a decrease. This descent is the steepest for smaller Kn. As expected, the influence of the ballistic component becomes more important as Kn is larger, whilst that of the diffusive component is dominant for smaller values of Kn. Moreover, although the ballistic component develops earlier, it appears that afterwards, the diffusion component is growing more over time than the ballistic one. This observation reflects the conversion of the ballistic internal energy into the diffusive one in time, but the reverse is not possible. We can also observe that for small Kn values, the temperature at the cold side is close to zero, the reference value. The fact that this temperature for Kn = 0.1 is not exactly zero, is due to the weak, non-zero, ballistic component that is responsible for well-known temperature jumps. These temperature jumps are more pronounced at larger Kn values, and are clearly seen for Kn = 10. Finally, it can be seen that Fourier's and Cattaneo's laws are qualitatively close to our model for small values of Kn (more dominant diffusive regime), but they differ considerably for larger values of Kn (higher ballistic contribution). up to a maximum (not so discernible for n = 10), after which one observes a decrease. This descent is the steepest for smaller n. As expected, the influence of the ballistic component becomes more important as n is larger, whilst that of the diffusive component is dominant for smaller values of Kn. Moreover, although the ballistic component develops earlier, it appears that afterwards, the diffusion component is growing more over time than the ballistic one. This observation reflects the conversion of the ballistic internal energy into the diffusive one in time, but the reverse is not possible. We can also observe that for small Kn values, the temperature at the cold side is close to zero, the reference value. The fact that this temperature for n = 0.1 is not exactly zero, is due to the weak, non-zero, ballistic component that is responsible for well-known temperature jumps. These temperature jumps are more pronounced at larger n values, and are clearly seen for n = 10.
Finally, it can be seen that Fourier's and Cattaneo's laws are qualitatively close to our model for small values of n (more dominant diffusive regime), but they differ considerably for larger values of Kn (higher ballistic contribution). The respective contributions of the ballistic, diffusive, and total temperatures are shown and compared to the ones obtained from Cattaneo's and Fourier's equations.

Temperature Distribution of the Fourier and Cattaneo Laws with an Effective Thermal Conductivity
Another way of studying heat conduction through nanomaterials is to use Fourier's or Cattaneo's law, but with an effective thermal conductivity, λ , that takes into account the size of the material. We propose to write λ in the form:

Temperature Distribution of the Fourier and Cattaneo Laws with an Effective Thermal Conductivity
Another way of studying heat conduction through nanomaterials is to use Fourier's or Cattaneo's law, but with an effective thermal conductivity, λ eff , that takes into account the size of the material. We propose to write λ eff in the form: where λ stands for the bulk thermal conductivity, and the contributions linked to the size-effects of the nanoparticles are described by the correction factor, f c . To determine f c , it is the idea to upgrade the heat flux and its higher order fluxes to the rank of independent variables at the same footing as the energy or the temperature. This procedure has been presented in Section 2. The strength of this development is that from Equations (8) and (9); not only constitutive equations can be obtained, but they can also be used to propose effective properties. In the case of heat transfer, this concerns effective thermal conductivity. So, continuing from Equations (8) and (9), we present how this can be obtained. The physical meanings of the phenomenological coefficients are also already discussed and obtained in Section 2. We assume now that the system is described by an infinite number of flux variables, i.e., N → ∞ . Applying the spatial Fourier transform q(k, t) = +∞ −∞ q(r, t)e −ik Inventions 2019, 4, x a comparison is made with Fourier's and Cattaneo's laws, adapted with a developed conductivity.

Methods: Constitutional Equation for Ballistic Heat Transport
For rigid material, we can define an entropy function . Let us assume that thi is depending on both the internal energy and the heat flux vector . The corre equation in a rigid conductor (i.e., ≡ ) at rest writes as: wherein and are measured per unit volume, 1 is a phenomenological coeffi be independent of , and a dot stands for the scalar product. The coefficient 1 w later on, where it will be shown to be related to the relaxation time and the ther . It should be noted that (3) does not take into account any terms that allow for c local effects, which are dominant at small length scales. To take them into account, introduce a hierarchy of fluxes (1) , (2) ,..., ( ) with (1) being identical to the h where (2) (a tensor of rank two) is the flux of , (3) the flux of (2) , and so on order flux, the generalized Gibbs equation can then be written as: ( , , (1) wherein the symbol ⨂ denotes the inner product of the corresponding tensors, a Moreover, the time evolution of the entropy is governed by a general balance equati written as: wherein stands for the entropy flux and for the rate of entropy production which is positive definite to fulfil the second law of thermodynamics. In orde expression for , we use (4) for , and we need a constitutive relation for in t variables established earlier. It is natural to expect that is not simply given expression −1 [25], but that it will depend on all the higher order fluxes up to or with designating the phenomenological coefficients. We substitute in (5) the ex and by (4) and (6), respectively. We then replace derivative ∂ t u by the energy c = −∇ • , which leads to the expression for the entropy production: The above is a bilinear expression in fluxes and forces (the quantities betwe which suggests the following hierarchy of linear flux-force relations: It is the purpose to replace the set of relations (8)-(9) by one single constituti takes into account all the ℎ order fluxes. For the sake of clarity, we will make the to the fourth order, and generalize afterwards. The fourth order flux equation ( = 4 r dr to relations (8) and (9), with k as the wave-number vector and r as the position vector, one is led to the following evolution equation of the Fourier-transformed heat flux: where λ(k) is wavelength-dependent thermal conductivity, taking the form of a continued fraction expansion [24,28,31], so that the correction factor (this continued fraction divided by λ) finally becomes: where l n is the correlation length of order n defined by l 2 n = β 2 n ν n ν n+1 = 2 , following the discussion under (15). This is reasonable if one wants to obtain constitutive relations, which was the case in Section 3 (even the first order of fluxes is, in most cases, sufficient). Furthermore, it is quite necessary to limit the order of fluxes, so that analytical solutions can be found, or numerical simulations can be performed. However, when one wants to use Fourier's or Cattaneo's law, it is necessary to take into account an infinite number of higher-order fluxes in the framework of an effective thermal conductivity. In that case, we need to use kinetic theories [32,33]. In case we consider heat transfer in one direction (such as through thin films), kinetic theories [32,33] stipulate that the correlation lengths should be selected as l 2 n = a n+1 2 , with a n = n 2 / 4n 2 − 1 , and is identified as the mean free path, independent of the order of approximation. Generally, kinetic theories [32,33] also hold that the relaxation times τ n (n > 1), corresponding to higher order fluxes, are limited by the first-order one, i.e., τ 1 . In the limit of the latter, we take τ n = τ 1 = τ. In the case of nanoparticles or nanofilms, there is only one dimension, namely the radius a p,s of the spheres, or the thickness of the nanofilms. It is therefore natural to define k = k ≡ 2π/a p,s . With these results in mind, and taking the stationary solution, the continued fraction (34) reduces to an asymptotic limit [32,34], leading finally to the following expression for f c : where the Knudsen number is defined as Kn = /a p,s . Combining (34) and (37) leads to the effective thermal conductivity: With (38) instead of λ, we can define effective Fourier and Cattaneo laws. Using these effective laws in combination with (19) with r = 0, we finally obtain in dimensionless form: Figure 3 compares the total temperature from our model (from Figure 2) to the ones obtained from (39) and (40).   Figure 3 shows that, although no exact correspondence is observed, for small numbers and large times, the effective Cattaneo and Fourier results are qualitatively not far from our results. The temperature jumps are not represented by these effective models, and for large Kn numbers and small times, even the qualitative comparison does not hold. In general, the effective model, presented in this section, is reasonable for cases where the number is not too high ( ≤ (1)) and the system is close to stationarity (at large times, moving slowly, and being close to their quasi-stationary state).   Figure 3 shows that, although no exact correspondence is observed, for small Kn numbers and large times, the effective Cattaneo and Fourier results are qualitatively not far from our results. The temperature jumps are not represented by these effective models, and for large Kn numbers and small times, even the qualitative comparison does not hold. In general, the effective model, presented in this section, is reasonable for cases where the Kn number is not too high (Kn ≤ O(1)) and the system is close to stationarity (at large times, moving slowly, and being close to their quasi-stationary state).

Discussion
The work started by showing principles, using Extended Non-Equilibrium Thermodynamics, in order to obtain constitutive equations that are applicable to systems at nanoscale. Since the subject of this paper is focused on heat transport, a generalized equation for heat transport is used for the purposes of considering the problem of heat conduction through a rigid nanofilm. The central assumption of the present work is that, contrary to previous approaches, the set of variables, namely the internal energy and the energy flux, is split into contributions of diffusive and ballistic natures. Heat transport is represented here as a two-component diffusion-reaction process, where ballistic particles can convert irreversibly into diffusive ones. The latter are obeying a Cattaneo equation, while the behaviors of the ballistic phonons are dictated by the aforementioned generalized Equation (16) and (22). This choice stems from the premise that non-local effects are dominating in ballistic collisions. Often, Fourier and/or Cattaneo laws are used, and are adapted by effective material properties, in order to take into account the mentioned non-local effects. Not expecting exact correspondence, it is interesting how close such an approach can be to our two-temperature model. An effective thermal conductivity is derived from our theory, and introduced into the Fourier and Cattaneo laws, for the purposes of comparison with our model. The results suggest that the comparison is reasonable for Kn ≤ O(1) and large times.
It should be noted that one of the objectives is to demonstrate the versatility of Extended Non-Equilibrium Thermodynamics, proposing relatively simple solutions to delicate and complicated phenomena at the nanoscale level. However, it is noteworthy to mention that questions may be raised about the definitions of diffusive and ballistic temperatures at nanoscales. This is the reason for why in our development, we considered the two temperatures as measures of internal energies. As such, the diffusive and ballistic components are to be viewed as quasi-temperatures that stand for the amount of internal energy that is of either diffusive or ballistic natures. However, the total internal energy is a well-defined property, and it directly representative for the temperature.
Our approach in this work can be extended beyond the linear regime, considering the non-linear contributions of the heat flux. Nonetheless, it should be realized that this formulation forms the necessary framework for further investigations. A set of variables has been defined, but the most appropriate set of variables is not easy to define: should non-linear terms in the heat flux and their higher orders be included as well, and if so, would this contribute significantly to the correct values of the heat flux? Answering this asks for a more complex mathematical formulation that is not too complicated to formulate, but not easy to solve, nor to interpret physically. Finally, the choice of the boundary conditions is still a delicate issue, and calls for more study.
All in all, the present model presents a solid framework of investigations in heat conduction through nanoscopic materials, and it is shown to be applicable to a transient problem of heat conduction through a nanofilm. Other works, based on the theory put forward here, have proven to be validated against other works [25,30,35].

Conclusions
The aim of this work is to show how insight can be obtained into the ballistic and diffusive components of phonon heat transport through a nanofilm, by using extended developments from thermodynamic principles. Our model shows good agreement with other sophisticated models in the literature [25,30,35]. The importance of the ballistic and diffusive components depends on the Kn number. For low Kn numbers, the ballistic component is almost absent, and only persists close to the boundary on the hot side of the nanofilm. This is understandable, since the continuous heat flux passing from the heat source into the nanofilm causes a flux of incident phonons, which are transformed irreversibly almost immediately into diffusion ones. For Kn numbers that are higher, these incident phonons evolve further into the nanofilm, and much less is converted into diffusion phonons. This causes completely different temperature profiles. For small Kn numbers, the result is a temperature profile that is close to the ones predicted by the classical laws. For higher Kn numbers, the temperature profiles are completely different from those from the classical Fourier and Cattaneo laws. One may also observe pronounced temperature jumps on both sides of the nanofilm. For comparison, and also because it is used rather often, the classical laws with effective thermal conductivities are also compared to our complete model. It appears that the comparison is rather satisfactory for small Kn numbers, but not at all for higher Kn numbers. Finally, one may state that for higher Kn numbers, it is necessary to use non-local effects, which are explained by our model to come from the ballistic phonons that are only slightly transformed into diffusive ones.
Funding: This research was funded by BelSPo.