Steady State and Modulated Heat Conduction in Layered Systems Predicted by the Phonon Boltzmann Transport Equation

Based on the phonon Boltzmann transport equation under the relaxation time approximation, analytical expressions for the temperature profiles of both steady state and modulated heat conduction inside a thin film deposited on a substrate are derived and analyzed. It is shown that both steady state and modulated components of the temperature depend strongly on the ratio between the film thickness and the average phonon mean free path, and they exhibit the diffusive behavior as predicted by the Fourier law of heat conduction when this ratio is much larger than the unity. In contrast, in the ballistic regime when this ratio is comparable to or smaller than the unity, the steady-state temperature tends to be independent of position, while the amplitude and the phase of the modulated temperature appear to be lower than those determined by the Fourier law. Furthermore, we derived an invariant of heat conduction and a simple formula for the cross-plane thermal conductivity of dielectric thin films, which could be a useful guide for understanding and optimizing the thermal performance of the layered systems. This work represents the Boltzmann transport equation-based extension of the Rosencwaig and Gerko work [J. Appl. Phys. 47, 64 (1976)], which is based on the Fourier law and has widely been used as the theoretical framework for the development of photoacoustic and photothermal techniques. This work might shed some light on developing a theoretical basis for the determination of the phonon MFP and relaxation time using ultrafast laser-based transient heating techniques.

