Thermal pulse propagation beyond the Maxwell–Cattaneo theory: Application to one-dimensional nanosystems

A non-local and non-linear thermodynamical model of heat transfer at nanoscale beyond the well-known Maxwell–Cattaneo theory is derived. The compatibility of the proposed model with second law of thermodynamics is proved. The model is subsequently used to investigate the propagation of a heat pulse in one-dimensional nanosystems in the linear case. The predicted results are compared with those arising from the Maxwell–Cattaneo theory in order to point out the possible influence both of the non-local effects, and of the relaxation effects of the higher-order fluxes. Some problems related to initial data and boundary conditions are also discussed.


Introduction
In the general context of continuum mechanics, a very important role is played by the so-called balance laws. From the theoretical point of view, a balance law is just the statement of the causes of the contents time-changes of a certain extensive physical quantity over a given domain. Those causes can be only of the following two different types: (i) Flux. The physical quantity varies owing to possible exchanges with the surroundings through the boundary of the domain. (ii) Production. The physical quantity varies since some (thermodynamic) process may occur inside the domain. a great importance deserves the (local) balance equation of the specific entropy s (per unit volume), namely, with J (s) i being the specific-entropy flux, and σ (s) the specific-entropy production (per unit volume). In Eq. (1) the term J (s) i, i (which accounts for the rate of local-entropy exchanged with the surroundings) may be zero, positive, or negative; according to the second law of thermodynamics, instead, in Eq. (1) the quantity σ (s) (deriving from processes occurring inside the system) has to be always non-negative (namely, it is zero at equilibrium, or for reversible transformations, whereas it is larger than zero for irreversible transformations) in any point and at any time, whatever the thermodynamic process is [1,2]. Since σ (s) points out the sources of irreversibility occurring in any thermodynamic process, it should be (directly, or indirectly) related to all quantities which characterize the irreversible process at hand. In non-equilibrium thermodynamic theories, in fact, the physical constraint σ (s) ≥ 0 ( 2 ) plays a very important role: since all solutions of the balance equations have to comply with the condition in Eq. (2), in those theories that unilateral constraint is mainly used to determine the forms either of the constitutive relations [3], or of the evolution equations of the state-space variables [2]. In pursuing these goals, as it will be briefly shown below in the case of two very well-known theories of heat transfer in a rigid body (just for the sake of simplicity), besides the correct choice of the state-space variables, a crucial point is also the right setting of the relation between the specific-entropy flux J (s) i and the local heat flux q i [4][5][6].

The Fourier theory
In the case of the Fourier (F-) theory of the heat transfer in a rigid body, the specific internal energy (per unit volume) u is the only state-space variable: in the F-theory, therefore, only the evolution equation of u is needed. The time variations of u, indeed, are stated by the balance law of energy, the local formulation of which reads According to the point of view of the balance law, as a consequence, from Eq. (3) one may straightforwardly note that q i is just the flux of u: the heat flux, therefore, is the only cause of time variations of the internal energy, in the absence of a source term (here and in what follows neglected just for the sake of simplicity). Since in this theory s = s (u), moreover, the use of Eq. (3) yields once the usual relation d s is employed, 1 with θ being the non-equilibrium temperature [2,7,8]. If we recall that the classical constitutive assumption for the entropy flux states that the specific-entropy flux only depends on the heat flux and the non-equilibrium temperature [1], i.e., J (s) then Eq. (4) can be also rewritten as A direct comparison between Eqs. (1) and (7) finally allows to recognize that 1 Whenever the internal energy is no longer the only one state-space variable (as it will be through this paper), Eq. (5) properly reads ∂s ∂u = 1 θ .
In order to avoid any possible violation of the unilateral constraint (2), in the F-theory it is therefore sufficient to assume as the constitutive equation for the heat flux (which is worldwide known as the F-law), if it is tacitly understood that the positive-valued constant λ means the thermal conductivity. The coupling of Eqs. (8) and (9), in fact, lets that the form of the specific-entropy production reads from which it is easy matter to see that σ (s) = σ (s) (u) is always a positive-definite function.

