Heat flux induced blueshift of dominant phonon wavelength and its impact on thermal conductivity

The concept of dominant phonon wavelength is investigated in systems submitted to a heat flux at low temperatures. Using spectral energy distributions, a treatment of two-dimensional and three-dimensional structures is conducted in parallel. We demonstrate a significant reduction of the dominant phonon wavelength, up to 62%, due to a displacement of the phonon spectrum towards higher frequencies in presence of a heat flux. We name this phenomenon blueshift effect. A formula is provided to directly calculate the corrected dominant phonon wavelength. We illustrate the impact of the blueshift effect by showing that a temperature gradient of 10% at 4K yields a 20% reduction in the thermal conductivity. Therefore, ignoring the blueshift effect in a thermal model can notably alter the physical interpretation of measurements. The results suggest that an appropriate heat flux environment can improve thermoelectric device performances.

The concept of dominant phonon wavelength is investigated in systems submitted to a heat flux at low temperatures. Using spectral energy distributions, a treatment of two-dimensional and three-dimensional structures is conducted in parallel. We demonstrate a significant reduction of the dominant phonon wavelength, up to 62%, due to a displacement of the phonon spectrum towards higher frequencies in presence of a heat flux. We name this phenomenon blueshift effect. A formula is provided to directly calculate the corrected dominant phonon wavelength. We illustrate the impact of the blueshift effect by showing that a temperature gradient of 10% at 4K yields a 20% reduction in the thermal conductivity. Therefore, ignoring the blueshift effect in a thermal model can notably alter the physical interpretation of measurements. The results suggest that an appropriate heat flux environment can improve thermoelectric device performances. © 2017 Author(s). All

I. INTRODUCTION
The necessity to control heat propagation in ceaselessly smaller devices has been sustaining the growing interest in nanoscale thermal transport for the past two decades. [1][2][3] Today, the electronic technologies are based on the Silicon in the frame of three-dimensional (3D) thermal physics. Meanwhile, two-dimensional (2D) materials are being intensively studied and exhibit remarkable heat transport properties. [4][5][6] In general, thermal conduction properties depend on temperature and phonons frequencies. At a given temperature, a broad phonon frequency spectrum, ranging from GHz to THz, 7,8 is excited, leading to phonon wavelengths ranging from µm to nm, respectively. Except by simulations, [9][10][11] it is not possible to experimentally track the specific contribution of each of these modes yet. However, it is common to extract a global behavior by using the notion of a dominant phonon angular frequency ω d defined as where refers to the reduced Planck constant, k B the Boltzmann constant, T the temperature and α a dimensionless factor that we name the dominant coefficient. ω d is often converted into the dominant phonon wavelength (DPW) λ d which is used at many occasions. [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28] The DPW represents a crucial and convenient parameter to immediately assess the physics of phonons without having to conduct sophisticated and time consuming simulations. It is strongly correlated to the heat propagation properties in a material. Also, the DPW can be compared to the system dimensions in order to evaluate the presence of coherence effects 29,30 and to determine the specularity parameter for the phonon-surface roughness interactions. 18 Despite its important role, several definitions of the DPW coexist in the literature differing by the dominant coefficient α. For 3D materials, four values of α can be found. First, 12-15 α = 4.25, which finds its origin in measurements conducted by Pohl. 31 A second commonly [16][17][18] applied coefficient is α = 2.82, which is often attributed to Ziman. 16 In the context of thermal contact resistances with superfluid Helium, [19][20][21][22] the peak of the specific heat is utilized to find α = 3.83. Furthermore, some authors 23-26 omit the factor α and simply set ω d ≈ k B T , arguing that an order of magnitude accuracy remains relevant, which may certainly be true in some cases. For 2D materials, fewer studies are available up to now, but the usual value is consistent with 27,28 α 2D ≈ 1.6 given by Ziman (Chap 8 Sec. 5 of Ref. 32).
Over time, the dominant coefficients α have been employed rather loosely, casting confusion on the value to apply to determine the DPW. In the more recent context of nanoscale thermal physics, both α = 4.25 and α = 2.82 are used, without any reasons justifying why one or the other should be chosen, although these values differ by 50%.
Until now, the DPW has been determined by using only one temperature, generally the regulated cold reservoir temperature of the system. However, in order to get a heat flow, a second heat source is necessary. Therefore, the effective heat flux flowing from the hot reservoir to the cold reservoir depends on the temperatures of both reservoirs. Yet, the effect of the heat flux on the DPW has never been considered before.
Our investigation reveals that for systems subjected to a heat flux, the phonon spectrum is shifted to higher frequencies. As a result, the DPW becomes significantly smaller as the temperature difference between the thermal reservoirs is decreased, a phenomenon we refer to as blueshift effect. In the limit of small temperature gradient, we show that not considering the heat flux leads to 62% error in the DPW of 2D materials. We detail the calculations describing the behavior of the DPW in 2D and 3D systems between two thermal reservoirs and provide the corrected dominant coefficient as function of the temperature difference. In the light of our results, we discuss the different DPWs available in literature and show the advantages to work with the phonon energy spectrum. Finally, the dependence of the DPW on temperature differences is shown to notably reduce the thermal conductivity. Thus, the blueshift effect alters the physical interpretation of experimental results and can be applied to improve the performances of thermoelectric devices.