It has been shown that the phonon Boltzmann transport equation (BTE) is a more appropriate tool to describe the transport phenomena in nanostructured materials and during the ultrafast processes. 4,6 Though great progresses have been made in solving the BTE for micro/nanoscale heat conduction with significant efforts in recent years, the inherent difficulties associated with its solution have significantly limited the consideration of the size and transient effects. Under the steady-state conditions, the BTE has been solved numerically and applied to study the heat transport through a variety of layered systems and complex geometries. 1,4,5,11,12 These works showed the two main findings: i) a reduction of the effective thermal conductivity with respect to their bulk values, and ii) the temperature profiles differ significantly from those obtained using the Fourier's law, due to the ballistic behavior of the energy carriers. The transient heat conduction in thin films has also been revisited in some recent works using the transient phonon BTE. [13][14][15] Taking into account that the energy carriers travel ballistically without being deflected out of their propagation direction in a spatial scale in the order of one MFP, Chen proposed the ballistic-diffusive equations to study transient heat conduction from macro-to nanoscales. 13,14 Even though this model presents good agreements with the predictions of the BTE for the heat conduction, it cannot be easily implemented. More recently Ordonez-Miranda et al. 15 [20][21][22][23][24][25] Clearly there is a strong need of phonon BTE solutions to better understand and validate data reduction schemes in these experiments which at present are dominantly through fitting an effective thermal conductivity or conductance with the Fourier's law. 18 The objective of this work is to solve analytically the phonon Boltzmann transport equation for the temperature profile in a dielectric film deposited on a substrate, when excited with a modulated laser beam. The obtained results for the steady state temperature distribution depend strongly on relative thickness of the film with respect to the phonon MFP and they agree quite well with the numerical results previously reported in the literature. 1,6 The amplitude and the phase of the modulated temperature, on the other hand, depend strongly on the product of the modulation frequency and the phonon relaxation time, such that they reduces to the corresponding results predicted by the Fourier's law when this product approaches zero. This work represents the BTE-based extension of the work by Rosencwaig  Let us now consider a two-layer system shown in Fig. 1(a), where a dielectric film with a finite thickness is heated by a laser beam with modulated intensity, at the surface z 1 = 0 . This system represents a thin film deposited on a substrate in analogue of sample configuration for photo-thermal characterization of thermal conductivity. The surface heat flux q at z 1 = 0 by the external thermal excitation is given by 26 where ω is the angular modulation frequency, 2q 0 is the average intensity of the laser, t is the time and Re() stands for the real part of its argument. The steady state and modulated temperature fields within the layers, for the surface heat flux defined in Eq.(1) can be easily found using the Fourier's law and other revised Fourier-like models of heat conduction. 2,27,28 However, considering that the MFP l 1 of the energy carriers inside the film can be comparable to or even smaller than the thickness L (l 1 < L), these macroscopic models might not be valid. 6 We hereby develop analytical solutions using the transient phonon BTE.
For the sake of simplicity and to present the results analytically, this work will be developed assuming that the scattering processes inside the layers are described by the average values of the MFP ( l ) and relaxation time (τ ) of phonons, such that these parameters are independent of the phonon frequency. Assuming that the changes of temperature due to heating is small, the phonon BTE, under the single mode relaxation time approximation and in the intensity representation, can then be written as 1 0 , where µ = cos(θ ) is the cosine of the angle between the phonon propagation direction and the z + axis, I is the total phonon intensity, and 0 I is the "equilibrium" phonon intensity defined by 1,5 where ρ is the mass density and c the specific heat at high temperature. 6 Taking into account the definition and the temperature dependence of the specific heat c , 8,10 Eqs. (3) and (5) shows that 0 I can be written in the form of the second equality in Eq. (6), i.e., linear relationship between the equilibrium intensity and temperature, for both low and high temperatures. We will use this relationship that the equilibrium intensity is proportional to the temperature to express some results of the present work.
If the intensity I in the one-dimensional phonon BTE is solved, the heat flux q can then be determined as follows 15 For the surface thermal excitation defined in Eq.(1), the solution for the intensity I of Eq.
(2) has the following general form where s I and J are the stationary and modulated components of I , due to the first and second terms of the external thermal excitation, respectively. Given that Eq. (2) is a linear partial differential equation, Eq.(8) implies that the equilibrium intensity 0 I should have the same time dependence as the total intensity I . Therefore, the general form of 0 I is Furthermore, given the direct proportionality between the equilibrium intensity and temperature for small temperature changes as articulated earlier, Eqs. (6) and (9) establish that where T (z) and ψ (x) are the steady state and modulated components of temperature. After inserting Eqs. (8) and (9) into Eq.(2), the following two uncoupled differential equations can be written for s I and J , respectively: where x = z l and 1 i χ ωτ = + . Note that the frequency-dependence of J comes through the complex parameter χ . When 0 ω = , Eq.(10b) reduces to its stationary counterpart Eq. (10a).
Given that the typical values of the average relaxation time are τ < 10 −10 s for a wide variety of materials at room temperature, 6 the parameter χ → 1 , for any periodic heating with frequency f = ω 2π << 1 GHz . This value indeed covers almost all the operating frequency range of existing photo-thermal and photoacoustic techniques. For these reasons, we are going to derive the solutions for the modulated heat transport with the assumption χ ≤ 1. The solution of Eq.
(10b) is straightforward and is given by where x 0 is an integration constant. Due to the interface roughness, 6 usually diffuse interface scattering dominates and we will limit our discussion to diffuse interface scattering in this work.
For diffuse surface/interface scattering, the intensities leaving the surfaces at z = 0, L are uniform, 11 and therefore the coefficients J (x 0 ,µ) should also be independent of the direction µ .
To facilitate the evaluation of the boundary conditions, it is convenient to split Eq. (11) into two parts, one for each layer of the system shown in Fig. 1(a), as follows 5 where λ = L l 1 , χ n = 1+ iωτ n , the superscript ( ) shown in Fig. 2(b), the constants A 1 ± and A 2 + are determined by the boundary conditions, and the subscripts n = 1 and n = 2 stand for the first and second layers, respectively. The temperature field is determined by the principle of energy conservation, which in terms of the phonon intensities I 0 , can be written as 1 For the modulated components of the intensities J 0 and J , Eq. (13) yields After inserting Eqs. (12a)-(12d) into Eq. (14), the following equations for J 0 inside the first ( J 01 ) and second ( J 02 ) layers are obtained where E n () is the exponential integral function of order n. 29 Based on Eqs. (7) and (8), the modulated component q t of the heat flux in each layer can be written as The combination of Eqs. (12a)-(12d) and Eq. (16) yields the following modulated heat fluxes inside the first ( q t1 ) and second ( q t 2 ) layers The constant A 1 ± and A 2 + can be determined by imposing the energy balance at the interfaces of the layers. According to Eq. (7) and Fig. 2(b), the conservation of energy at the inner interface of the layers establishes that where ij r and ij t are the energy reflectivity and transmissivity of phonons coming from the i th layer toward the j th layer, respectively. Considering that under the diffuse scattering, the scattered phonons completely lose their memory, 5 these coefficients are direction-independent and satisfy the relations 5,11 On the other hand, the energy balance at the illuminated surface which is a particular case of Eq. (16). The combination of Eqs. (12), (18), and (20) yields the following system of equations for the constants A 1 ± and A 2 Under the boundary conditions in Eqs. (21a)-(21c), Eqs. (15) and (17) along with Eq. (4), determine fully the modulated heat transfer through a finite layer on a semi-infinite substrate.
The steady-state component of the equilibrium intensity I 0s and the heat flux q s can be easily determined by replacing J 0n → I 0sn , q tn → q sn , and χ n = 1 in Eqs. (15), (17), and (21). Below we analyze separately the steady state and modulated heat conduction problems.