The Maxwell-Cattaneo theory
The Maxwell-Cattaneo (MC-) theory has been the first attempt to go beyond the validity's limits of the Ftheory. Those limits are particularly evident when the heat transfer is a very fast and/or steep phenomenon (e.g., in the case of ultrasound propagation, the light scattering in gases, shock waves, etc.) [2,9,10], or whenever the relaxation time (i.e., the temporal lagging between the application of a temperature gradient and the appearance of a heat flux) becomes large enough to require a description of the transient regime [11].
According to the point of view of Extended Irreversible Thermodynamics (EIT) [2,7], in the MC-theory of heat transfer in a rigid body the state space is spanned by u and q i ; from the physical point of view, this is tantamount to suppose that both the independent variations of u and q i can, in principle, influence the consequent variations of the other thermodynamic fields, as for example the specific entropy.
Whereas u is still ruled by the balance law in Eq. (3), in this theory q i can be no longer given by a constitutive relation (as it has been done in the case of the F-theory), but an evolution equation is needed for it. In order to theoretically derive it, we may observe that in the MC-theory s = s (u, q i ) in such a way that in the framework of EIT the generalized Gibbs equation read [2,7] with the (positive) material function τ R being the relaxation time of q i , namely, the response time for the onset of the heat flow after a temperature gradient is suddenly imposed. That quantity is substantially the reciprocal of the frequency of resistive phonon collisions, i.e., the mechanisms of phonon scattering in which the phonons' energy is conserved, but not their momentum. 2 From Eq. (11) we have that the time derivative of s reads if Eq. (3) is employed. The use of the constitutive assumption in Eq. (6) allows then to rewrite Eq. (12) as A direct comparison of Eqs. (1) and (13) yields as the constitutive equation for the specific-entropy production. In order to guarantee the (thermodynamic) compatibility of the MC-theory with the unilateral constraint (2), it is still sufficient that σ (s) = σ (s) (u, q i ) is a positive-definite function. For this reason in Eq. (14) we may set which is also the required evolution equation of q i and it is worldwide known as the MC-equation [12]. When Eq. (15) is introduced into Eq. (14), in fact, the following local form of the specific-entropy production for the heat transfer in a rigid body is recovered (16) in the case of the MC-theory. It seems worth noticing that Eqs. (15) and (16) reduce, respectively, to Eqs. (9) and (10) whenever τ R → 0. From the above considerations, it arises that the MC-theory is compatible with the classical constitutive equation for J (s) i [1], too.

Aims and scopes
The increasing interest in nanoscale heat-transport sciences led modern non-equilibrium thermodynamic theories to depart both from the F-theory, and from the MC-theory. Dealing with the heat transfer at nanoscale is not properly an easy task; from the theoretical point of view, to tackle with it several interesting approaches can be found in the literature, as it can be seen, for example, in the outstanding review by Guo and Wang [13]. All those approaches lie either on macroscopic, or mesoscopic, or microscopic methods. Among those methods, one spreading over a large span of time is the phonon hydrodynamics model, the usual starting point of which is the so-called Guyer-Krumhansl equation [2,14], i.e., wherein is the mean-free path of phonons (i.e., quanta of the vibrational mechanical energy arising from oscillating atoms within a lattice) representing the heat carriers in dielectric solids. Equation (17) represents an interesting step beyond Eq. (15) (the latter recovered from the former whenever → 0, for example) since it is able to capture non-local effects in the heat transport, a topic of much interest at nanoscale wherein even small temperature differences may lead to high temperature gradients. In nanosystems, indeed, also non-linear effects may strongly influence the electronic and optical properties: it could be, therefore, important to deeply examine them by introducing generalized nonlinear and non-local heat-transport equations [15]. Whereas the possibility of dealing with a non-local and non-linear model is clearly an important step toward a detailed and more refined description of the heat transfer at nanoscale, it has to be however noted that the introduction of non-linear terms may especially yield additional difficulties among which here we underline the correct setting-up of initial data (ID) and boundary conditions (BCs), which is of itself already a difficult task in the linear case (i.e., when only non-local effects are taken into account). In any model, in fact, both ID, and BCs should be always assigned in such a way that they have to be physically meaningful, besides correctly depicting the experimental reality, and the consequent initial-boundary value problem (IBVP) is mathematically consistent.
Along with the aforementioned observations, in the present paper: • In Sect. 2, we derive a non-linear and non-local model of heat transfer in rigid nanosystems on macroscopic grounds and in agreement with the second law of thermodynamics. That model encompasses Eq. (17) as a special case. The aim of this section is to show how non-local and non-linear effects may naturally arise from Eq. (1). • In Sect. 3, the above model, in its linear formulation, is employed to analytically solve an IBVP which sketches the propagation of second sound in one-dimensional (1D) nanosystems perturbed by a thermal shock. Therein a comparison with the analytical results arising from the MC-theory is also made. The aim of this section is twofold: 1. to point out the possible influence of the relaxation times of higher-order fluxes, as well as the role played by non-local effects on second-sound propagation; 2. to point out a plausible mathematical way of assigning ID and BCs which are both physically consistent, and depict the exact experimental situation under consideration. • In Sect. 4 final comments and future perspectives are given.

