## Abstract

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.

## INTRODUCTION

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 (*1*–*4*). 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_{ss} is the thermal relaxation time, ρ is the mass density, *C _{p}* is the specific heat, and

*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}or ∂

^{2}

*T*/∂t

^{2}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 SrTiO_{3} (*9*); and most recently in highly oriented pyrolytic graphite (*10*). Several theoretical works have also recently addressed its occurrence in low-dimensional systems (*11*–*13*). 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 (*6*–*10*), was found to be τ* _{N}* < τ

_{exp}< τ

*; i.e., the typical experimental observation times (τ*

_{R}_{exp}) must be larger than normal phonon scattering times (τ

*) to allow momentum redistribution but smaller than resistive phonon scattering times (τ*

_{N}*) 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 (*

_{R}*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 (

*6*–

*10*). 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.

## RESULTS AND DISCUSSION

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 *R*_{spot} ≈ 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, ^{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*/*∂T*)Δ*T* + (*∂R*/*∂n*)Δ*n* (*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*/*∂T*)Δ*T* and (*∂R*/*∂n*)Δ*n*. 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} m^{3}, 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.

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 ^{exp} = 3 × 10^{−5} m^{2}/s, thus leading to a propagation velocity *20*–*22*) 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, *f*_{c} 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 *f*_{c} 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} is expected and indeed experimentally observed for lower temperatures. Wave-like effects are already present below *f*_{c}, 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.1*f*_{c}, 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 *f*_{c} relative to the frequency of the minimum originate from the relation between a, τ_{ss}, and *R*_{spot} (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

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: *t* is the temporal coordinate, and λ denotes each phonon mode. We note that the case with *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 *f* was then introduced into the linearized Boltzmann transport equation (BTE) to find the solution for _{ss}, in terms of individual phonon relaxation times (τ_{λ}) as shown in section S6_{λ} 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 *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, ∣Δ*T*∣_{max} < 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 *f*_{c}). 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.

## MATERIALS AND METHODS

### Samples

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/cm^{2}. 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/*e*^{2} 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 thirdorder.py (*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 MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/27/eabg4677/DC1

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.

## REFERENCES AND NOTES

**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.

- Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).