A. Steady-state heat conduction
Under steady state condition, Eqs. (15) and (17) can be normalized as follows with A 2 − = 0 . By taking the derivative of Eqs. (22c) and (22d) and comparing the results with the corresponding equations (22a) and (22b), it can be shown that Q s1 and Q s2 are constants. This is consistent with the principle of energy conservation, which establishes that q s1 = q s2 = q 0 . In what follows, we are going to obtain and discuss the analytical solutions of Eqs. (22a) and (22b).

• Semi-infinite layer
Based on the properties of the function E n () , 29 it can be seen that the unity is a particular solution of Eq. (22b). Therefore the general solution for Eq. (22b) can be written as In terms of the function G(x) , Eqs. (22b) and (22d) reduce to Given the direct proportionality between the equilibrium intensity and temperature, from the Fourier's law of heat conduction, we know that in the diffusive limit ( where a and b are constants. In this limit, we insert the function G in Eq. (20b) and evaluate explicitly the involved integrals to obtain b = 3Q s2 4 . Therefore, the general solution of Eq.
(24a) can be written as where A is a constant and the function g(x) → 0 for x → ∞ . The combination of Eqs. (24a) and (25) yields the following integral equation for g(x) which suggests that for a first-order approximation, g(x) is given by The constants A , B , and C can be calculated by evaluating Eqs. (24a) and (24b) at x = 0 : The required third equation can be obtained by expressing the right-hand side of Eq. (24b) as an exact derivative and integrating on both sides. The final result is: where δ is an integration constant. The evaluation of Eq. (28) in the diffusive limit ( x → ∞ ) By inserting Eqs. (25) and (26b) into Eqs. (27) and (29), the following system of equations for A , B , and C is obtained where the coefficients I nm are defined by which can now be calculated analytically. 29 The solution of Eq. (30) for five decimal figures is A = 0.71047 , B = −0.25082 , and C = 0.23526 . Hence, the approximate solution of Eq. (29a) can be written as By comparing Eq. (32a) with the solution for a semi-infinite heat conduction problem using the Fourier's law, it is clear that the non-Fourier (ballistic) contribution to the temperature is contained in p(x) , which is a positive function that increases monotonously with x = z l , as shown in Fig. 3. The minimum and maximum values for p(x) are p(0) = 1 3 and p(∞) = A , respectively. The first of these values has been obtained using the exact values of A , B , and C .
A second-order approximation for f (x) can directly be obtained by inserting Eq. (32) into the right-hand side of Eq. (24a). The mathematical expression of this iterated solution is much longer than Eq. (32a), but in numbers they differs in less than 0.1% for any x ≥ 0 . This fast convergence for the solution of f (x) is reasonable given that the exponential decay of the integral exponential functions E n (x) as x increases. The fast convergence of expansions involving the functions E n (x) was also found and verified in radiative heat transfer. 30 Hence, for applications of practical interest, Eq. (32a) provides a good approximate solution for the In absence of the layer with finite thickness in Fig. 2(a), and considering that the temperature at the illuminated surface of the semi-infinite layer is T 0 , Eqs. (6) and (32) establish that the explicit expression for the temperature field T inside this semi-infinite layer is which differs from the Fourier's law prediction by the function 1 For a position inside the semi-infinite layer with a distance much larger than the phonon MFP from the illuminated surface, i.e., x = z l >> 1, Eq. (33) reduces to the diffusive prediction, as should be.