Non-local and non-linear heat conduction in rigid solids
Whereas in a theoretical model (dealing with the heat-transfer phenomenon at nanoscale) the non-local effects may be introduced by the spatial derivatives of the state-space variables, the non-linear effects may be instead substantially understood in the following two different ways: (i) by letting the several material functions depend on the state-space variables [16,17]; (ii) by letting products between some state-space variables (or their spatial derivatives) directly enter the model equations [8,18].
Here we especially aim at analyzing the consequence of the second type of nonlinear effects; thus, throughout this paper we assume constant values for all material functions.
Besides the specific internal energy u and the heat flux q i , in our approach we also assume as independent (state-space) variable the flux of the heat flux Q i j . Referring interested readers to Refs. [2,7,13,19,20] for exhaustive theoretical considerations about the physical interpretation of the second-order tensor Q i j , here we only observe that within the kinetic theory of diluted gases one has is the distribution function (r being the position of the generic particle and m its mass), and C is the relative speed of the particles with respect to the barycentric velocity, with C being its modulus. 3 Since Q i j is a symmetric tensor we may split it into its deviatoric (traceless) part o Q i j and its bulk part Qδ i j , namely, in the present theory we postulate that the state space is The basic idea underlying the above selection is that those fluxes, which are genuinely non-equilibrium quantities, constitute a rather natural way of accessing to a system through its boundary.
In agreement with both EIT [2,7], and the Non-Equilibrium with Internal Variables [21], each of the above state-space variables has to be characterized by its own evolution equation. Recalling that the time rate of u is given by Eq. (3), in this section we therefore derive the evolution equations for the above (state-space) flux variables. To do this, we start to observe that the very general form of the constitutive equation of the specific entropy is owing to Eq. (18). For this reason, in agreement with EIT (see Part III-Subsec. 9.3.2 in Ref. [2]), in our approach we generalize the Gibbs equation (11) as wherein the two (positive-valued) material functions τ 1 and τ 0 mean the relaxation times of o Q i j and Q, respectively, and account for the relaxation effects of the higher-order fluxes [2,7]. In some situations, as for example in diluted gases, the relaxation times of higher-order fluxes are of the same order as that of the heat flux itself, in such a way that when q i is taken as an independent variable, all higher-order fluxes should also be taken into account. Further considerations about τ 1 and τ 0 are postponed to Sect. 2.2.
At the present stage, we leave undetermined the other two quantities η 0 and η 1 entering Eq. (19). They however will be clearly identified in Sect. 2.1. For our scope it is enough to claim that they are positive quantities.
The properties of o Q i j and Q allow to suppose that the ratios τ 1 η 1 and τ 0 η 0 remain constant [2], as it will be shown in Sect. 2.2. As a consequence of this from Eq. (19) we may formally obtain the following time derivative of s: If we go beyond the classical constitutive assumption in Eq. (6) and therefore further suppose according to EIT [2], by direct calculations Eq. (20) can be also rewritten as By comparing Eqs. (1) and (22) we may then point out the following constitutive equation for the specificentropy production: Note that Eq. (23) encompasses Eq. (14) as a special case: in fact, the latter form of σ (s) is obtained by the former form of σ (s) whenever o Q i j → 0 and Q → 0. In order to guarantee that the unilateral constraint (2) is always fulfilled, it is sufficient to make assumptions letting the function σ (s) be positive definite. Since Eq. (23) substantially displays the bilinear form wherein A, B and C are strictly negative coefficients, whereas the following quantities are the so-called thermodynamic forces [1,2], a possible way of fulfilling Eq. (2) in any situation is to assume for example. The assumptions in Eqs. (26) straightforwardly lead to

