Research ArticlePHYSICS

Observation of second sound in a rapidly varying temperature field in Ge

See allHide authors and affiliations

Science Advances  30 Jun 2021:
Vol. 7, no. 27, eabg4677
DOI: 10.1126/sciadv.abg4677


Second sound is known as the thermal transport regime where heat is carried by temperature waves. Its experimental observation was previously restricted to a small number of materials, usually in rather narrow temperature windows. We show that it is possible to overcome these limitations by driving the system with a rapidly varying temperature field. High-frequency second sound is demonstrated in bulk natural Ge between 7 K and room temperature by studying the phase lag of the thermal response under a harmonic high-frequency external thermal excitation and addressing the relaxation time and the propagation velocity of the heat waves. These results provide a route to investigate the potential of wave-like heat transport in almost any material, opening opportunities to control heat through its oscillatory nature.


The study of heat transport beyond Fourier’s regime has attracted renewed interest in recent years. Great efforts have been performed to unravel the physical properties of thermal waves, as well as the experimental conditions that are necessary for their observation. Applications based on such concepts have been envisioned and discussed extensively already in many recent publications (14). The spatiotemporal propagation of the temperature field in the form of waves is known as “second sound,” a term that was adopted in analogy to “first sound” (or simply “sound,” i.e., mechanical lattice vibrations). As pointed out in (5), first sound and second sound are described by a similar equation where the variables have a different physical meaning, i.e., pressure and temperature, respectively.

The simplest differential equation that describes wave-like heat transport from a mesoscopic perspective is the hyperbolic heat equation (HHE) according to Maxwell, Cattaneo, and Vernotteτss2Tt2+Ttα2T=1ρCp(S(r,t)+τssS(r,t)t)(1)where α is the thermal diffusivity, τss is the thermal relaxation time, ρ is the mass density, Cp is the specific heat, and S(r,t) is an external power heat source. In our context, the system is in local equilibrium, and it is well characterized by a local temperature T (see discussion in section S8). The previous equation describes the propagation of a temperature wave with a damping term given by ∂T/∂t and a propagation velocity υss=α/τss. The solutions of this equation lead to different heat transport regimes, depending on the temporal and spatial length scales under investigation. The key to unlock the different regimes is the magnitude of the first term on the left-hand side of Eq. 1, thermal inertial term; i.e., if τss or ∂2T/∂t2 is sufficiently large, the spatiotemporal distribution of temperature field will exhibit wave-like behavior.

Second sound in solids was first experimentally observed in solid He (6); later in NaF (7), Bi (8), and SrTiO3 (9); and most recently in highly oriented pyrolytic graphite (10). Several theoretical works have also recently addressed its occurrence in low-dimensional systems (1113). In all these experimental observations of second sound, the dominance of momentum conserving phonon scattering (Normal processes) with respect to resistive phonon scattering (Umklapp processes) was found to be the key mechanism leading to its observation. Second sound was observed almost exclusively in the very low temperature regime (T < 5 K), with the exception of a recent example (10) at higher temperatures (125 K) for samples with low resistive phonon scattering.

A condition for the experimental detection of second sound, based on these experimental observations (610), was found to be τN < τexp < τR; i.e., the typical experimental observation times (τexp) must be larger than normal phonon scattering times (τN) to allow momentum redistribution but smaller than resistive phonon scattering times (τR) to avoid decay of the phonon wave packet into the phonon equilibrium distribution. The theoretical foundations of second sound were set in the 1960s by M. Chester, R. J. Hardy, C. P. Enz, and co-workers (5, 14, 15), who predicted the existence of “drifting” second sound, i.e., a type of wave-like heat transport that is triggered by the dominance of normal phonon scattering events. This type of second sound was experimentally confirmed in (610). Furthermore, the existence of a “driftless” or “high-frequency” type, as well as “other types” of second sound, was also envisioned. In this case, the dominance of normal scattering events is not a necessary condition for the existence of wave-like heat transport, but the key requirement is instead the slow decay of the energy flux, as predicted by different solutions to the Boltzmann transport equation (5). These different flavors of second sound are all contemplated in the mesoscopic HHE (Eq. 1), because different microscopic mechanisms can lead to the propagation of thermal waves.