• Thin Film Layer
Here we derive an explicit expression for the steady state temperature inside a thin film of thickness L. This will allows us comparing our analytical approach with previous numerical results reported for this system. 1,6 Based on Eq. (A19) of the appendix, the solution of Eq. (22a) for a thin film layer is given by To find the parameters β and γ , we evaluate Eqs. (22a) and (22c) at x = 0 : The combination of Eqs. (34) and (35) yields the following system of equations where (37) Figure 4 shows the behavior of the parameters β , γ and of the normalized heat flux Q s1 as a function of the normalized layer thickness λ = L l 1 . Note that the values of these three parameters β , γ Q s1 are bounded, such that β and γ reach their maxima when λ = L l 1 is close to 1. For a semi-infinite layer ( λ → ∞ ), β → p(∞) , γ → 1, and Eq. (34a) reduces to Eq.
(32a) of the semi-infinite layer. On the other hand, for a very thin layer with λ << 1, the normalized heat flux reduces to Q s1 = 1, which agrees with Eq. (35b). Assuming that the layer surfaces x = 0,λ are at temperatures T 1 and T 2 (< T 1 ) , then this last condition renders that the heat flux across the layer is where σ is the analogous Stephan-Boltzmann constant for phonons as defined in Eq. (4). The results at the ballistic limit is very different from the prediction of Fourier's law in the diffusive regime and coincide with the result reported by Swartz and Pohl,9 in absence of phonon scattering within a film. Furthermore, for an arbitrary λ > 0 and assuming that the temperature difference T 1 − T 2 is small enough, Eqs. (6) and (34b) establishes that the heat flux across the layer is given by has the form of the Fourier's law but with a modified thermal conductivity where k 0 = ρcvl 3 is the bulk thermal conductivity of the layer. For λ << 1, β → 2 3 and Eq.
(35) reduces to the numerical result reported by Majumdar 1 and Chen. 6 Given that Eq. (38) has been rigorously derived from an analytical approach, this result represents an accurate extension of those previous results for the cross-plane thermal conductivity of a thin film. • Two-layer system For the two-layer system shown in Fig. 2(a), according to Eq. (34a), the equilibrium intensity I 0 s1 inside the finite layer can be written as Furthermore, based on Eqs. (23a) and (32a), the intensity I 0 s 2 inside the semi-infinite layer is given by The remaining constant A 1 + is set by the specification of the temperature at the other side ( x → ∞ ) of the semi-infinite layer. Given that this value is usually not known explicitly, it is convenient to replace this condition by the temperature T 0 at the illuminated surface ( x = 0 ) of the finite layer. Therefore A 1 + = ρ 1 c 1 ν 1 T 0 4 and Eqs. (39a) and (39b) can be written in terms of the temperature as where we used the fact that δ 12 = ρ 1 c 1 ν 1 ρ 2 c 2 ν 2 = r 12 t 12 , as established by Eq. (19). Equations (41a) and (41b) describe fully the temperature inside the finite and semi-infinite layers shown in Fig. 2(a), by taking into account the effects of the film thickness and of the thermal mismatch between the layers ρ 1 c 1 ν 1 ρ 2 c 2 ν 2 .

B. Modulated heat conduction
In this subsection, we are going to solve Eqs. (15) and (17) for the modulated heat conduction inside the finite and semi-infinite layers shown in Fig. 2(a). For simplicity, we are going to start with the semi-infinite layer in absence of the finite one, as we did for the steadystate problem. •

Semi-infinite layer
Based on Eqs. (15b) and (17b), the equilibrium intensity and heat flux inside the semi-infinite layer can be normalized as where V 2 (x) = π J 02 (x) A 2 + . In the diffusive limit ( x → ∞ ), the solution of the Fourier's law indicates an exponential decay for the temperature V 2 (x) = χ 2 aexp(− χ 2 bx) and the heat flux where the constants A 0 and B 0 can be calculated by evaluating Eqs. (42a) and (42b) at x = 0 : The combination of Eqs. (43) -(45) yields the following system of equations 2 − I 12 1− I 13 where the coefficients I nm are defined in Eq. (31).
The solution of Eq. (46) for five decimal figures is A 0 = −3.97433 and B 0 = 9.17085 .
Thus, the solution of Eq. (42a) in terms of J 02 (x) = A 2 and the corresponding modulated temperature (see Eq. (6)) and heat flux (see Eq. (42b)) are which coincides with its usual value provided that the thermal conductivity k = ρcvl 3.