The model equations
Provided the physical meanings of η 0 and η 1 are given, Eqs. (27) are proper the evolution equations of the (state-space) flux variables we're searching for. In order to make those identifications, indeed, we may start wherein the second-order tensor can be substantially meant as an effective thermal conductivity (which could be related to the compelling problem of flux limiters [22][23][24]). Whenever the non-equilibrium-temperature approximation [2,16] holds, namely, if one has with θ eq being the (constant value of the) local-equilibrium temperature, instead of Eqs. (19) and (21), respectively, the non-linear terms proportional to o Q i j θ , j and Qθ , i vanish and the effective thermal conductivity above reduces to the bulk thermal conductivity. In that case, moreover, it is easy matter to observe that Eqs. (17) and (28) coincide provided the following identifications hold: Equations (30) turn out the expressions of η 1 and η 0 in terms of clear physical quantities. If the assumptions in Eqs. (30) are made, from the local balance of energy (3) and Eqs. (27) we finally have which are our model equations.

Comparison with the 9-moment theory
When the non-equilibrium-temperature approximation holds Eqs. (31) become once the classical thermodynamic relation is used, with c v being the specific heat at constant volume per unit volume (always larger than zero, as experimental evidences point out). We may indeed recognize that the model in Eqs. (32) is closely related to the model which characterizes the 9-moment theory [25], i.e., to the model based upon wherein c means the modulus of the phonons' Debye speed (defined as some mean value of longitudinal and transversal speed of sound), and withh being the Plank's constant, k being the wave vector (the modulus of which is k ), the function f (x, t, k) means the phase density. The relaxation time τ appearing in Eqs. (34) is such that where τ N means the reciprocal of the frequency of normal phonon collisions, i.e., the mechanisms of phonon scattering in which both the phonons' energy, and their momentum is conserved. By comparing Eqs. (32) and (34), indeed, it is possible to set the following parallelisms between the field variables involved in those models, the only difference being the form of the isotropic part of the two symmetric second-order tensors Q i j and N i j , which in the latter case is only related to u, as it follows from Eq. (35b). The parallelism between Eqs. (32) and (34), besides giving a common root (very well defined in the kinetic theory) to our model and the 9-moment theory, allows us to give a deeper physical insight into the relaxation times τ 1 and τ 0 : both of them, in fact, should be related to the relaxation time τ defined by Eq. (36) which is tantamount to claim that the relaxation times of higher-order fluxes in our model are related both to the resistive scatterings, and to the normal scatterings. The above considerations justify the assumption of constant values for the ratios τ 1 η 1 and τ 0 η 0 that we have made in Sect. 2. If Eqs. (30) hold, in fact, the above parallelism allows us to state that both τ 1 η 1 , and τ 0 η 0 are proportional to c 2 , which can be surely supposed to be constant.

Second-sound propagation
In the present section we investigate the second-sound propagation predicted by the theoretical model in Eqs. (31). To this end, for computational needs, we introduce the following dimensionless variables: with L being the characteristic length of our system, and T 0 a suitable reference temperature (e.g., the initial temperature of the system). For our goal, we also make use of Eq. (33) such that Eqs. (32a) can be rewritten in terms of the time variation of non-equilibrium temperature θ . In this way Eq. (31) become wherein we set τ 0 = ατ R , τ 1 = βτ R , and Kn = L stands for the Knudsen number. 4 According to what previously said, therefore, in the present analysis the non-dimensional coefficients α and β account for the relaxation effects of the higher-order fluxes (related both to the resistive scatterings, and to the normal ones), whereas Kn introduces non-local effects.
In closing this part, we note that in deriving Eqs. (39) we estimated the resistive relaxation time of phonons as with V 0 standing for the average phonon speed (measured at the reference temperature-value T 0 ), and we used the well-known Ziman limit [26] for the thermal conductivity, i.e., as it is usually assumed in the hydrodynamic theory of second-sound propagation in insulators. Although the assumptions in Eqs. (40) and (41) may lead to a quantitative inaccuracy, it is clear that the results obtained below are not critically dependent upon those assumptions.