We show that it is possible to observe a type of high-frequency second sound in natural bulk Ge by driving the system out of equilibrium with a rapidly varying temperature field. Our concept is based on taking advantage of the second-order time derivative in the HHE, Eq. 1, in a frequency-domain experiment. As the driving frequency increases toward the hundreds of megahertz range, the relative weight of this term with respect to the damping term (first-order time derivative) increases proportionally to the frequency upon a harmonic excitation, hence making the observation of wave-like heat propagation possible. We show that this approach is robust enough to expose second sound independently, to a certain extent, of the phonon scattering rates of the studied material, as well as of temperature. Although heat transport in Ge is dominated by resistive phonon scattering processes, which partly originate from its large isotopic diversity, we show that it is still possible to observe second sound in the high-frequency limit.


Our experiments are based on a frequency-domain optical reflectance pump-and-probe approach based on two lasers with different wavelengths (λpump = 405 nm and λprobe = 532 nm) focused onto the surface of a Ge sample to a spot size with radius Rspot ≈ 5.5 μm. The studied samples are pieces of a substrate of natural Ge. Further details are provided in the Supplementary Materials (section S1). The pump laser (thermal excitation) is modulated between 30 kHz and 200 MHz with a sinusoidal power output waveform, leading to a dynamic modulation of the optical reflectivity of the surface of the sample, which is also well described by a harmonic waveform. A frequency-dependent phase lag gradually develops, defined as a phase difference between the harmonic thermal excitation,  S(r,t), and the response of the sample,  T(r,t) (sections S2 and S7), which can be modeled using Eq. 1 in the present experimental conditions. The choice of Ge as a candidate for the observation of second sound is not arbitrary, and it is mostly based on the large optical absorption coefficient of this material for the wavelengths used in this experiment. The optical penetration depth of the pump and probe lasers is δ405nm = 15 nm and δ532nm = 17 nm, respectively (section S3). These particular conditions make Ge an ideal material for this study, because the small penetration depth of both lasers ensures that the measured phase lag is local and, thus, accurately describes the oscillations of the thermal waves. The penetration depth of the pump and probe lasers should be compared with the thermal wavelength υss/f. At room temperature and the highest excitation frequency considered, υss/f = 840 nm, thus considerably larger than the optical penetration depth. We note that wave-like effects are not observed in the presence of a metallic transducer (section S9), because the thermal interface between the transducer and the Ge substrate, as well as the transducer itself, dominates the system response in the frequency range where wave-like effects are expected.