• Thin Film Layer
According to Eqs, (A19) and (A31) in the Appendix, the modulated equilibrium intensity J 01 and modulated heat flux q t1 inside the finite layer shown in Fig. 2(a) are given by where a 1 ± = (4η 1 3q 0 ) A 1 ± and the functions M and M 1 are defined in Eqs. (44) and (A32), respectively. The constants A, B, and C can be determined in terms of a 1 ± by evaluating Eqs.

Two-Layer System
The modulated components of the temperature and heat flux inside the two-layer system shown in Fig. 2(a) are derived in this subsection. These results are then compared with those predicted by the Fourier's law for the diffusive regime. According to our previous analysis for a thin film layer in the semi-infinite limit, the modulated equilibrium intensity J 02 and modulated heat flux q t 2 inside the semi-infinite layer in thermal contact with the finite layer shown in Fig.   2(a) can be written as follows where a 2 + = (4η 2 3q 0 ) A 2 + . After combining Eqs. (48a) and (52a) with the boundary conditions in Eqs. (21a)-(21c), the following system of equations is obtained for the constants a 1 where we have used the fact that a 2 + = τ 2 τ 1 a 1 − δ 12 (Eq. (21a)), and the parameters a ± , b ± , and c ± are given by Eq. (50), such that A = a + a 1 + + a − a 1 − , (50)  where ψ n (x) = 4π J 0n (x) ρ n c n v n .

III. RESULTS AND DISCUSSIONS
In this section, both the steady state and modulated temperature profiles are analyzed as a function of the film thickness and the modulation frequency. For the two-layer system shown in Fig. 2(a), calculations were done with the data reported in Table I. and (41b)), respectively. As expected, when the distance from the excited surface z = 0 of the semi-infinite layer increases to values much longer than the phonon MFP, the temperature predicted by the BTE reduces to the straight line as predicted by the Fourier's law. By contrast, as this distance takes any smaller values than a phonon MFP, the temperature determined by BTE tends to the constant value T B = T 0 − 3q 0 ρcv , which differs from the predictions of the Fourier's law. This is reasonable given that the emitted phonons at z = 0 travel ballistically within an average MFP distance, and therefore the temperature does not change remarkable within a phonon MFP distance. 1,6 The fact that T 0 > T B indicates that the temperature imposed at the outer surface z = 0 + is higher than the emitted phonons at the inner surface z = 0 − . Similar to that in the semi-infinite medium in Fig. 5(a), a temperature jump is also observed in Fig. 5(b) at the boundaries of a finite layer. These jumps increase as the layer thickness reduces. This is reasonable considering that the "hot" phonons emitted from the surface z = 0 heat up the "cold" phonons emitted from the opposite surface z = L and vice versa. In the extreme case of ballistic limit ( L << l 1 ), the temperature inside the layer is given by (T 1 + T 2 ) 2 .

A. Steady-state temperature profiles
It is clear from Fig. 5(b) that the analytical results of the present work are in very good agreement with the dotted lines obtained through numerical simulations of phonon Boltzmann transport 1,6 The advantage of using the analytical approach over the numerical one is the simple description of the temperature through Eq. (34a), which is able to provide physical insights without much of numerical efforts. As the layer thickness is scaled down to values comparable to the phonon MFP, the temperature tends to become independent of the position, which is the signature of the ballistic heat conduction. The temperature profiles and the jumps at the interfaces are very similar to those observed for photon radiative transfer in a plane-parallel medium. 31 The combined behavior of the temperature inside the semi-infinite and finite layers arises in the two-layer system, as shown in Fig. 5(c).  thickness, and (c) the two-layer system shown in Fig. 1(b). Dotted lines in (b) correspond to the numerical predictions reported by Chen,6 and the calculations in (c) were performed with the data of Table I for a film of Si of thickness L = l 1 deposited on a semi-infinite Ge substrate. The predictions of the Fourier's law in Fig. 5(c) were calculated using the interface thermal ) , which is appropriate for the diffusive interface scattering. 5 The cross-plane thermal conductivity defined in Eq. (38) for a single layer is shown in Fig. 6, as a function of its normalized thickness. This result agrees well with the work by Majumdar 1 in the pure ballistic ( L << l 1 ) regime, but they differ in the intermediate diffusive-ballistic regime.
This is because the parameter β involved in Eq. (38) is different than its ballistic value 2/3, within this thickness range, as shown in Fig. (4). Given that our formula for the thermal conductivity has been derived rigorously from the analytical solution of the phonon BTE, its predictions are expected to be an extension of Majumdar