The mathematical formulation
Here we put ourself in the 1D case for the sake of simplicity. From the physical point of view this means that we are only interested in the thermal transport within 1D (indefinitely long) nanostructures (e.g., nanowires, nanofibers, nanotubes, and nanorods) which represent one of the best-defined and controlled classes of nanoscale building blocks. We refer to Refs. [27,28] and the references cited therein for some comprehensive reviews about the use of these material systems. Dealing with Eqs. (39) is not properly an easy task, at least if one is searching for analytical solutions which can depict the experimental evidences. One of the main problems in dealing with Eqs. (39) lies, for example, on the question regarding ID and BCs and the way they have to be correctly assigned, as it has been previously observed in Sect. 1. For this reason, as a first attempt to a more general analysis, to reduce our computations to a simpler level, here we put ourself in the non-equilibrium-temperature approximation and therefore we neglect the terms proportional to o Hi j T , j and HT , j in Eq. (39b). By indicating with x the sole (non-dimensional) Cartesian coordinate, as a consequence of the assumptions above, Eqs. (39) reduce to The above system of coupled PDEs will be solved by prescribing the following ID which, from the physical point of view, mean that our rigid conductor is initially at equilibrium and at the reference temperature T 0 . Moreover, we also assign the following BCs wherein H is the Heaviside unit step function. This is tantamount to mathematically depicting the practical situation in which at the time t = 0 + the domain boundary at x = 0 is only affected by a step-jump in temperature. The neglect of non-linear terms has allowed us to not assign any BCs for the other unknown Hxx (x, t) and H (x, t): we let, in this way, the model furnish the values of those functions also at the boundary. Otherwise, we would have had to assign those BCs, thus incurring the problem of the physical meaning of the prescribed values for the unknown functions at the boundary, a problem which is still open. 5 The IBVP in Eqs. (42)-(44), therefore, is able to describe the behavior of an infinite 1D nanosystem which is perturbed from its initial equilibrium state by a thermal pulse (i.e., a thermal shock) applied in its left-hand side end. From the practical point of view, in Eq. (44) T • ∈ R + stands for a dimensionless constant which captures the energy amount of the thermal shock perturbing the system at hand. Further comments about the ID in Eqs. (43) and the BCs in Eq. (44) will be made in Sect. 3.3.

The form of the basic fields
The approximate analytical solution of the IBVP in Eqs. (42)-(44) is the following: wherein we have set just for the sake of a compact notation. Interested readers can refer to "Appendix" at the end of this paper for a sketch of the solution method used to obtain the form of the basic fields in Eqs. (45); therein considerations about the degree of approximation of Eqs. (45) are also made. Here, instead, we aim at drawing the attention on some consequences arising from Eqs. (45). In Fig. 1, we plot the theoretical behavior of the non-dimensional temperature T arising from Eq. (45a) when α = β = 1, T • = 1 and Kn = 1 (just for the sake of computational needs). That figure clearly displays the propagating temperature wave (so pointing out the hyperbolic behavior of Sys. (39) in the 1D linear case): in fact, only the region before the front wave is perturbed (i.e., therein T = 0) whereas the region ahead the front wave is still in equilibrium (i.e., therein T = 0). Owing to dissipative effects, from Fig. 1 it is also possible to see how the energy amount of the thermal wave progressively shrinks for increasing time.
In the case of Eq. (45a), it seems also interesting to observe that the dissipative effects are influenced both by the relaxation times of Q and o Q i j (accounted by α and β, respectively), and by the phonon mean-free path (accounted by Kn).
In particular, we may infer that the temperature-wave energy amount • decreases as steeper as α and β reduce. This can be seen, for example, in Fig. 2 wherein the temperature T behavior, as a function of the distance x from the thermal-shock source, is investigated for different values of α and β. • decreases as steeper as Kn reduces. This can be seen, for example, in Fig. 3 wherein the temperature T behavior, as a function of the distance x from the thermal-shock source, is investigated for different values of Kn. From the practical point of view, we stress that both ID, and BCs should arise from experimental evidences and have their own physical clear meaning. The determination of the initial conditions is generally of no concern either in steady-state situations, or in transient problems, since it is usually admitted that the system at hand is initially in a thermodynamic equilibrium. This is the main reason why ID have been given by Eqs. (43). The right setting-up of BCs, instead, is in general more delicate and controversial. It has to be also noted that in the case of Eqs. (42) the above problem becomes still more acute since therein the four unknown fields Hxx (x, t) and H (x, t) are coupled. As a consequence, the following questions may naturally arise: is the IBVP in Eqs. (42)-(44) physically consistent? Since in Sect. 3.1 only the value of the unknown function T (x, t) at the boundary has been given, are the BCs in Eq. (44) admissible? Whereas the results contained in Sect. 2 make us sure about the physical validity of the theoretical model in Eqs. (42), to answer the above questions in the present section we analyze the consequences of the unilateral constraint (2). In particular, if we introduce the following further dimensionless quantities from Eqs. (21) and (23), respectively, we have By neglecting the non-linear terms proportional to HT , j and o Hi j T , j in Eq. (48b), in the 1D case we have In Fig. 4, we plot the theoretical behavior of the non-dimensional specific-entropy production (s) arising from Eq. (49b) when α = β = 1, T • = 1 and Kn = 1 (just for the sake of computational needs). As it was expected, Fig. 4 shows that only the (out-of-equilibrium) region before the front wave is characterized by (s) values different from zero (but strictly positive, according to second law of thermodynamics). In the region ahead the front wave, which is instead still at equilibrium, it obviously holds (s) = 0. The behavior of (s) is monotonically decreasing: this is logical since in time the temperature wave progressively loses its energy amount, as observed in Sect. 3.2.
The above considerations are enough to claim the validity of the physical setting of the phenomenon at hand, i.e., Eqs. (43) and (44).

