Mesoscopic Hydro-Thermodynamics of Phonons

A generalized Hydrodynamics, referred to as Mesoscopic Hydro-Thermodynamics, of phonons in semiconductors is presented. It involves the descriptions of the motion of the quasi-particle density and of the energy density. The hydrodynamic equations, which couple both types of movement via thermo-elastic processes, are derived starting with a generalized Peierls-Boltzmann kinetic equation obtained in the framework of a Non-Equilibrium Statistical Ensemble Formalism, providing such a Mesoscopic Hydro-Thermodynamics. The case of a contraction in first order relevant variables is worked out in detail. The associated Maxwell times are derived and discussed. The densities of quasi-particles and of energy are found to satisfy coupled Maxwell-Cattaneo-like hyperbolic equations. The analysis of thermo-elastic effects is done and applied to investigate thermal distortion in silicon mirrors under incidence of high intensity X-ray pulses in free electron laser (FEL) facilities. The derivation of a generalized Guyer-Krumhansl equation governing the flux of heat and the associated conductivity coefficient is also presented.


I. INTRODUCTION
It has been noticed [1] that the ceaseless innovation in semiconductors design creates a demand for better understanding of the physical processes involved in the functioning of modern electronic devices, which operate under far-from-equilibrium conditions. The main interest is centered on the behavior of "hot" carriers (electrons and holes) but the case of "hot" phonons is also of interest, particularly in questions such as refrigeration of microprocessors and heat transport in small systems with constrained geometries [[2]- [4]].
These questions belong to the area of non-equilibrium phonon-dynamics [5], more precisely to the subject of phonon hydrodynamics associated to non-equilibrium (irreversible) thermodynamics [ [6], [7]].
It may be noticed the fact that kinetic and hydrodynamics of fluids are intimately coupled. A first kinetic-hydrodynamic approach can be considered to be the so-called classical (or Onsagerian) hydrodynamics, which gives microscopic (mechano-statistical) foundations to, for example, the classical Fourier's and Fick's diffusion laws. But it works under quite restrictive conditions, (see for example [8]) and, hence, more advanced approaches are required to lift these restrictions, mainly, because, as noticed, the requirements for the analysis of situations involved in the nowadays advanced technologies. Several improved approaches were introduced as, for example, the Burnett approximation of hydrodynamics [9], the one associated to Extended Irreversible Thermodynamics [6], and, recently, the so called Mesoscopic Hydro-Thermodynamics (MHT for short) dealing with heat transport [10] covering, in principle, all kinds of motion, that is, involving intermediate to short wavelengths and ultrafast motion. A complete MHT for classical fluids is reported in Ref. [11].
In this communication we describe the construction of a MHT of phonons, presenting an in depth study of phonon hydro-thermodynamics in the framework of a nonlinear quantum kinetic theory based on a non-equilibrium statistical ensemble formalism [ [12]- [17]], referred to as NESEF for short.
In section II derivation of a phonon higher-order hydrodynamics is presented. The In section III we fully derive a description said of order one, where we retain only the two densities and the first flux of each density.
In section IV we present an analysis of some characteristics of these equations, namely, (a) the Maxwell times [18]  In section V we consider the simpler case when thermo-elastic effects are neglected and the equations for the two densities decouple.
In section VI, on the contrary, thermo-elastic effects are analyzed and used to investigate transient thermal distortion of an optical substrate illuminated by a single intense ultra-short X-ray pulse as currently available in FEL facilities.
Section VII contains our concluding remarks. Details of the calculation are given in several Appendices.

II. PHONON MESOSCOPIC HYDRO-THERMODYNAMICS
We consider a system of acoustic phonons in a semiconductor, in anharmonic interaction among them, in interaction with what we call a thermal bath of other degrees of freedom of the system, and in the presence of an external source capable of driving them out of thermal equilibrium. Moreover, the system is in contact with an external thermostat at temperature T 0 . Since the divergence of a transverse field vanishes, and since in the absence of vortices the rotational also vanishes, in what regards the hydrodynamic motion we consider the longitudinal acoustic (LA) phonons only.