B. Modulated temperature profiles
The normalized amplitude and phase of the modulated temperature at the illuminated surface z = 0 of a semi-infinite layer are shown in Figs. 7(a) and 7(b), respectively. We have considered the frequency interval ωτ ≤ 0.1 , because the analytical result in Eq. (47b) is only valid for ωτ << 1. In the low frequency regime (ωτ → 0 ), the amplitudes and phases predicted by the BTE and Fourier's law agree with each other, as expected. Given that the average relaxation time of phonons in most materials is shorter than hundreds of picoseconds (τ < 10 −10 s ), 8

this fact
indicates that the Fourier's law holds when the frequency f = ω 2π << 1 GHz . On the other hand, for frequencies close to but smaller than ω = 0.1τ −1 , the predictions of the BTE and Fourier's law differ, especially on the phase signal. As the frequency increases, the phase predicted by the Fourier's law remains independent of the frequency, while the one predicted by the BTE decreases. Taking into account that the BTE is much more general than the Fourier's law, the modulated heat conduction at high frequency (f > 1 GHz) is expected to be better described by the BTE.  z → ∞ , which create a zone of ballistic phonons ( 0 < z ≤ l ), whose modulated temperature is slightly smaller than the one they would have in equilibrium.   Table I.  Table I.

IV. CONCLUSIONS
Steady state and modulated components of the temperature inside a dielectric film deposited on a semi-infinite substrate have been analyzed by solving analytically the phonon Boltzmann transport equation under the relaxation time approximation. It has been shown that both the steady state and modulated components of temperature depend strongly on the ratio between the film thickness and the phonon mean free path, and they reduces to their corresponding diffusive values predicted by the Fourier's law as this ratio increases. By contrast, in the ballistic regime in which this ratio is comparable to the unity, the steady state temperature tends to be independent of position, and the amplitude and phase of the modulated temperature display lower values than those determined by the Fourier's law, such that their corresponding difference increases with modulation frequencies higher than 1 MHz. Furthermore, an invariant of heat conduction and a simple formula for the cross-plane thermal conductivity of dielectric thin films are obtained, which could be useful for understanding and optimizing their thermal performance. This work could serve as the theoretical basis for the determination of phonon mean free path using transient thermoreflectance method.