A. Frequency dependent energy spectrum
The lattice energy associated with a thermal reservoir at a temperature T is given by the general where ω is the phonon angular frequency, f BE = (e ω/(k B T ) − 1) −1 represents the Bose-Einstein distribution and g m (ω) = ω 2 /(2π 2 c 3 m ) denotes the density of states per unit frequency for each phonon propagation mode m characterized at low temperature by a single group velocity c m = (dω/dk) m . The zero-point energy, represented by the 1/2 term in Eq. (2), does not play any role in the heat transfer because it is temperature independent. We conduct the calculations within the framework of the Debye theory for phonons so that ω max corresponds to the Debye angular frequency ω D and consider a linear dispersion relation ω = c m k.
The temperature dependent terms under the integral in Eq. (2) define the spectral radiance or spectral energy density per unit frequency: B m (ω, T ) = f BE (ω, T )g m (ω) ω. The difference between 2D and 3D systems lies in the density of states. The shape of the spectral energy density is given by B m (ω, T ) ∝ f BE ω 3 and B m (ω, T ) ∝ f BE ω 2 for 3D and 2D, respectively.
The dominant angular frequency ω d is associated to the peak of B m (ω, T ). The conversion to wavelengths yields the DPW λ d = 2πc m /ω d = hc m /(αk B T ), which is equivalent to Eq. (1). Letting Fig. 1(a). The peak of the spectral energy density is determined numerically with high accuracy to find the value α = 2.821 for 3D.
Along the same lines, we derive the 2D energy density spectrum is given Fig. 1(b). Compared to the 3D spectrum, the 2D spectrum is clearly more concentrated at the low frequencies. The dominant coefficient is then significantly reduced to be found for α = 1.594.

B. Specific heat distribution
Following the works by Pohl,12,31 the effective phonons at a given temperature are those that contribute predominantly to the specific heat. Hence, dominant phonon wavelengths can also be determined from the specific heat, which is the temperature derivative of the lattice energy defined in Eq. (2). Therefore, we can define a specific heat density spectrum by M(ω, T ) = ∂B(ω, T )/∂T . A similar normalization as for the energy leads to the normalized specific heat spectrumM(X) as shown in Fig. 1(c) and (d) for 3D and 2D cases, respectively.
The 3D specific heat spectrum ( Fig. 1(c)) is broader than the energy spectrum ( Fig. 1(a)) and peaks at higher frequencies with a dominant coefficient α = 3.830 while α = 2.576 is found for 2D.
The dominant coefficient we have calculated with the energy spectrum and the specific heat spectrum are identical to those given in the introduction and therefore explain the origin of the different DPWs found in the literature. In the case of 3D materials, the experimental value α = 4.25, while being a benchmark, seems to point in favor of the specific heat spectrum theory which gives α = 3.83. However, the energy spectrum theory, with α = 2.82, has also been often chosen. For 2D materials, in the absence of experimental data, the energy spectrum theory has been exclusively used.

III. DPW UNDER HEAT FLUX
The calculations we made in the previous section considered that the thermal properties in a system, 3D or 2D, are given by only one phonon reservoir. However, a heat flow exists only in presence of a temperature difference. We represent this situation by surrounding the system by two reservoirs (see inset Fig.2) and observe the impact of the heat flux on the DPW.