Comparison with the MC-theory
The limits of F-theory, which are well known since decades, have led to several interesting generalizations of it. Those generalizations are obviously not unique: they generally differ in the underlying physical principles and in modeling capabilities. In Ref. [21], for example, two theoretical models are discussed and compared: one based upon Rational Extended Thermodynamics, and the other one based upon Non-Equilibrium Thermodynamics with Internal Variables. In Ref. [29], instead, the Grad-13 equations [30] have been embedded into the framework of GENERIC [31]. Since the comparison between different theories is surely one of the major stimuli both to enhance previous models, and to find new others, by means of the second-sound propagation in this part we compare the results of Sect. 3 with those arising from the MC-theory. The MC-theory is considered here only for the fact that it can be recovered (as special case) from any theory beyond the F-law. We postpone to future works a more complete and refined comparison between the different models proposed in the literature.
Since in the MC-theory the two basic fields in such a way that one would have the following analytical solutions instead of Eqs. (45a) and (45b). A sketch of the solution method used to obtain the form of the basic fields in Eqs. (51) can be found in "Appendix," too. In Fig. 5, we compare the behaviors of the analytical solutions in Eqs. (45a) and (51a) when t = 0.5, α = β = 1, T • = 1, and Kn = 1 just for the sake of computation. It clearly points out that the temperature wave predicted by Eq. (45a) (i.e., the dashed red curve therein) is faster than the temperature wave predicted by Eq. (51a) (i.e., the solid blue curve therein): in the former case, in fact, the front wave (characterized by the jump discontinuity) is placed at a distance larger than that at which the front wave in the latter case occurs. This result was indeed expected since the second-sound speed predicted by Eqs. (31) is larger than that predicted by the MC-theory [2,32], whatever the values of the parameter α, β and Kn are.  In the MC-theory the dimensionless 1D specific-entropy flux and production, respectively, read In Figs. 6 and 7, we compare, instead, the behaviors of the analytical solutions by Eqs. (49b) and (52b) when t = 0.8, α = β = 1, T • = 1, and Kn = 1 just for the sake of computation. From that figures, it follows that both the predicted specific-entropy productions tend to a plateau limit which is attained before the front of the moving temperature wave.