Appendix: Analytical Solution of the BTE Obtained using the Discrete-Ordinates Method
To obtain the solution of the phonon BTE, we are going to use the discrete-ordinates method, 33 which has widely been used to solve numerically the Boltzmann transport of photons, neutrons and most recently phonons. In this appendix, we are going to show that this method can also provide analytical solutions for both the steady state and modulated components of the phonon temperature. This method is based on the Gaussian quadrature where the coefficients a n = a − n are known for n = ±1,±2,... ± N . Thus, the quadrature method consists of splitting the interval −1 ≤ µ ≤ 1 to 2N symmetrical directions ( µ − n = −µ n ). For • Steady-state heat conduction problem According to Eqs. (10a), (13) and (A1), the steady-state BTE for the phonon intensity I s is given by (A3) yields where 2D = a n n ∑ C n is a constant independent of n . Inserting again the exponential solution into Eq. (A3) and using Eq. (A4), the values of the exponent δ are determined by the following characteristic equation, According to Eq. (A2) for i = 1 , Eq. (A5) can be conveniently rewritten as a n µ n 1− µ n δ = a n µ n which is polynomial of degree N − 1 in δ 2 . Therefore, the 2(N − 1) roots of Eq. (A6) have the form δ = ±δ j for j = 1,2,..., N − 1 . Thus, the first 2(N − 1) solutions of the system in Eq. (A3) are determined by where D j and D − j are arbitrary constants. Given that the system of equations in Eq. (A3) is of order 2N , we require two more independent solutions. Note that Eqs. (A2) and (A5) can be written as a n n=1 N ∑ = 1 = a n 1− (µ n δ ) 2 n=1 N ∑ , which clearly indicate that δ = 0 is a double root of the characteristic equation. Based on the theory of differential equations, the intensity corresponding to this double root is given by I s (x,µ n ) = α n ' (x + β n ' ) , where α n ' and β n ' are constants. The substitution of this expression into Eq. (A3) yields α m ' = a n α n ' 2 n ∑ = α ' and β m ' + µ m = a n β n ' 2 n ∑ = β ' . The general solution of the system in Eq. (A3) can therefore written where the constants D j and D − j appearing in Eq. (A7) have been redefined to absorb α ' . Note that Eq. (A8) is valid for any −1 ≤ µ n ≤ 1 and therefore it can also be written with the substitution µ n → µ . Equation (A8) shows clearly that the steady-state intensity has both diffusive and ballistic contributions provided by its linear and exponential dependence on the position, respectively. For phonon heat conduction across a layer with diffuse scattering at its boundaries, Eqs. (11) and (12) establishes that for 0 < µ ≤ 1 , the boundary conditions of I s are given by where λ = L l is the ratio between the layer thickness and the MFP of phonons inside the medium. After inserting Eq. (A8) into Eqs. (A9) and (A10) and summing the obtained results, the following is obtained The intensity I s is now expressed as Based on Eqs. (7), (13)   Eq. (A16) can be considered as the series expansion of a given function and therefore they can be written as follows where γ is a constant. Equation (A18) can then be expressed in a more meaningful way in term of the normalized equilibrium intensity U = π I 0s − A − is obtained to be valid for any position 0 ≤ x ≤ λ . Equation (A22) establishes that the sum of the equilibrium intensities at two equidistant points from the external surfaces of a finite layer, is an invariant of heat conduction. For temperatures smaller than the Debye temperature, this indicates that if we impose the temperatures T 1 and T 2 at the external surfaces of a finite layer as shown in Fig. A1, then A + = πσ T 1 4 , A − = πσ T 2 4 , and T 4 (x) + T 4 (λ − x) = T 1 4 + T 2 4 . In particular, T 4 (λ 2) = (T 1 4 + T 2 4 ) 2 is obtained at the center of the layer, which is the average of the fourth powers of the temperatures at the external surfaces. This feature of temperature is an exclusive prediction of the phonon BTE, which cannot be predicted by the Fourier's law. An analogous behavior for temperature was found in the radiative heat transfer. 30 1 A1. Temperature distribution inside a finite layer.

• Modulated heat conduction problem
Similar to the steady-state case, we start with the numerical form of the phonon BTE for the modulated component of the temperature The characteristic equation of Eq. (A22) can be derived by following a similar procedure as the one done for the steady-state problem. The final result is a n χ − µ n σ = 2 χa n χ 2 − (µ n σ ) 2 = 2 As discussed earlier, we are only interested in the case of ωτ << 1 for practical applications. Hence for an approximation up to 3iωτ = 3( χ − 1) , Eq. (A24) reduces to Eq.
(A6), which has been derived for the steady-state case. The roots of Eq. (A24) are then real and given by δ = ±δ j for j = 1,2,..., N − 1 . To find the two additional roots, note that Eq. (A24) suggests that there exist a constant µ n = µ 0 , such that (µ 0 δ ) 2 = χ(χ − 1) . Note that in this case, the second equality in Eq. (A22) becomes J (x,µ) −1 1 ∫ dµ = 2J (x,µ 0 ) , which is nothing more than the well-known mean value theorem of integrals, applied to the modulated intensity J. Therefore, the existence of the µ 0 is supported by this theorem, the remaining two solutions are therefore δµ 0 = ± χ(χ − 1) ≈ ± iωτ , and the general solution for J is where η = iωτ µ 0 and A, B, C j and C − j are arbitrary constants. By applying the boundary conditions of diffuse scattering for the modulated intensity J (see Eqs. (A9) and (A10)), it can be shown that The intensity J can now be expressed as The sum on each exponential term in Eq. (A29) represents the series expansion of a given function and therefore they can be written as follows Finally, the combination of Eqs. (A15) and (A28) yields the following modulated heat flux where M 1 is the negative integral of M defined in Eq. (38). Note that the derivative of q t is related to the equilibrium intensity in Eq. (A30) through the simple relation ′ q t (x) = −4π (χ − 1)J 0 (x) . This result can also be derived using the integral Eqs. (15a) and (17a), which shows the consistence of Eqs. (A30) and (A31).