A. Heat flux formalism
First, we consider a 3D thermal reservoir at temperature T. It emits an energy density flux given by: For a medium subjected to a temperature gradient due to a hot reservoir at temperature T h and a cold reservoir at T c (with T c < T h ), the effective energy density flux flow through the medium is where the hot and cold reservoirs emit a phonon density flux q h and q c , respectively. Taking the heat conduction in the direction of the temperature gradient, the orthogonal projection of the velocities on the x axis leads to an effective heat flux given by: where B eff m (ω, T h , T c ) is the effective energy distribution per mode and is defined by: From Eq. (5) and Eq. (6), it is immediately visible that for T h = T c there is no heat flux in the system. This situation is similar to the one we describe in Sec. II. As a supportive example, we represent in Fig. 2 the different terms in Eq. (6) with T h = 4.4K, T c = 4.0K and by setting c m = 6000 m.s −1 . We clearly observe the displacement of the peak frequency of the effective energy distribution B eff m (ω, T h , T c ) to a higher frequency compared to that of the hot B m (ω, T h ) and cold B m (ω, T c ) energy distributions. Therefore, phonons contributing predominantly to the effective heat transport have higher frequencies compared to phonons contributing to the peak frequencies of the individual thermal reservoirs. We refer to this shift to higher frequencies as blueshift effect. This effect is all the more substantial for smaller heat fluxes, when T h approaches T c . Since the cold energy distribution is always overlapped by the hot one, the impact of the cold energy distribution is gradually suppressed as the temperature difference increases between the two  reservoirs. Consequently, the blueshift effect is damped until saturation to a constant value for very large temperature differences, as shown later.
To quantify the blueshift effect, we determine the peak of the effective energy distribution Eq. (6) by solving the equation (see details in appendix) whereθ = ω/k B is a reduced temperature.
The blueshift effect can also be observed by considering the specific heat. As in Eq. (5), an effective global specific heat in the system is defined as: where M eff m (ω, T h , T c ) is the specific heat density spectrum per unit frequency and per mode given by: A similar expression as Eq. (7) can then be established by using Eq. (9). Eq. (6) and Eq. (8) are valid for any dimensionality of the system; only the phonon density of states has to be changed accordingly. Therefore, the concepts we have developed up to now are also applicable to 2D materials.
Actually, we generalized Eq. (7) to encompass the four cases that are the 2D and 3D energy peaks and the 2D and 3D specific heat peaks by taking into account the dimensionality δ of the system that is considered and the order of derivative n with report to the energy distribution. The general equation is given by (see Appendix):  Fig. 3(a) shows the spectrum of α =θ d /T h values on a plot of the temperatures of the hot and cold reservoirs, after solving Eq. (10) for n = 0 and δ = 3, that is, for a 3D effective energy density distribution. The color scale depicts the evolution of α as a function of T h and T c . The linear contour lines passing through the origin displays equal values of α for different T h /T c ratios. On the diagonal, T h → T c (q eff → 0 + ) and α has a maximum value of α 0 = 3.830. As ∆T = T h − T c increases and tends to T h T c , α decreases to attain a constant value of α ∞ = 2.821. In this case the effective heat flux is at its maximum so that q c is negligible compared to q h and the situation with a single thermal reservoir is again retrieved. Fig. 3(b) showsθ d as a function of T h for three different T h /T c ratios. Calculations were done with n = 0 in Eq. (10); and the δ = 3 and δ = 2 cases are treated for comparison. All datasets indicate a perfect linear dependence betweenθ d and T h for a constant T h /T c ratio. The slope of each curve in Fig. 3(b) represents the dominant coefficient α for a given T h /T c ratio. Clearly, α decreases and stabilizes as the temperature differences increase. Fig. 4 summarizes the behavior of α as a function of for 2D and 3D systems respectively. For each of these cases, the calculations are performed with the effective energy density distribution and with the effective specific heat density distribution for comparison. All the curves in Fig. 4 have the same exponentially decreasing trend, which we express as:

B. DPW dependence to ∆T
where α 0 corresponds α values when ∆T → 0 + and α ∞ represents either ∆T = 0 or ∆T → ∞ because these two cases depict the single thermal reservoir situation as mentioned earlier. The γ parameter represents a characteristic decay constant of α from α 0 to α ∞ as ∆T/T c increases. From Table I, the γ values for the 3D case are larger than the γ values for the 2D case. Consequently, in 3D systems, the evolution of α is completely damped at ∆T/T c ≈ 2, that is T h → 3T c . And, in 2D systems evolution of α persists up to ∆T/T c ≈ 3, that is T h → 5T c . All values of α 0 and α ∞ are displayed in Table I. It is interesting to note that the α ∞ values are identical to the α values given in Section II. These results imply that for temperature differences larger than those leading to a constant α value (plateau region in Fig. 4), the contribution of the cold thermal reservoir to the net heat flux can be neglected. Finally, for large temperature differences, we retrieve the configuration of a single reservoir, as discussed in the first part of this paper. On the other hand, in the limit of small temperature differences, we remark that α 0 for B eff (ω, T h , T c ) is identical to α ∞ for M eff (ω, T h , T c ) in Table I. This result is logical since the definition of the specific heat is the temperature derivative of the energy. Indeed, a simple Taylor development as ∆T → 0 + yields B eff = M eff ∆T . Therefore, the peak of B eff (ω, T h , T c ) when ∆T = 0 + corresponds to the peak of M eff (ω, T h , T c ) when ∆T 1. From Eq. (1) and (11), we finally define the general The relative difference between α 0 and α ∞ compared to α ∞ is ∼62% in the 2D case and ∼36% in 3D case. These percentages represent the maximum error that can be committed by ignoring the effective heat flux. The increase of α demonstrated Eq. (11) has a direct consequence to lower λ d , and therefore to induce a frequency blueshift. This blueshift is all the more important in thermal experiments because the relative temperature difference is generally maintained under 10% in order to avoid non-linear effects. 34 This region of small ∆T/T c is precisely where the blueshift is the strongest because of the exponential decay form. For example, if we consider the energy spectrum in Fig. 4(b) for the 3D case, we see that a ∆T /T c = 0.1 leads to α ≈ 3.66 which is ∼30% above the widely used value of 2.82 corresponding to the peak of the energy density distribution.

C. Comment on the use of the energy spectrum
The only experimental value of the DPW was given by Pohl who deduced α = 4.25 thanks to low temperature thermal conductivity measurements. 31 Taking into account the 10% experimental error, this value was then attributed to the peak of the specific heat spectrum 12 which theoretically yields α = 3.83.
Thermal conductivity measurements requires the presence of a heat flux in the material. The experimental value of α can therefore be also interpreted as a manifestation of the blueshift effect. Indeed, in the small heat flux limit, we have demonstrated that the energy spectrum also leads to α = 3.83 which is in agreement with the measurements. The parallel between phonon and photon makes possible their common physical treatment as part of the boson family. As for photons, only the energy spectrum is used. 35,36 Finally, as the blueshift effect reconciles the experiment and the energy theoretical value, we would recommend to use the energy spectrum to determine the DPW.

D. Blueshift effect on the thermal conductivity
In the boundary scattering regime, the mean free path of phonons is strongly affected by phononsurface roughness interactions. 11 The nature of these interactions is determined by a specularity factor which represents the fraction of total number of phonons undergoing specular scattering; and it is given, according to Ziman's formula, 32 by p(σ, λ) = exp(−16π 2 σ 2 /λ 2 ), where σ is the surface roughness height. A surface for which σ λ is defined as perfectly rough and phonon scattering is fully diffusive, that is p(σ, λ) = 0. Casimir 37 defines the thermal conductivity under these conditions to be written as: κ diff = C p cΛ diff /3, where Λ diff is the mean free path of the diffusive phonons.
When both specular and diffuse scattering are present, the effective mean free path is given as wherep(λ d ) is the average specularity parameter. Ziman showed that it depends on λ d following the integral 32p The distribution of roughnesses P(σ) we use in Eq. 14 was determined experimentally by Heron et al 18 as P(σ) = e −σ/σ 0 where σ 0 is the root mean square surface roughness. Inserting Eq. (12) into Eq. (14), the blueshift effect decreases the average specularity parameter and so reduce the effective mean free path. The impact of the blueshift effect on the thermal conductivity is assessed by using Eq. (13) which directly leads to the κ eff /κ diff ratio. Fig. 5 shows this ratio as a function of temperature for a 3D structure with a root mean square roughness σ 0 = 4nm. For a given cold temperature T c , the κ eff /κ diff ratio has a maximum value when the system is at thermal equilibrium ∆T/T c = 0 (which also corresponds to the single thermal reservoir case). This situation is represented by the dashed line. The lower limit, pictured by the full line, is obtained when the blueshift is maximum which is attained FIG. 5. κ eff /κ diff calculated by using the DPW for the 3D energy distribution and a roughness σ 0 = 4nm as a function of temperature for different ∆T/T c . The inset shows the relative difference of the R = κ eff /κ diff ratio at a given ∆T/T c , compared to the R 0 ratio at thermal equilibrium.