In the present experiments, the changes ΔR of the optical reflectivity upon the pump laser excitation are driven by the lattice temperature variation, hereafter referred to as ΔT, and by the modulation Δn of the carrier concentration in the conduction band. The resulting ΔR at the energy corresponding to the probe wavelength is provided by the first-order expansion ΔR = (∂R/∂TT + (∂R/∂nn (16). Interestingly enough, the contribution provided by the variation of the carrier concentration is expected to dominate the way the reflectivity is affected in experiments, like the present ones, involving pulsed laser sources, i.e., for high electronic excitation densities (17, 18). However, in our excitation conditions, the electronic contribution to the optical reflectivity can be neglected; hence, the optical reflectivity is dominated by the temperature of the lattice (ΔT) for all excitation frequencies. We demonstrate this by studying the relative magnitude of (∂R/∂TT and (∂R/∂nn. In particular, we take advantage of the fact that for bulk Ge, ∂R/∂T ≈ 0 at T = 220 K for 532 nm of probe wavelength (see section S3), whereas ∂R/∂n is expected to be almost independent of temperature between 220 K and room temperature. Furthermore, ∂R/∂T exhibits a sign inversion at this temperature, which is reflected in a change of the phase lag by an angle of π. Figure 1A displays the temperature dependence of ΔR/R, as well as the phase lag for low (30 kHz) and high (100 MHz) modulation frequencies, and for a constant pump power of ≈10 mW. The reflectivity signal exhibits a minimum (19) around ≈220 K, followed by a slow signal recovery at lower temperatures. In the lower temperature range (T < 100 K), the lattice and electronic contributions to the phase lag have opposite signs. Hence, the measured phase lag at high frequency in this range cannot be explained with the electronic contribution to the optical reflectivity. On the other hand, the phase lag is shifted by π at ≈220 K (140°→−40°). The optical reflectivity extracted from ellipsometry experiments corresponding to ΔT = 1 K is also shown for relative comparison. A similar behavior is observed independently of the excitation frequency, which resembles the ellipsometry data for which Δn = 0. From the reflectivity signal ratio between the minimum observed in Fig. 1A and the value of the reflectivity at room temperature, we can estimate an upper boundary for (1/R)[∂R/∂n] ≈ −2 × 10−29 m3, assuming that the residual signal at ≈220 K has a purely electronic origin. The inset of Fig. 1A displays calculations of the electronic (solving the electron recombination diffusion equation) and lattice (solving HHE) contributions to the optical reflectivity as a function of frequency (section S3). Whereas the lattice contribution to the optical reflectivity is at least 20-fold larger than the electronic component at the higher frequencies, for lower frequencies, this ratio is substantially increased.

Fig. 1 Optical reflectivity, experimental and simulated phase lag versus frequency, and thermal penetration depth for a Ge sample.

(A) The top panel displays the optical reflectivity change as a function of temperature for 30 kHz (red) and 100 MHz (blue). The black dots were obtained from ellipsometry measurements to a temperature rise of 1 K (full symbols→dR/dT < 0, open symbols→dR/dT > 0). The bottom panel displays the phase lag of the signal with respect to the pump excitation. The phase lag is directly obtained from the measurements and must be corrected by ±180° to normalize it to the [0;±45°] as in (B). The inset displays calculations accounting for the thermal and electronic contributions to the reflectivity at 300 K. (B) The experimental phase lag as a function of the pump excitation frequency is shown in black open symbols. The inset displays numerical experiments using NEMD. In dashed red line, we display the prediction based on Fourier’s law. The solutions based on the 3D HHE are shown in a black line with a resulting fitted τss = 500 ps. (C) Schematic illustration of the geometry used for the NEMD numerical experiments. (D) Frequency-dependent thermal penetration depth calculated using the solution of the HHE (ʌHHE), the diffusive case (ʌdiff), and the penetration depth (ʌss) obtained in the high-frequency limit.

Figure 1B displays the experimental phase lag as a function of frequency between 30 kHz and 200 MHz at room temperature. The complex thermal response of the specimen was, at first, computed numerically within Fourier’s model, solving the parabolic approximation to the three-dimensional (3D) HHE (diffusive case), which is obtained when the first term of Eq. 1 can be neglected. In the range between 30 kHz and 1 MHz, the agreement between Fourier’s solution and the experimental data is excellent, although deviations are already observed around 1 MHz. Above 30 MHz, the difference between the experimental phase lag and Fourier’s predictions is evident. For the higher-frequency range, the experimental data show that the phase lag (absolute value) decreases with increasing frequency. This trend cannot even be qualitatively reproduced by Fourier’s model, which predicts that as frequency increases, the phase lag approaches −π/4 and even lower values (see section S7 and fig. S7-3). The full 3D solution of the HHE based on the finite-element method (FEM) was used to fit the experimental data through the entire frequency range, and it is shown in Fig. 1B. A detailed description of the fitting procedure is presented in section S7. We obtained τssexp = 500 ps and αexp = 3 × 10−5 m2/s, thus leading to a propagation velocity υssexp=250 m/s, all at room temperature. These experimental observations are in qualitative agreement with the computational experiments (see schematic illustration in Fig. 1C) by nonequilibrium molecular dynamics (NEMD) and are shown in the inset of Fig. 1B (see also section S5 and fig. S5-1). Although a quantitative agreement cannot be expected because of differences between the experimental and the computational setups (reduced size of the sample and purely 1D heat flux in NEMD), the NEMD results show a notable similarity with the experimental results, with the phase lag initially decreasing, hitting a minimum, and then recovering. These results are of particular interest because, within the NEMD approach, no assumption is made regarding the heat transport regime. Thus, the numerical experiments are an independent confirmation of the appearance of second sound in the high-frequency limit. The observed large reduction of the phase lag with respect to the Fourier prediction (see Fig. 1B) indicates the emergence of a more efficient heat transport mechanism as compared to diffusion at high frequencies. Ballistic behavior of phonons (2022) cannot be invoked to rationalize the experimentally observed dependence of the phase lag with frequency. As discussed previously in (21), a reduction of the thermal conductivity at high frequency would lead to an increase in the phase lag (see fig. S7-4), rather than a decrease (absolute value), as we observe. Conversely, the observed frequency-dependent phase lag is well reproduced by the HHE, i.e., considering the contribution of wave-like heat transport.

The frequency window where second sound is expected can be estimated comparing the thermal penetration depth of the diffusive and wave-like regimes. Figure 1D displays the penetration depth, ΛHHE, calculated using the exact solution of the HHE (see section S7), as well as the diffusive and wave-like limits, Λdiff=α/(πf) and Λss=2ατss, respectively. A critical frequency fc is obtained when Λdiff = Λss, thus providing an estimation of the frequency for which the diffusive and wave-like contributions to heat transport are similar. The temperature dependence of τss, υss, and fc was studied between room temperature and 7 K, and it is shown in Fig. 2A (full set in section S2). As temperature decreases, the ratio between the penetration depth of the wave-like and diffusive contributions is Λss/Λdiff=4πfτss(T), which implies that lower temperatures favor the spatial propagation of the thermal waves because larger τss is expected and indeed experimentally observed for lower temperatures. Wave-like effects are already present below fc, as can be observed comparing the experimental data with the corresponding fits using the HHE to the Fourier predictions, as shown in Fig. 2A. The onset frequency can be estimated as f > 0.1fc, corresponding to a phase lag difference between the experimental data and the Fourier solution >5°. The onset of wave-like effects is also evidenced by the deviations between Λdiff and ΛHHE as shown in Fig. 1D. The minimum observed on the phase lag curves and the position of the critical frequency fc relative to the frequency of the minimum originate from the relation between a, τss, and Rspot (see discussion in section S7). The frequency-dependent phase lag in Fig. 2A was fitted at each temperature using the HHE, as described for the room temperature case, and the results for τssexp and υssexp are shown in Fig. 2B and table S2. We note that the accuracy of the fits reduces for lower temperatures due to the frequency-dependent temperature rise induced by the pump laser.

Fig. 2 Phase lag curves versus temperature, relaxation time, and velocity of second sound (theory and experiment), and simulated spatial propagation profiles for the diffusive and wave-like cases.

(A) Phase lag versus frequency for the higher-frequency range as a function of temperature. We plot three points at 300, 100, and 15 K with the corresponding fits to the data point using the 3D HHE. In dashed lines, we display the prediction based on Fourier’s law at each temperature. (B) The experimental relaxation times (τss), as well as the propagation velocity (vss), are shown as a function of temperature. The dashed lines are guides to the eye. The full lines are the predictions based on the expansion of the perturbed phonon distribution function, f=fλeq+βλq+γλ(q/t), combined with DFT simulations as described in sections S6 and S4, respectively. (C) Finite-element simulations of the spatial distribution of the temperature field at a function of temperature in the direction perpendicular to the surface of the sample at the highest excitation frequency of ≈300 MHz for an arbitrary time. The parabolic and hyperbolic solutions are shown in dashed and full lines, respectively.

To understand the origin of these observations, we have developed a rather simple model (see derivation in section S6) based on the expansion of the perturbed phonon distribution function (23), i.e., the intermediate state assumed by the nonequilibrium distribution before it decays to the equilibrium one by means of dissipative resistive processes, as: fλ=fλeq+βλq+γλ(q/t), where fλeq is the equilibrium phonon distribution function, q is the heat flux, βλ and γλ are mode-dependent functions to be determined, t is the temporal coordinate, and λ denotes each phonon mode. We note that the case with γλ=0 leads to the same propagation velocity for second sound as in (5) (section S6). In our case, however, the presence of rapidly varying temperatures leads to rapidly varying q, suggesting that the expansion of the perturbed phonon distribution function in terms of q/t is a reasonable assumption. The previous ansatz for f was then introduced into the linearized Boltzmann transport equation (BTE) to find the solution for βλ and γλ. It can be shown that within this framework, it is possible to derive the HHE, with a corresponding explicit expression for τss, in terms of individual phonon relaxation times (τλ) as shown in section S6τss=λωλνλ2τλ2fλeqTλωλνλ2τλfλeqT(2)where ℏ is the Planck constant, ωλ is the phonon energy, and νλ is the phonon group velocity. We computed ωλ, τλ, and νλ from the solution of the BTE based on density functional theory (DFT) interatomic force constants (IFCs) (section S4 and table S4). We restricted ourselves to the relaxation time approximation (RTA) after verifying that the full iterative BTE picture does not alter the prediction of the theory. Using the RTA has the additional benefit of providing a comparison on equal footing with previous theoretical descriptions (5) and allowing unambiguous definition of the relaxation times (24). We observe that the corrections to the RTA provided by a full iterative solution of the BTE are very small (within 4% in the thermal conductivity) for Ge at temperatures as low as 50 K (fig. S4-1). The resulting values were inserted into Eq. 2, which yielded τsstheo and υsstheo=αtheo/τsstheo as a function of temperature, as shown in Fig. 2B. The agreement of the predicted values with those obtained from the experiments is remarkable for T > 100 K, considering that the values of ωλ, τλ, and νλ are evaluated within a fully ab initio scheme. We note that the model leading to Eq. 2 is expected to be valid for ∣ΔT∣<< T, where ∣ΔT∣ is the amplitude of the laser-induced thermal oscillations and T is the absolute temperature as set by the cryostat. In our experiments, ∣ΔTmax < 20 K; thus, the observed deviations between the theoretical predictions and the measured values at very low temperatures are expected (see section S8 for details). We note that similar experiments (25, 26) were explained in terms of nonlocality as described by the Guyer-Krumhansl equation, which reduces to Eq. 1 in the absence of nonlocal effects. The present experiments, however, are beyond the applicability of this equation because the heating region is much smaller than the nonlocal length (27). Future work should address the attenuation of nonlocal effects under the present conditions and its influence in the propagation of thermal waves.

The spatial dependence of the temperature field in the parabolic (diffusive) and the hyperbolic (wave-like) cases was simulated using FEMs in the direction perpendicular to the surface of the sample at an arbitrary time. Figure 2C displays the normalized temperature profiles for 15, 100, and 300 K at the highest experimental excitation frequency of ≈300 MHz (see fig. S7-5 for simulations of the temperature field at fc). As expected, the wave-like behavior of the temperature field exhibits a strong temperature dependence. The observed propagation depth, particularly at lower temperatures, is especially interesting if considering the possibility of high-frequency modulated thermal interference. We think that the present approach could open interesting possibilities for the experimental observation of wave-like heat transport in other materials and lead to the development of strategies to control heat transport.



All samples were pieces cleaved from a 2-inch-diameter nominally undoped Ge wafer with (100) crystallographic orientation, high resistivity (>40 ohm·cm), and etch pit density <3000/cm2. The wafer was purchased from International Wafer Service Inc. (USA). We have studied three different types of samples: (i) clean Ge (no native oxide), (ii) Ge + native oxide + 60 nm of Au, and (iii) a 60-nm-thick Au transducer was evaporated onto the surface of a Ge piece similar to (i). Additional details on the etching and evaporation processes of the samples are given in section S1.

Amplified frequency-domain thermoreflectance

We developed a low-noise custom-built frequency-domain thermoreflectance setup to measure the thermal response of the samples. The thermal excitation was provided by a pump laser diode with a wavelength of 405 nm. The probe laser used was a continuous-wave (CW) laser of 532-nm wavelength and with a maximum output power of 100 mW. The output power was controlled with neutral density filters to ≈50 μW (CW) for the probe and ≈20 mW (RMS) for the pump laser. A 30-mm achromatic lens was used to focus both Gaussian beams onto the same spot, whose size was measured using the knife-edge method to a 1/e2 radius of ≈5.5 μm. A portion of the probe reflected light was sent to an avalanche photodiode detector. Two notch filters were inserted into the optical path with the purpose of individually blocking the laser components. First, the probe component was blocked, and the phase of the pump laser was measured using a high-frequency lock-in amplifier. After performing this phase calibration step, the notch filter mount was mechanically displaced to block the pump laser component, thus allowing us to measure the harmonic signal arising from the probe laser. An extended description of the experimental setup is given in section S2.

Finite-element modeling

The 3D electron diffusion-recombination equation is solved using FEMs with COMSOL Multiphysics to estimate the electronic contribution to the photoreflectance signal (see section S3C). The 3D HHE is also solved using FEM to calculate the phase lag between the harmonic laser excitation and the temperature response of the system (see section S7). The laser energy deposition is restricted to a region defined by the Gaussian function of the pump beam in the radial direction and an exponential decay in the cross-plane direction with the characteristic length of the optical penetration. The temperature oscillations correspond to a weighted average across the surface of the sample computed using the Gaussian function of the probe beam as the weight.

Ab initio calculations

The thermal conductivity and the heat flux relaxation time expressions in terms of the velocities and the relaxation times of the phonon modes are obtained by solving the BTE (see derivation in section S6). The required microscopic properties are calculated within DFT using VASP code (28). The harmonic IFCs were calculated from finite differences in a 5 × 5 × 5 supercell, while we used a 4 × 4 × 4 supercell for the anharmonic ones, limiting the interactions to fourth neighbors. The inequivalent displacements needed to obtain the IFCs were obtained with the phonopy (29) and (24) codes (more details can be found in section S4).

Nonequilibrium molecular dynamics

Computational experiments within NEMD used a 5 × 5 × 145 supercell of the eight-atom cubic Ge conventional cell. The interatomic interactions were described by a bond-order potential of the Tersoff type (30). All the NEMD simulations were performed using the LAMMPS code (31). The equations of motion were integrated with a time step of 1 fs, and temperature control was obtained by Nosé-Hoover thermostatting, while equations of motions have been integrated by the velocity-Verlet algorithm (more details can be found in section S5).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge M. Campoy-Quiles and Sergey Sobolev for fruitful scientific discussions and a critical reading of the manuscript and K. Esfarjani for discussions during the development of the BTE model. Funding: We acknowledge financial support from the Spanish Ministry of Economy, Industry, and Competitiveness through the “Severo Ochoa” Program for Centers of Excellence in R&D (SEV-2015-0496 and CEX2019-000917-S), MAT2017-90024-P (TANGENTS)-EI/FEDER, 2020AEP141, grant no. RTI2018-097876-B-C22 (MCIU/AEI/FEDER, UE), and the Generalitat de Catalunya under grant nos. 2017-SGR-1506 and 2017-SGR-00488. We thank the Centro de Supercomputación de Galicia (CESGA) for the use of their computational resources. M.L.-S. was funded through the Juan de la Cierva programme. Author contributions: This work was conceived and led by J.S.R. The experiments were done by L.A.P., M.I.A., and J.S.R. NEMD simulations were conceptualized, executed, and interpreted by M.L.-S., C.M., and L.C. Finite-element modeling was done by A.B. and F.X.A. Ab initio simulations were carried out by R.R. Modeling based on BTE was done by A.B., L.S., J.C., J.B, and F.X.A. All authors contributed to drafting the manuscript and participated in the scientific discussion. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article