Final comments and future perspectives
The raise of nanotechnology has required (and, indeed, is still requiring) several efforts to better understand the thermal-transport properties of nanodevices, the performance and reliability of which could be much influenced by relaxation (i.e., memory), non-local and non-linear effects [2,7,[33][34][35]. For this reason, in Eqs. (31) of this paper we proposed a non-linear and non-local thermodynamical model of heat transfer at nanoscale in a rigid body beyond the usual classical theories. That theoretical model, which encompasses both the F-theory, and the MC-theory as special cases, has been derived in framework of EIT [2,7] in Sect. 2 and essentially lies on the concept of flux of the heat flux Q i j (which has been indeed decoupled into its deviatoric and bulk parts), the physical root of which arises from the kinetic theory [2].  31) is indeed also strictly related to Eq. (23) which substantially displays the bilinear form (24) in fluxes and forces. Although experimental observations confirm that thermodynamic forces and fluxes are interwoven [1,2], the exact relations among them are indeed not known a priori. In our approach we considered linear force-flux relations and, therefore, we finally set Eqs. (25). Non-linear fluxforce relations could have been also considered in principle [36][37][38][39], indeed, in such a way that a different non-linear and non-local model of heat transfer was obtained. In the present paper, however, we did not purse the approach of non-linear flux-force relations since it is generally related to theoretical questions concerning the interpretation of the second law (see Part I-Sec. 2.4 in Ref. [2], for example) which are surely conceptually interesting, but not properly in the focus of this paper.
By means of the so-called non-equilibrium-temperature approximation (i.e., whenever the difference θ −1 − θ −1 eq → 0), in Sect. 3 the non-linear terms in Eqs. (31) have been neglected and the propagation of thermal pulses in 1D nanosystems has been consequently analyzed. That analysis, although simplified, could be however useful in practical applications. For example, in current days the very high circuit densities yield modern chips run so hot that heat is becoming one of the major problems in the design of integrated circuits. As a consequence, one of the aspects which is much being studied at nanoscale is the isolation of very small areas which could be heated by an external source (either static as a hot spot due to a very miniaturized working device, or dynamic as a fast laser pulse): a continuous control of the temperature's rise, in fact, could be a crucial point for the correct operation of modern devices. In Sect. 3 we have tried to exactly tackle with that problem in 1D nanosystems. The analysis performed therein points out that one should pay attention to the role played by relaxation and non-local effects: as it follows from Fig. 2 (for the relaxation effects) and Fig. 3 (for the non-local effects) they may strongly influence the dissipative effects (i.e., one of the main-natural-tool for heat isolation of different components on a chip), thus not ensuring the right setting-up of modern devices which currently remains a very compelling challenge. That challenge could be also complicated by non-linear effects which may either enhance, or reduce, the influence of memory and non-local effects [15].
Another step toward a more detailed investigation will be surely the analysis of a pulse propagation once in Eqs. (31) no terms are neglected, or when non-linear flux-force relations are considered. We note that if Eqs. (31) are considered in their complete version, the interesting problem of flux limiters [22][23][24] also comes out: in that case, in fact, the form of the unknown fields will have to however guarantee that the effective thermal conductivity in Eq. (29) is positive defined.
The analysis performed in Sect. 3 could be, indeed, also used to indirectly compute a first rough estimation of the relaxation times of the high-order fluxes: whereas different models have been developed in order to estimate τ R , very scant information about τ 1 and τ 0 (i.e., information on τ N according with the observations in Sect. 2.2) can be found in the literature. Since the wave front depends on the ratios between those relaxation times and τ R , as Fig. 2 clearly shows, by checking the wave front in practical applications one could also estimate whether τ 1 and τ 0 are larger (or smaller) than τ R , in principle, for example.
At the very end, it seems worth noticing that, in general, it is not so simple to claim which is the best (or, alternatively said, the right) proposal among the different theories beyond the F-law. To this end, practical applications in which different theories clearly predict different results could be suitable benchmarks. This is, for example, just the case of the heat-wave propagation we analyzed in the present paper: in fact, the results contained in Sect. 3 can be used to check, between the MC-theory and the model given by Eqs. (31), which one better fits with experimental evidences.
At the end of this paper, just for the sake of completeness, we give some details about the technique used to solve the above IBVPs, the corresponding results of which are contained in Sect. 3.2.
The inverse (temporal) Laplace transforms of Eqs.

Numerical solution
Although we refer to Eqs. (45) as the analytical solution of the IBVP in Eqs. (42)-(44), it seems worth noticing that they are not exact solutions since we neglected infinitesimals of order lower than ν −1 in deriving Eq. (60). The small-time approximation is surely logical in very fast phenomena (as the pulse propagation at hand), but a priori the orders of infinitesimals which may be neglected are not known.
To be sure that the approximation we used does not significantly affect Eqs. (45), we also perform a comparison with a numerical solution. In Fig. 8, we show how effective is the approximation used to obtain the analytical expression for the temperature field in Eq. (45a) (which is fundamentally the real benchmark in the present approach). In particular therein the solid (blue) curve refers to Eq. (45a), whereas the dotted (red) curve has been obtained by inverting numerically the exact Laplace transform of the temperature field through the modified Tzou series [40] F (x, t) ≈