015017-9
Ramiere, Volz, and Amrit AIP Advances 7, 015017 (2017) when the temperature difference tends to zero ∆T/T c → 0 + , without being strictly equal to zero. Among all the intermediate relative temperature differences, we plot the special case ∆T/T c = 10%, which corresponds to the conventional upper limit to avoid non-linear heat flux effects, with a blue dotted line. The results clearly indicate a decreasing effective thermal conductivity due to a reduction of the DPW caused by the blueshift effect. The inset in Fig. 5 shows the relative difference in percentage between the R = κ eff /κ diff ratio at a given ∆T /T c and R 0 = κ eff /κ diff for ∆T /T c = 0 at thermal equilibrium. For ∆T /T c = 10%, the observed reduction exceeds 20% below 4K showing the importance of the blueshift effect at low temperature. As the temperature increases the relative difference decreases because of the transition from the semi-ballistic regime to the diffusive regime. Consequently, the impact of the blueshift is below 5% above 20K.
High performance thermoelectric modules are important for cooling electronic circuits in quantum computing and superconducting devices. 38 The performances of thermoelectric devices can be improved by two ways thanks to the blueshift effect by choosing the correct thermal environment. Firstly, thermoelectric devices working according to incoherent thermal transport principles, e.g. rough nanowires, 39 should be placed under low heat flux conditions. Following the above discussion, a smaller temperature difference yields a smaller thermal conductivity, thereby increasing the thermoelectric figure of merit. Secondly, phononic crystals, where partially coherent thermal transport occurs, should be placed under high heat flux. Indeed, it has been recently shown that the periodic pattern of phononic crystals creates interference of coherent phonons by Bragg diffraction, thus hindering heat propagation. 40 The higher the heat flux, the higher the DPW which leads to a higher proportion of coherent phonons finally providing a lower thermal conductivity. However, increasing the heat flux in a structure has also an impact which decreases the DPW because of the temperature increase. The temperature difference in the structure, while relatively large, has therefore to remain limited. As reported in Fig. 4, maximum performances of thermoelectric phononic crystals should be attained for ∆T /T c ∼ 3 in 2D materials and ∼ 2 in 3D materials. Finally, we note that in Ref. 11 we give the relationship between sample size and temperature for boundary scattering to predominate.

IV. CONCLUSION
Having established and consolidated the definitions of DPW, we then studied its evolution induced by a heat flux imposed by two thermal reservoirs. We show that the effective energy and specific heat density distributions describing the system undergo a shift to higher frequencies that is controlled by the temperatures of the two thermal reservoirs. Consequently, the blueshift effect now leads to the dominant coefficient α to vary as a function of the temperature difference. We have established that the DPW is reduced as the heat flux decreases, following an exponential decay law. This reduction is substantial as it can reach up to ∼62% for 2D systems and ∼36% for 3D systems.
The impact of the DPW shift on the thermal conductivity has been studied via the phonon mean free path which depends on the average specularity parameterp. The results clearly show that the relative change in the thermal conductivity κ eff /κ diff decreases by more than 20% at 4K for temperature differences as small as 10%. It indicates that an adequate heat flux environment can improve the performances of thermoelectric devices. Besides, the misjudgment on the thermal conductivity by ignoring the blueshift effect can significantly alter the physical interpretation of experimental results.
Our calculations provide a ready-to-use formula to determine the dominant phonon wavelength. It serves as a useful and immediately available reference especially for nanostructured materials where effects such as confinement or thermal rectification may operate. Experimental measurements of the decay constant γ in the expression for α may help to corroborate our predictions of the blueshift effect in presence of a heat flux.

ACKNOWLEDGMENTS
This work received support of the French Ministry of Education via ED 543 MIPEGE of the University of Paris-Sud.

APPENDIX
This appendix presents the details of the derivation of the DPW expressions Eq. (7) in the 3D energetic case and its generalization to 2D and 3D with the specific heat given Eq. (10).
The 3D spectral density of energy in the Debye approximation can be written as which yields the effective spectral density of energy By definition of the peak, Eq. (A.4) is equal to zero whenθ =θ d , so We separate the terms with the hot temperature T h from those with T c to find a symmetrical form identical to Eq. Then the derivative yields It is now straightforward to write the equation to find the maximum
The complete generalization Eq. (10) is made by noticing that the 3D spectral density of energy expression Eq. (A.1) can be extended to incorporate the 2D case by introducing the dimensionality parameter δ as follow The same derivation line is then implemented for 2D systems.