The system Hamiltonian quantum mechanical operator is
where H OS is the Hamiltonian operator of the free LA system, namely, H OS = q ℏω q a + q a q + 1/2 (2) with ω q being the frequency dispersion relation, the wavevector q runs over the 1st Brillouin zone and a q , a + q are annihilation and creation operators for LA phonons in mode q. The H SB term accounts for the interaction of the LA phonons with what we have called the thermal bath, namely, the anharmonic interaction with the TA phonons, the effect of impurities elastic deformation and interaction with the electrons, or holes, they being, bonded or itinerant, the effect of imperfections (mainly dislocations), and we write for it where R q and R + q indicate the corresponding operators associated with transition processes involving the interactions mentioned above; Λ q is the strength coupling in each case. For illustration let us consider the anharmonic interaction operator between the LA and TA phonons, namely where the first contribution on the RHS describes the processes LA(q) → T A(q ′ ) + T A(q ′′ ), the second the processes LA(q) → LA(q ′ ) + T A(q ′′ ), the third LA(q) + LA(q ′ ) → T A(q ′′ ) and the Hermitian conjugate describe the inverse processes. M qk are the appropriate matrix elements. The operators a q b + q−k b + −k lead to the appearance, in the kinetic equation of evolution, of a linear term characterized by a relaxation time τ an q , while the operators a q a + q−k b + −k and a q a q−k b + k produce non-linear terms sometimes referred to as Lifshits [20] and Fröhlich [21] contributions, respectively. It is worth mentioning that these non-linear terms produce, when the system is driven (by the external pumping source) sufficiently far from equilibrium, a so called "complex behavior" [ [22], [23]]. In the case of phonons, it consists in the emergence of a type of Bose-Einstein condensation, propagation of long-lived solitons and a kind of "laser" action [ [24]- [32]]. In what follows we shall take the case of not too intense excitation and disregard non-linear contributions, restricting the analysis to the linear (Onsagerian) regime, when the emergence of "complex behavior" is restrained according to Prigogine's theorem of minimum entropy production in local equilibrium [33].
Moreover, H OB is the Hamiltonian of the subsystems involved in the thermal bath, and H SP is the interaction of the LA-phonons with the external pumping source, to be specified in each practical application.
Next, to apply NESEF it is required first of all to specify the basic dynamic variables used to characterize the non-equilibrium ensemble [ [14]- [17]]. A priori, when the system is driven away from equilibrium, it is necessary to include all observables of the system, which leads to the introduction of many-particle dynamic operators [ [34], [35]]. Here we use the single-particle dynamic operator in the second quantization representation in reciprocal space, with q and Q running over the Brillouin zone.
The two-phonon and higher-order dynamic operators can be ignored because of Bogoliubov's principle of correlation weakening [ [36], [37]]. Moreover, since phonons are bosons, it would be necessary also to include the creation and annihilation operators a q and a + q because their eigenstates are the coherent states [38], and also the pair operators a k a k , ( a + k a + k , ) because the number of quasi-particles is not fixed [39]. However, we disregard them because they are of no practical relevance for the problem at hands. In Appendix A we describe the corresponding non-equilibrium statistical operator. The energy of the thermal bath H OB is also a basic microdynamical variable and then we have the basic set where, Q = 0 corresponds to the occupation number operator ν q = a + q a q describing a homogeneous phonon population and those with Q = 0, namely ν qQ = a + q+Q/2 a q−Q/2 , account for changes in space of the non-equilibrium phonon distribution function.
The average, over the non-equilibrium ensemble, of the microdynamical variables in the set of Eq.(6) provides the variables which characterize the non-equilibrium macroscopic state of the system. Let us call them where Q = 0 and < ν q |t >= T r{ ν q ̺ ε (t)}, etc..., that is, the average of the microdynamical operators of the set (6) where 1/β 0 = k B T 0 .
Going over to direct space we introduce the space and crystal momentum dependent distribution function V is the sample volume which shall be taken as one in what follows, where it may be noticed that ν q (t), Q = 0, refers to the global properties and those with Q = 0 reflect the changes in direct space. Hence, the hydrodynamic variables that describe the hydrodynamic movement consist of two families, namely with ℓ = 2, 3, . . . , describing the movement of the quasi-particles, which we call the n-family, and where is the density of the quasi-particles, is the first (vectorial) flux of the density (current), is the rank-ℓ, with ℓ > 1, (ℓ = 2, 3, 4 . . .) tensorial flux of the density, where meaning tensorial product of ℓ-times the vector ∇ q ω q , rendering a tensor of ℓ-rank.
On the other hand, we have again with ℓ = 2, 3, . . ., describing the movement of energy of quasi-particles, which we call the h-family, and where h(r, t) = q ℏω q ν q (r, t), which are, respectively, the density of energy, its first flux, and the higher order fluxes.
By deriving on time t, the two sides of the general equations (11) to (13) and (16) to (18), there follow, written in compact form, the hydrodynamic equations of motion where p is n or h (for the corresponding families), and and, Evidently, all the hydrodynamic equations of motion are dependent on the equation of motion for the single quantity ν q (r, t). Returning, for practical convenience, to reciprocal Q space, ν qQ (t) satisfies the evolution equation which is the average over the non-equilibrium ensemble (described by the statistical operator with where ν 0 q is the equilibrium LA phonon distribution at temperature T 0 , I ext. q describes the rate of change due to external sources/sinks acting on the system, Γ q is a reciprocal lifetime given in Eq.(25), Π q is the self-energy correction of the LA phonon frequency giving ω q of Eq. (27), which, if in Eq.(3) we consider only the linear anharmonic interaction, that is, the first term on the right of Eq.(4) and its Hermitian conjugate, are given by where Ω k is the k-dependent frequency of TA phonons whose distribution is ν TA k , and P.V. stands for principal value.
Inserting Eq.(23) in Eq. (19), the hydrodynamic evolution equations for the n and h families become Next, a closure for the set of Eq.(29) must be introduced. This means, that we must express the ν q (r, t) appearing on the RHS in terms of the hydrodynamical variables belonging to the set. First thing to notice is that we are dealing with an enormous set of coupled integro-differential equations linking densities and fluxes of all orders. In reference [40] (equation 7) an equivalent representation, in direct r-space, of these equations is given where L p´( r´, t) stands for the contraction in all ℓ´tensor indices. Equation (30) is a continuity equation, having on the RHS a collision integral which accounts for sinks (system relaxation effects) and sources (pumping, driving the system out of equilibrium). It can be considered an extended Mori-Heisenberg-Langevin equation [41].
It represents a quite cumbersome set of coupled equations, of unmanageable proportion. To proceed further it is then necessary to introduce a contraction of description [42], that is to say, a reduction in the number of basic quantities, retaining only a few fluxes. Hence, we must look, in each case, on how to find the best description using the smallest possible number of variables. In other words to introduce an appropriate contraction of description: This contraction implies in retaining the information considered as relevant for the problem in hands, and to disregard nonrelevant information [43].
Elsewhere [42] we have considered the question of the contraction of description (reduction In other words, since MHT implies in describing the motion when governed by smaller and smaller wavelengths, or larger and larger wavenumbers, accompanied by higher and higher frequencies, in a qualitative manner we can say that, as a general "thumb rule," the criterion indicates that a more and more restricted contraction can be used when larger and larger are the prevalent wavelengths in the motion (changes smoother and smoother in space and time). Therefore, in simpler words, when the motion becomes more and more smooth in space and time, the more reduced can be the dimension of the basic macrovariables space to be used for the description of the nonequilibrium thermodynamic state of the system.
It can be conjectured a general criterion for performing contractions, namely, a contraction of order r (meaning keeping the densities and their fluxes up to order r ) can be introduced, once we can show that in the spectrum of wavelengths, which characterize the motion, predominate those larger than a "frontier" one, λ 2 (r,r+1) = v 2 θ r θ r+1 , where v is of the order of the thermal velocity and θ r and θ r+1 the corresponding Maxwell times [42].
In the next section we consider the situation when it is possible to use a contracted description including only the densities and their first fluxes.

III. ANALYSIS OF THE PHONON MHT OF ORDER 1
Let us consider a contracted description including the densities of quasi-particles n(r, t) and energy h(r, t) plus their first fluxes I n (r, t), I h (r, t) only. For practical convenience it is better to work in reciprocal Q space. These quantities are then given by It is recalled that where (see Appendix A) with, in this contracted description, Using the generalized Peierls-Boltzmann Eq.(B6), the resulting equations of motion Next, we need to proceed to the closure of the equations, that is to express both ν qQ (t) and the 2 nd order fluxes in terms of the four variables n, h, I n and I h .
Let us now analyze the contents of these equations. In Eq.(53) (for the density n(r, t) of quasi-particles) and Eq.(54) (for the density of energy h(r, t)), the first term on the RHS is the one of conservation, that is, the divergence of the corresponding flux; however, they are modified by the presence of contributions a [2] 13 (t) and a [2] 24 (t) arising from the self-energy correction. The third term on the RHS corresponds to a relaxation-type contribution; the corresponding coefficients b 11 (t) and b 22 (t) are minus the inverse of Maxwell characteristic times which we call θ n (t) and θ h (t). They will be discussed in Section IV. The second and fourth terms are cross-contributions accounting for thermo-elastic effects coming both from phonon relaxation effects (coefficients b 12 (t) and b 21 (t)) and phonon energy renormalization (coefficients a [2] 14 (t) and a [2] 23 (t)). The last terms I         44 (t), similarly to the cases of n(r, t) and h(r, t), are relaxation-type contributions corresponding to minus the inverse of (tensorial) characteristic Maxwell times, which we call θ [2] In (t) and θ [2] Ih (t), see Section IV. As in the equations for n(r, t) and h(r, t), we also find in the equations for the fluxes crossterms proportional to off-diagonal a [2]´s (due to self-energy corrections) and b [2]´s (due to relaxation effects), and external driving forces/sources. Furthermore, we call the attention to the fact that in Eq.(55) for the flux I n (r, t) of quasi-particles, the first term on the RHS (roughly proportional to ∇n) plays the role of a thermodynamic force analogous in classical fluid hydrodynamics to Fick's Law. In the same way, the first term on the RHS of Eq.(56) for the heat flux I h (r, t) (the term here roughly proportional to ∇h) is analogous to Fourier´s Law.
As a matter of fact, Eq.(55) and Eq.(56), under steady-state conditions and after neglecting both external forces/sources and thermo-elastic cross-terms, reduce to The D's play the role of rank-2 tensor diffusion coefficients.
In the next section we analyze several other aspects of this mesoscopic phonon hydrothermodynamics of order 1.

IV. CHARACTERISTIC MAXWELL TIMES AND MAXWELL-CATTANEO-LIKE HYPERBOLIC EQUATIONS
We consider here some additional characteristics which can be derived from the results of the previous section.

A. Characteristic Maxwell Times
In Eqs.(45)-(48) the four coefficients b 11 , b 44 are minus the reciprocal of the so-called Maxwell times [18], [47], namely In (t) In (q, t) /τ q , where Γ q is given in Eq. (25), b j (q, t) for j = 1, 2, b j (q, t) for j = 3, 4 and b ij (t) are given in Appendix D. 1/τ q = Γ q is the relaxation rate towards the equilibrium phonon population ν 0 q in mode q, which depends on q and ω q .
Eq.(59) to (62) tell us that the Maxwell characteristic times are given by a Mathiessen-like rule involving all the relaxation times in each mode, τ q , weighted by different kernels which are normalized, i.e., q w [2] In (q, t) = q w [2] Ih (q, t) = 1 [2] , where 1 [2] is the unit diagonal tensor.
As a consequence, if we assume all τ q to be independent of q (not a possible physical situation, see below) then the Maxwell characteristic times are all equal. On the other hand, if we apply the mean-value-theorem of calculus, taking outside the q a suitable mean-value In = 1 [2] τ In = 1 [2] /Γ (In) , θ [2] Ih = 1 [2] τ Ih = 1 [2] /Γ (Ih) . The quantity Γ q , given in Eq. (25), vanishes unless ω q = Ω k+q + Ω k due to the δ-function; then it can be rewritten identically as The matrix element |M kq | 2 behaves as ∼ |k + q|kq. The frequency ω (q) of LA phonons vanishes at the center of the Brillouin zone and is maximum at the boundary of the zone. Then inspection of Eq.(65) allows us to estimate that Γ q is an increasing function of |q| over the Brillouin zone, therefore, τ q = 1/Γ q is a decreasing one with q increasing.
The weighting functions in Eq.(59) to (62), and the time-dependent coefficients A ij and ∆ −1 ij below, are given in Appendix D, w [2] In (q, t) = ∆ −1 In a strictly Debye model it follows that where x D = βℏsq D , q D is the Debye wavenumber and s is the sound velocity, and It must be stressed that the equalities θ n = θ In and θ h = θ Ih are a consequence of using a In + θ −1 n ∂n(r, t)/∂t + −∇·B [2] 1 (t) · ∇ + θ −1 In θ −1 n + b 34 b 21 n(r, t) In these equations, the contributions associated with self-energy corrections were neglected and S n and S h account for sources and/or external forces.
With exactly the same procedures, we get for h(r,t) an equation of the same form, but where all n are replaced by h: It can be noticed that these hyperbolic equations resemble the so-called telegraphist equation in electrodynamics. Regarding the solution of Eq.(82) and (83) we can state that whenever the effective source has no spectral components at the frequencies of the hydrodynamic modes, there is a unique "particular" solution obtainable with a Green´s function.
Morse and Feshbach [45] give such Green function for the simpler case of a strongly localized effective source. However, if the effective source has spectral components at the frequencies of the hydrodynamic modes, "particular" solutions may still be found but they will in general diverge as |t| → ∞. for the phonon density n(r, t), the damping (terms with ∂n(r, t)/∂t) is related to sound attenuation in the material medium, which is very small in most rigid material media. In this case, the damping terms are a small correction on the undamped wave equation and can frequently be ignored.
Regarding Eq.(83) for the density of energy, suppose conditions of quasi-equilibrium have been attained. Then we can define a local quasi-temperature T * (r, t) such that h(r, t) = C V T * (r, t) where C V is a constant specific heat per unit volume. Then, for steady-state processes in ordinary solid material media, it is empirically established that ∂ 2 h(r, t)/∂t 2 is negligible in comparison with the damping terms ∂h(r, t)/∂t. In this case, Eq.(83) reduces to Fourier´s equation for T * (r, t) [46].
Let us consider heat motion in more detail. It is governed by Eq.(78) and (79).
Deriving Eq.(79) in time and taking the spatial gradient of (78), we eliminate h(r, t) and get the following hyperbolic equation for the heat flux I h (r, t): We recall that ∇ (∇·I h (r, t)) = ∇ 2 I h (r, t)+∇ × ∇ × I h (r, t). In that way, Eq.(84) is an extended version of the Guyer-Krumhansl equation, which in the steady-state and assuming ∇ × I h (r, t) = 0 reads as where, with ℓ 2 h having dimension of length. Notice that in the Debye model, C 2 is one third of the square of the sound velocity.

VI. THERMO-ELASTIC EFFECT
When a pulse of energy, well localized in space and time, is absorbed by a solid, the solid will heat and thermally expand, creating local time-dependent stresses and strains.
As time elapses, the absorbed energy will spread away from the "hot spot" and the solid will eventually reach a state of thermal equilibrium with a uniform temperature, as well as mechanical equilibrium with vanishing stresses and strains.
Of immediate practical relevance in this analysis are the scalar field T * (r, t) describing a local quasi-temperature, and the vector field u(r, t) describing time -and positiondependent displacements of the atomic nuclei.
Consider a homogeneous equilibrium solid with uniform density ξ 0 . If now the solid is locally excited in some way and the nuclei in some neighborhood undergo time dependent displacements, the density, ξ (r, t), will also become time -and position-dependent and for small displacements can be written as [47] ξ (r, t) It is clear that all motions of the nuclei can be expressed in terms of the normal modes of vibration in the solid, or, equivalently, in terms of the phonon density, be it crystalline or vitreous. Hence we argue that, upon a time -and space-localized excitation in a solid, the change in material density is, within a multiplicative constant, the same as the change in n(r, t), the phonon-density discussed in the previous sections. Therefore, we write ∂n(r, t)/∂t = −ξ 0 ∇ · ∂u(r, t)/∂t, ∂I n (r, t)/∂t = −∇·I [2] n (r, t) − θ −1 The five equations (88) through (92) can be combined to produce the hyperbolic equation ξ 0 ∂ 2 u(r, t)/∂t 2 + (ξ 0 /θ In ) ∂u(r, t)/∂t = −∇·I [2] n (r, t)+I [1]ext. n (r, t) .
To close this equation, we need to express ∇·I [2] n (r, t) in terms of a mixed representation, what is described in Appendix C. We end up with Moreover, as shown in Appendix C, the evolution equation for the nonequilibrium temperature (quasitemperature T * (r, t)) is Where, B n a h a n = q ℏω q ∇ q ω q · ∇ q ω q δν q /δF q q , ℏω q , δν q , /δF q , q δν q /δF q , (97) with ν q (t) given in Eq.(C11).
In the evolution equation for u(r, t) Eq.(94) we find in the RHS a term which is proportional to ∇T * (r, t) in agreement with the phenomenological equation of elasticity [47]. In the application to be considered here, there is no external source of "particle flux" and we set I [1]ext. n (r, t) = 0. In addition, the term ∂u(r, t)/∂t describes sound attenuation, which will be negligible in the illustration, hence we drop it.
In the evolution equation for T * (r, t), Eq.(95), there is an external source feeding energy in the solid medium, namely, a radiant energy pulse, so we keep the term I to be further specified below. Except at very short delays after excitation, the flow of heat is diffusive, which means that ∂ 2 h(r, t)/∂t 2 ≪ θ −1 h + θ −1 Ih ∂h(r, t)/∂t, therefore we drop the 2 nd order time-derivative. The effect of the 2 nd order time-derivative at short delay times was investigated in [50]. The term linear in h(r, t) describes time-relaxation effects which we lump together with the effect of the 1 st order time-derivative. We also neglect the term with n(r, t), because n is almost uniform and heat transport is weakly affected by small fluctuations in material density; in other words, the measured macroscopic diffusion coefficient D h already contains in itself the effect of density fluctuations which will be always present.
Given these considerations, we take as practical equations the following: where D h is the diffusion coefficient, C V is the specific heat per unit volume, w(r, t) is the power density being transferred from the radiant pulse into the medium, c s is the speed of sound, The pulses produced at FEL facilities are strongly localized in 3-dim space and in time. When such a pulse of X-rays reaches the interface vacuum/solid at the surface of an optical element, some fraction of the radiant energy is reflected and some fraction penetrates the solid to a depth δ 0 and is absorbed, generating heat and thermal distortion which impacts the optical performance of the device. We want to estimate the seriousness of the surface distortion in a time scale of pico-to nano-seconds after incidence of a single X-ray pulse lasting only a few femto-seconds. The analysis below will show that, although the incident femto-second pulse is gone long before the optical surface has time to distort, the next pulses in a pulse-train can be badly affected.
Consider the interface between vacuum at z < 0 and an infinite slab of silicon at z > 0, extending indefinitely in the x and y directions. Crystalline silicon is a very popular material for optical substrates, because it is available in large sizes, it has high thermal conductance and low coefficient of thermal expansion (hence thermal distortion is minimized), accepts state-of-the-art polishing (RMS surface roughness of only a few Angstroms) and is relatively cheap.
A short Gaussian-shaped pulse of X-rays (duration t F EL , radius r 0 ) centered at wavelength λ F EL and propagating along the z axis comes from z = −∞ and hits the interface at t = 0. This scenario is cylindrically symmetric about the z axis. Tables I and   II describe the FEL pulse and the solid medium.  The power density w(r, t) to be used in Eq.(99) is if z < 0, for any time −∞ < t < ∞, if z > 0, for any time −∞ < t < ∞, For simplicity we take the limit of "short" z 0 and t F EL , with (z 0 /t F EL ) = c, the speed of light in vacuum. Integrating P over all time and all space we get the total energy W F EL = U 0 π 3/2 r 2 0 z 0 carried by the pulse, and find that U 0 has the meaning of average energy density in the radiant pulse. µ is the light absorption coefficient in the solid medium, which, for a given medium, depends on the wavelength λ F EL of the radiation. Notice that the calculation is linear and all results scale linearly with the amount of absorbed energy, which in this application is 0.16 mJ.
The procedure here will be the following. First, one solves Eq.(95) for the temperature field T * (r, t). Next, one uses ∇T * (r, t) as source in Eq.(94) for u(r, t).
Green's functions, analytical expressions for the solutions T * (r, t) and u(r, t), and a discussion of numerical methods are given in Ref. [48]. Here, let us recall that u(r, t) = u P art (r, t)+u F ree (r, t) where u P art is a particular solution which depends on the given source, while u F ree is any arbitrarily chosen solution of the associated homogeneous equation (no sources), chosen to satisfy boundary/initial/asymptotic conditions, as the case may be. The condition to be met here is that all normal stresses at the "free" surface z = 0 of the solid medium be vanishing: σ zi (z = 0) = σ iz (z = 0) = 0. These conditions lead to coupled Fredholm integral equations of 1 st kind, which we have solved only approximately using truncation, see Ref. [48].
The result of these thermo-elastic calculations is that, on absorption of 0.16 mJ of energy from a Gaussian X-ray photon pulse at ℏω = 5000 eV , with r 0 = 2.0 mm, lasting 10 f sec, there is an outwards surface bulge several nm high that comes after the FEL pulse, with a delay of several hundred nanoseconds.
Detailed results are shown in Figures 1 to 5. Figure 1 shows the temperature T * versus depth z inside the Silicon slab, at selected times after incidence of the X-ray pulse, assuming T * = 0 o C as initial temperature. The thick vertical line is the "causal cut-off" at t = 100 f sec; namely, for t = 100 f sec, the light  penetrates only as far as ct = 2.998× 10 4 nm, hence the temperature at |z| > ct is still identically zero. The cutoff for the other is off-range in this figure. However, as t increases, the X-ray pulse penetrates deeper and deeper till it is depleted by absorption. Then, further heating of Silicon layers far away from the surface proceeds by diffusion only.   the profile at 400 nsec we find a surface figure error (maximum slope of the surface) of about 2 µrad, which is not small, and takes microseconds to decay. Notice that the first X-ray pulse hitting "cold" silicon surface suffers no adverse effects, because at the time the surface bulges out, the pulse is already long gone. If, however, the experiment envisages a train of FEL pulses, and the spacing in the train is a few hundred nsec, the pulses following the first will be defocused by the heat-induced surface bulge. Furthermore, the bulge can be resonantly enhanced in a disastrous way.   We can compare the present results to the previously studied case of vacuum ultraviolet light [48], [49], [50] where the photon energy was ℏω = 100 eV , absorbed energy per unit area W 0 = 180 µJ/cm 2 , absorption coefficient 0.128 (nm) −1 , peak displacement u M ax = 0.02 nm and "specific displacement" u M ax /W 0 = 1.3 × 10 −4 nm/(µJ/cm 2 ).
For X-rays, present calculation, we find u M ax /W 0 = 8.4 × 10 −3 nm/(µJ/cm 2 ), which is 65 times larger, even though the absorbed energy is comparable. However, if µ is the wavelength-dependent absorption coefficient, the penetration depth is of order 1/µ. This is the thickness of the layer that is strongly driven out of equilibrium by absorption of the radiant energy in the light pulse. It is 7.8 nm for vacuum ultra-violet light (ℏω = 100 eV ), but 1.8 × 10 4 nm for X-rays (ℏω = 5000 eV ). We conclude that deep penetration of light wrecks havoc with the surface of optical substrates.

VII. CONCLUDING REMARKS
We have presented the description of a complete Mesoscopic Hydro-Thermodynamics of phonons, which can also be referred to as Higher-Order Nonlinear Generalized Hydrodynamics of phonons. It significantly extends the standard hydrodynamics of phonons by introducing a complete description in terms of the densities of phonon quasi-particles and of the energy, accompanied with their fluxes of all orders [Cf. Eqs. (10) and (15)]. That is, as indicated in the main text, we can talk of the motion of two families, the one associated to the motion of the quasi-particle density, together with the evolution equations for its fluxes of all orders, and the one for the motion of energy density, together with the evolution equations for its fluxes of all orders [Cf. Eqs. (29) and (30)], coupled together by cross-correlations describing thermo-striction effects.
As already noticed in the introduction MHT allows to cover all kinds of motion, in that it includes those characterized by short wavelengths and ultrafast time evolution. The system of coupled equations of motions is extremely cumbersome, in principle of unmanageable proportions. For handling it is required, depending on each case being considered, to introduce a contraction of description implying in retaining only a few number of fluxes, neglecting those that become negligible in time intervals smaller than the experimental resolution time.
In other words [43], the contraction of description implies in retaining the information considered as relevant for the problem in hands, disregarding nonrelevant information. In the main text (Section 2) it has been discussed a criterion for deciding on the order of contraction. It has been considered in full detail the case of a MHT of order 1, which corresponds According to NESEF ( [12], [14], [19] with a short overview given in [15]), the nonequilibrium statistical operator in terms of the basic nonequilibrium variables in sets (5) and (7) is given by where with ̺ (t, 0) being the auxiliary statistical operator (also called "instantaneous quasiequilibrium operator") and where t ′ is the dependence on time of the non-equilibrium thermodynamic variables F and the dynamic microvariables, in the Heisenberg representation, depend on (t ′ − t). Moreover, ̺ B is the canonical distribution of the bath of acoustic phonons in equilibrium at a temperature T 0 , and φ(t) ensuring the normalization condition plays the role of the logarithm of a nonequilibrium partition function.
We recall that the second term in the exponent in Eq.(A2) accounts for historicity and irreversibility in the non-equilibrium state of the system. The quantity ε is a positive infinitesimal that goes to zero after the trace operation in the calculation of averages has been performed. We also recall that i.e., it has an additive composition property, with a contribution of the instantaneous quasiequilibrium statistical operator plus the one of ̺ ′ ε (t) which contains the historicity and produces irreversible evolution.
Next, consider a change of description consisting into going from the one in terms of the set (5) to one in terms of the (hydrodynamic-in-character) basic microdynamical variables n(r, t), I n (r, t), ..., I [ℓ] n (r, t), ..., h(r, t), I h (r, t), ..., I whose average values, taken over the non-equilibrium ensemble, of these operators, are the macrovariables in sets (10) and (15), which we have dubbed as the n-family and the h-family respectively. In order to calculate the averages we use Eq.(A.3) once we write where ∇ q ω q is the group velocity of phonons with crystal momentum q and ∇

[ℓ]
q ω q is given in Eq. (14). The symbol ⊙ stands for fully contracted product of tensors.
"Contraction of description" in a given order, say k, is done by taking as null the quantities F h (r, t) for all ℓ > k. In the main text we have introduced a study of MHT of order 1, i. e., keeping only ϕ n (r, t), ϕ h (r, t), F n (r, t) and F h (r, t), and the closure of the evolution equations was done using the Heims-Jaynes [44] method as described in Appendix C.

Appendix B: Generalized Peierls-Boltzmann Equation
As indicated in Eq. (22), the evolution equation for the single particle distribution Applying the NESEF-based kinetic theory we find where, Here, δ stands for functional derivative. We stress that this expression corresponds to an approximation where the interactions H SB and H SP are retained only up to second order (memory and vertex renormalization effects are neglected). After performing the calculations it follows that for Q = 0, and for Q = 0, where ν 0 q is the distribution in equilibrium at temperature T 0 . Moreover Γ q , the reciprocal relaxation time of the phonons in mode q, is given in Eq. (25) in the main text, and Π q , the self-energy correction of the energy ω q , is given in Eq. (28).
For the sake of simplicity, we introduce an approximation of long wavelength, i.e., we consider Q << Q B , Q B being the Brillouin radius, in the form Then, when going over to direct space, through the Fourier transform of variable Q into r, there follows Eq.(23) which has a form resembling the standard Peierls-Boltzmann equation.

Appendix C: Application of Heims-Jaynes Formalism
Writing for the auxiliary statistical operator (cf. Eq.(A3) where that is, a separation in terms of the homogeneous contribution, A, and the departure from it, B, and introducing the homogeneous contribution, according to Heims-Jaynes perturbative expansion for averages [44] keeping only the first order contribution in the inhomogeneous part of B, and using the contraction in MHT of order one, i. e. [see Eq. (37)] where with F q (t) = ϕ n (t) + ϕ h (t) ℏω q + F n (t) ·∇ q ω q + F h (t) ·ℏω q ∇ q ω q , and we notice that ν q (t) [1 + ν q (t)] = δν q (t)/δF q (t) .
Using the nonequilibrium equations of state, after going over direct space we arrive at the expressions for the divergence of both second-order fluxes given in Eqs. (49) and (50), and then to the closed system of Eqs.(53) to (56).
On the other hand, introducing the concept of nonequilibrium temperature, better called quasitemperature T * (r, t) in the form we can obtain an evolution equation for it starting with the evolution equation for the energy in the form of the hyperbolic Maxwell-Cattaneo equation, Eq.(75), from which together with the nonequilibrium thermodynamic equation of state, Eq.(C14), we have that q (ℏω q ) 2 ν q (t) [1 + ν q (t)] 43 (t) 44 (t) · ∇ϕ h (r, t) and, after introducing the heat capacity where T 0 is the temperature in equilibrium in this linear treatment, we arrive at Eq.(95).