Introduction

The phonon dispersion of a crystal, which defines the frequency of normal lattice modes, is of fundamental importance because it determines thermal and mechanical properties1,2. Phonon dispersion in bulk crystals can be measured by inelastic neutron or X-ray scattering, whereby phonon frequency and momentum information is calculated from the energy and momentum exchange relations of neutron or X-ray beams3,4. However, these scattering approaches usually require large research facilities such as neutron sources, which are difficult to apply to nanoscale crystals due to their small scattering cross sections.

To examine the interface and fine structure hidden inside thin film heterostructures, picosecond acoustic measurements of photoexcited strain waves with high temporal resolution have been utilized5,6,7. Transverse acoustic (TA) wave propagation along the surface can be monitored by tracking the position of acoustic waves both spatially and temporally, but with a limited frequency range below a few GHz due to the spatial resolution allowed by the diffraction of light8,9,10. However, photoacoustic strain waves generated in a predefined nanoscale region, such as in a thin metal film transducer, can provide precise depth information of hidden nanostructures. Acoustic waves created in atomically narrow spatial regions, which convert to a short pulse duration in the time domain, however, comprise broad frequency components, and thus suffer from spatial and temporal broadening when propagating in a dispersive medium. Although the broadening may worsen the resolution of strain-wave nanoimaging, it may be useful for investigating the frequency-dependent behavior of crystal lattices.

Among various types of nanoscale crystals, two-dimensional (2D) van der Waals (vdW) materials have attracted considerable interest due to their unique and versatile functional properties, including the wide-ranging material selectivity of graphene, transition metal dichalcogenides, boron nitride (BN), and black phosphorus (BP)11,12. With atomically clean material interfaces, vdW heterostructures provide an excellent nanoscale platform to study strain wave behaviors at THz frequencies. For example, coherent lattice vibrations of optical and acoustic phonon modes have been reported in 2D vdW materials13,14,15,16,17,18,19. Notably, if a specific layer of atomic thickness is photoexcited, phonon modes with large momentum extending to the Brillouin zone edge can be accessed with the help of dimensional momentum in inverse proportion to the thickness20,21,22,23.

In this paper, we report a concept for measuring the phonon dispersion relation in nanoscale 2D vdW materials. This is demonstrated through the generation of broadband coherent acoustic phonons extending over 3 THz in few-layer BP crystals and their propagation along the stacking direction to adjacent hexagonal boron nitride (hBN) regions heterostructured with the BP. In a time-domain study utilizing a probe energy that the BP material is sensitive to, a strain wave is reflected at the air/hBN interface, generating echo signals that can be detected when they return to the BP layer. The strain wave group velocity dispersion (GVD) is obtained from spectral-temporal analysis of the photoacoustic strain wave. The results are readily explained by a linear chain model (LCM) calculation simplifying the interlayer coupling of BN as a spring constant connecting one-dimensional (1D) particle chains; i.e., the system is visualized as a 1D system of beads connected by springs (cf. Newton’s cradle). By measuring the frequency-dependent time-of-flight (TOF) of the strain waves, we demonstrate that the proposed method can be a versatile tool for investigating the phonon dispersion and nanoscale structure of 2D vdW materials.

Results

Optical pump–probe experiment

Figure 1 shows a schematic of the vdW heterostructure and the principles of time-resolved acoustic strain propagation induced by optical pump excitation employed in the experiment. The vdW heterostructure comprises a few layers of BP encapsulated between hBN layers (top and bottom hBN, abbreviated t-hBN and b-hBN) loaded onto a quartz substrate (Fig. 1a). Two samples were prepared with different t-hBN layer thicknesses of 12 nm and 29 nm (see the Methods section). Figure 1b illustrates the time-resolved evolution of a photoacoustic wave induced by ultrafast photoexcitation, which proceeds as follows. An optical pump pulse with a wavelength of 760 nm ( = 1.63 eV) excites carriers in the BP layer (step 1). A photoacoustic strain wave is generated due to inhomogeneous photoexcitation of the 760 nm pump pulse in the opaque BP layer24,25 (step 2). The photoacoustic strain pulse propagates into the adjacent hBN layers and is reflected back at the t-hBN/air interface, generating echoes that can be detected in the BP layer (steps 3 & 4). It should be noted that the waves traveling into b-hBN are either transmitted into the quartz substrate or scattered due to the surface roughness of SiO2 with minor contributions, as opposed to those reflecting from the t-hBN/air interface with a large impedance mismatch and a good flatness26. The dynamics of the photoacoustic strain pulse relating its round trips crossing the b-hBN/BP/t-hBN heterostructures is investigated by using an ultrafast time-resolved pump–probe method.

Fig. 1: Schematic of the van der Waals (vdW) heterostructure and time-resolved acoustic pulse propagation induced by optical pump excitation.
figure 1

a Schematic of the vdW heterostructure. Few-layer black phosphorus (BP) was encapsulated between bottom and top hexagonal boron nitride (hBN) layers. b Schematic diagrams of photoexcited lattice strain generation in BP and strain pulse propagation into hBN regions over time. Red arrows depict displacements of BP and hBN atoms induced by photoexcited lattice strain.

Full size image

A schematic of the pump–probe experiment used to measure the photoacoustic strain pulse in the vdW heterostructure is shown in Fig. 2a. Femtosecond pump pulse excitation induces a photoacoustic strain pulse in the BP layer. The optical transmission change (ΔT) of the probe pulse is measured with variation of the time delay with respect to the pump pulse24,25,27. The probe pulse selectively interacts with the BP layer because the probe energy is resonant with optical transitions in BP but transparent in the hBN layers. The time-resolved transient differential transmission (ΔT/T) obtained in the hBN/BP/hBN heterostructure with t-hBN thickness of 12 nm is shown in Fig. 2b. The positive value of ΔT near zero delay results from a reduced probe absorption due to the photobleaching effect28. The zoomed-in plot of ΔT/T in the inset of Fig. 2b exhibits a transmission change originated from the lattice vibrations in the BP layer29,30,31,32. When the exponential and slowly varying components are removed from the pump–probe signal (Fig. 2c), the signal exhibits strong strain pulse excitation at t = 0, and subsequent first and second echoes are observed within a 50 ps window.

Fig. 2: Time-resolved measurement of strain wave generation and propagation in the vdW heterostructure.
figure 2

a Schematic of the pump–probe method for detecting the photoacoustic strain pulse and its round trips in the vdW heterostructure of b-hBN/BP/t-hBN. b Photoinduced transmission change as a function of time delay for the heterostructure with a t-hBN thickness of 12 nm. The inset shows a zoomed-in plot of the strain wave signal superimposed on the pump–probe signal. c Coherent acoustic waves exhibiting multiple echoes extracted from the time-resolved transmission change. The inset shows a schematic of the strain wave propagation and round trips in the t-hBN region.

Full size image

The group velocity of the strain pulse can be obtained from the TOF measurement of round-trip echoes. Once the strain pulse TOF is spectrally resolved, it is possible to determine the acoustic phonon dispersion of the traveling medium from the group velocity. For this purpose, sliding window Fourier transform (SWFT) analysis of the time-domain strain wave was conducted with a Gaussian-shaped temporal window of 1.8 ps33. The resultant 2D contour plot shows the temporal evolution of the strain pulse spectrum (Fig. 3a). The strain signal is distributed over a wide spectral range from 0 to 3 THz, enabling experimental determination of the group velocity dispersion of hBN phonons. A large frequency dependence is apparent from the time-resolved FT intensity compared at two different frequencies (Fig. 3b): the TOF of the 1.5 THz strain wave (9.2 ps) is approximately 2 ps shorter than that of the 2.2 THz strain wave. As depicted in Fig. 1b, the strain wave is reflected at the air interface and travels through the t-hBN layer before making echoes. Because the optical probe at 760 nm is transparent in the hBN layers, the echo signals appear only when the strain wave is passing through the BP region. Thus, the strain wave group velocity can be assessed based on the echo times, and its dispersive characteristics in the frequency domain lead to the GVD of hBN.

Fig. 3: Frequency-resolved time-of-flight and group velocity dispersion of photoacoustic strain waves in hBN.
figure 3

a Sliding window Fourier transform (SWFT) spectrum of the heterostructure with 12 nm t-hBN (sample A). The solid lines show the frequency-dependent strain wave arrival time at the BP region predicted by linear chain model (LCM) calculations. b Temporal evolution of the FT intensity at two different frequencies, 1.5 and 2.2 THz, marked with colored horizontal dashed lines in Fig. 3a. c (left panel) Atomic structure of bulk hBN and its first Brillouin zone, denoted by Γ (0,0,0) and A (0,0,0.5) in reciprocal space, where the strain pulse propagates along the c axis. (right panel) Longitudinal acoustic (LA) phonon dispersion of hBN obtained from the LCM and density functional theory (DFT) employing the local density approximation (LDA). d Group velocity dispersion calculated from the LCM and DFT calculations (LDA and optB86b-vdW functionals).

Full size image

To verify the experimental observations in more detail, a 1D LCM was developed for the out-of-plane longitudinal acoustic (LA) mode of the hBN layers (see the Methods section for details). The numerically obtained arrival times of the strain pulse after one, two, and three round trips are delineated by the solid white lines in Fig. 3a. Excellent agreement with the experimental data indicates that the LCM model accurately describes the out-of-plane LA phonon dispersion of bulk hBN. Figure 3c shows the lattice structure and phonon dispersion of hBN along the c axis (the stacking direction of hBN layers), obtained from the LCM and density functional theory (DFT) simulations with the local density approximation (LDA) as the exchange-correlation functional. In Fig. 3d, the GVD values obtained from the LCM are compared with the DFT calculations employing the LDA or optB86b-vdW functional. Among the approximate exchange-correlation functionals of the DFT calculation, LDA exhibits the best match with the LCM results34, while optB86b-vdW shows a relatively good agreement with the LCM (see Supplementary Fig. 2 for benchmark calculations of exchange-correlation functionals for GVD predictability). According to a simple relation of λ = v/f, the wavelength λ of the strain wave at f = 3 THz is as small as 0.7 nm for the group velocity of v = 2 km/s in hBN. Thus, it is in principle possible to measure the TOF for the case of few nm thin hBN layers.

Analogous to the case of sample A in Fig. 3a, the SWFT plot of the time-domain strain wave for sample B (Fig. 4a) with t-hBN thickness of 29 nm is consistent with the GVD relation predicted by the LCM. As a light pulse confined in a laser cavity exhibits a frequency comb in the laser spectrum, the multiple echoes due to the phonon cavity in our 2D heterostructure are expected to modify the spectral features of the strain wave. The FT spectrum of the strain wave generated in sample B is presented in Fig. 4b, in which the t-hBN region constitutes the phonon cavity. The spectrum exhibits pronounced frequency combs, resulting from the sequential repetition of constructive and destructive interference of the strain wave echoes in the frequency domain. Because constructive interference occurs when the time interval between echoes coincides with a multiple of the strain wave period, its frequency (fm) can be described by the relation m/fm = 2d/vg, where m is an integer and d is the cavity length. In this case, vg can be determined from the frequency difference (∆f = fm+1 – fm) between adjacent comb modes, i.e., vg ≈ ∆f × 2d, with the reciprocal of ∆f being equal to the round-trip time in t-hBN. We note that the appearance of two types of combs in Fig. 4b is the consequence of multiple phonon echoes gradually decreasing with round trips, where the position of each comb peak is related to the round-trip spacing. As shown in the fm versus m plot in Fig. 4c, the slope of the fm curve decreases as the frequency increases, resulting in smaller ∆f at high frequencies (Fig. 4d). The group velocity deduced from ∆f in the frequency comb assuming d = 28.6 nm agrees well with the GVD from the LCM calculation. This result indicates that the frequency combs based on time-resolved strain wave analysis successfully reproduce the phonon dispersion in the thin film heterostructure.

Fig. 4: Frequency-resolved time-of-flight and frequency comb signal for sample B with a t-hBN thickness of 29 nm.
figure 4

a Sliding window Fourier transform (SWFT) spectrum of sample B with 29 nm t-hBN. The solid lines indicate the frequency-dependent strain wave arrival time at the BP region predicted by the linear chain model (LCM). b FT spectrum of the time-domain coherent acoustic phonon signal of sample B. The insets show expanded views of the shaded regions highlighting the frequency combs. The numbers in the insets are the relative indices of the combs, and the red indices indicate the presence of two peaks within one comb width. c Frequency versus relative index plot. The frequencies of the red indices in Fig. 4b are plotted as red squares, while those of the black indices are plotted as black squares. d Group velocity dispersion calculated from the frequency difference (Δf) of adjacent combs. The solid line denotes the LCM calculation.

Full size image

First principles simulation

Although the photoacoustic strain pulse generation and propagation have experimentally led to the phonon dispersion of bulk hBN, a microscopic explanation of the strain generation and propagation in the hBN layers is necessary. We performed ab initio simulations to investigate the quantum origin of the photoacoustic strain pulse generation and obtain an atomistic description of the strain pulse propagation. The electronic band structure of t-hBN/BP/b-hBN calculated by DFT is shown in the left panel of Fig. 5a. A finite-sized heterostructure of 6L-hBN/4L-BP/6L-hBN was considered for mechanical and quantum interactions between BP and hBN layers to ensure a reasonable computational cost for practical ab initio simulations. The band structure was calculated with the optB86b-vdW functional, which reveals a direct band gap at the Γ point, as reported in previous studies for bulk and few-layer BP27,35,36. The optB86b-vdW functional was utilized because it can reliably predict both the electronic structure and phonon dispersion of the hBN/BP/hBN heterostructure and bulk hBN, whereas the LDA functional severely underestimates the band gap of 2D vdW materials (Supplementary Fig. 3 for benchmark calculations of the electronic band structure and phonon dispersion). The charge density profile of the band edge states at the Γ point is displayed in the top-right panel of Fig. 5a. An occupation constraint was applied in the delta self-consistent field (ΔSCF) calculation to include photoexcited hot carrier generation24 (see the Methods section). In contrast to the conduction band minimum (CBM), the valence band maximum (VBM) exhibits an inhomogeneous charge distribution in the inner two BP layers. As shown in the bottom-right panel of Fig. 5a, photoinduced charge depletion from the inner BP layers results in a strain wave through atomic relaxations. This sudden lattice strain upon photoexcitation can generate an acoustic pulse in the BP region on a sub-picosecond timescale, which then propagates into adjacent hBN layers and creates multiple echo signals during successive round trips. After confirming the electronic origin of the photoacoustic strain pulse generated in the BP, a real-time atomistic scale description of the acoustic pulse propagation crossing the BP/hBN interface was calculated by molecular dynamics (MD) simulation based on a machine-learned force field (MLFF) (see the Methods section for a detailed description). Figure 5b presents the simulated atomic displacements at different times, demonstrating strain pulse propagation into the hBN layers from BP. Real-time travel and reflection of the strain pulse in hBN were simulated for the heterostructure b-hBN/4L-BP/t-hBN (12 nm). During the time period t = 0–10 ps, the strain pulse split into several pulses of different durations. As displayed in the phonon dispersion and GVD plots in Fig. 3c, d, the group velocity decreased as the phonon frequency increased. This dispersive behavior induced a splitting of the strain pulse into several wider (faster) and narrower (slower) wave packets associated with smaller and larger phonon momentum, respectively. Real-time propagation of the strain pulse in the t-hBN layer viewed on the atomistic scale is provided as a Supplementary Movie. The phonon power spectra of the hBN and 4L-BP are provided in Supplementary Fig. 6, which indicates that the acoustic phonon in the photoexcited 4L-BP transfers efficiently to the hBN along the stacking direction.

Fig. 5: Density functional simulation of photoacoustic strain generation and real-time travel in the hBN/BP/hBN heterostructure.
figure 5

a (Left panel) Band structure of 6L-hBN/4L-BP/6L-hBN and the delta self-consistent field (ΔSCF) occupation constraint applied to the valence band maximum (VBM) and the conduction band minimum (CBM) states at the Γ point. (Right panel) Plane-averaged charge densities of the VBM and CBM states and photoinduced strain resulting from lattice relaxation after photoexcited carrier generation. b Simulated atomic displacements over time generated by the photoacoustic strain pulse propagation through the 12 nm t-hBN/4BP/b-hBN heterostructure. c Experimental ΔT/T signals and simulated real-time phosphorus atomic displacements obtained for two heterostructures with different t-hBN thicknesses of 12 nm and 29 nm.

Full size image

The real-time strain wave amplitude simulations were then compared with the experimental results of ∆T/T for two different t-hBN thicknesses: 12 nm (sample A) and 29 nm (sample B) (Fig. 5c). In the simulation, 4L-BP instead of 6L-BP was adopted for sample B to avoid spurious energy-gap closing of the semi-local functional. The atomic displacements of BP obtained from the MLFF simulation exhibit quantitively good agreement with the experimental differential transmission regarding the oscillatory behavior as well as the first phonon echo for both samples, indicating that the differential transmission reflects the time-transient atomic displacements (i.e., lattice strain) of the BP layer37,38. Note that the overall decay of the acoustic phonon signal with time is caused by acoustic wave loss due to either energy penetration into the quartz substrate or phonon scattering at the heterointerfaces. In both simulation and experiment, the echo times of the strain pulse linearly correlate with the t-hBN layer thickness, e.g., the first echo occurred at 7.8 ps and 17.0 ps for sample A and sample B, respectively. This linearity indicates that time-resolved detection of photoexcited strain pulses may provide a nondestructive tool for inspecting nanosized vdW heterostructures of several atomic layer thickness.

Discussion

Time-resolved measurement of photoacoustic pulse propagation in vdW b-hBN/BP/t-hBN heterostructures was conducted in tandem with frequency-dependent time-of-flight (TOF) analysis of the propagation and echoes through the t-hBN layer to determine the full group velocity dispersion (GVD) of LA phonons along the stacking direction. The material-specific optical response of the probe transmittance, being sensitive to the strain wave originating from the BP, enabled observation of the arrival of the strain wave back at the BP after round-trip travel within the heterostructure. The GVD of the out-of-plane LA acoustic wave mode was obtained from the phonon transit time as well as the energy separation in the frequency comb of the Fourier-transformed spectrum. The reliability of the phonon dispersion data deduced from the photoacoustics was confirmed by comparison with linear chain model simulations, in which the hBN layers were simplified as a linear mass chain connected by springs. We also revealed the microscopic origin of photoinduced strain pulse generation and its propagation from first principles calculation. The time-resolved THz acoustic wave analysis demonstrated in this work may provide a nondestructive tool for investigating the hidden interfaces and lattice properties of nanoscale thin film heterostructures. It is expected that our method utilizing photoexcitation within a several atomic layer thick sensing layer can be extended to various nanoscale systems, including van der Waals heterostructures and semiconductor quantum wells.

Methods

Sample preparation

Encapsulated BP heterostructures were prepared in an argon-filled glove box (oxygen and water content below 0.5 ppm) by sequentially transferring bottom hBN (b-hBN), BP, and top hBN (t-hBN) flakes onto quartz substrates with a well-established Gel-Pak dry transfer method. Two t-hBN/BP/b-hBN samples were fabricated with thicknesses of 12 nm/1.3 nm/30 nm (sample A) and 29 nm/2.3 nm/27 nm (sample B). The thickness of hBN was estimated from thin film interference in a reflectance spectrum, whereas the thickness of BP was measured from the surface morphology determined by atomic force microscopy (AFM).

Pump–probe measurement

To study the generation and propagation of coherent acoustic phonons in the hBN/BP/hBN heterostructures, a pump–probe experiment was performed using ultrashort pulses of 60 fs duration supplied by a 1 MHz repetition-rate cavity-pumped Ti:sapphire laser operating at a center wavelength of 760 nm. A beam splitter was used to divide the pulse into a relatively strong pump pulse and a weak probe pulse (Fig. 1a). Because hBN layers are transparent in the near-infrared, the pump pulse excites carriers and coherent phonons selectively in the BP layer. The probe pulse was used to monitor the pump-induced optical modulation in the BP layer as a function of the time delay relative to the pump pulse. To maximize the light–matter interaction in anisotropic BP, the polarization of both pump and probe beams was set to be parallel to the armchair crystal axis of BP. Using a 20× objective lens, 1 nJ pump and 0.5 nJ probe pulses overlapped in a 5 μm diameter spot in the sample. The time delay between the pump and probe pulses was scanned at 13 Hz, and the signal was acquired until the noise was sufficiently suppressed to obtain a clear oscillation signal of coherent acoustic phonons.

Linear chain model (LCM) for the out-of-plane longitudinal acoustic (LA) mode

Based on symmetry considerations, thermal expansion of the BP layer upon photoexcitation was expected to generate a strong out-of-plane LA mode. The macroscopic lattice behavior of the hBN layers was simplified with a 1D LCM, where each layer was modeled as a mass particle and the interlayer interaction as a spring force39. According to the LCM, the phonon dispersion is given by ω(q) = ω0 × sin(q/2c), where ω0, q, and c are the maximum angular frequency at the Brillouin zone edge, the wavevector, and the interlayer spacing of hBN (c = 3.33 Å), respectively. The group velocity, vg, is calculated as vg(q) = v0 × cos(q/2c), where v0 = ω0/2c is the sound velocity in the long wavelength limit. The strain wave echo time follows the relation Te = 2d/ vg, where the effective t-hBN thickness d is used to include the effect of finite BP thickness. The value of v0 was assumed to be v0 = 3.44 km/s based on literature from inelastic X-ray scattering40.

Electronic band structure calculation

Ab initio phonon calculations were carried out for bulk hBN, few-layer BP, and their vertically stacked structures using DFT with the optB86b-vdW exchange-correlation functional as implemented in the VASP package41,42,43,44. Phonon dispersions, harmonic eigenmodes, and phonon group velocities were calculated using the phonopy code45. The optB86b-vdW functional46,47 was used to calculate both the finite band gap and group velocities to compare with experimental results. Benchmark calculations of exchange-correlation functionals for predictions of electronic structures and group velocities are provided in the Supplementary Information.

Photoinduced strain calculation

To model the light-induced strain in the hBN/BP/hBN heterostructure, atomic displacements induced by carrier excitation were calculated using an occupation-constrained ΔSCF approach for a slab model consisting of hBN (6 layers)/BP (4 layers)/hBN (6 layers) by imposing a photoexcited non-Aufbau hole and electron occupation at the Γ point with fixed occupation of the valence and conduction bands, as illustrated in Fig. 4b48,49,50,51,52. The excited carrier imposed by the occupation constraint was estimated as 6.36 × 1012 cm−2, which is a moderate photoexcitation density49,53,54.

Machine-learned force field (MLFF) construction

An MLFF was trained by on-the-fly molecular dynamics (MD) simulations of an isothermal-isobaric molecular ensemble at 300 K for 100,000 steps with a 2 ps time step using the Langevin thermostat55,56, as implemented in the VASP package. A 4 × 4 × 4 supercell of bulk hBN and 6L-hBN/4L-BP/6L-hBN atomic geometries were subsequently used in the MLFF training. Details of the on-the-fly construction of MLFFs can be found in the literature57,58,59,60,61. The present study employed cutoffs of 5 Å and 8 Å for the radial and angular descriptors of the Gaussian approximation potential62,63, respectively, and a weighting factor of 1000 was imposed on the ionic forces. Validation of the MLFF-predicted results compared to ab initio counterparts is presented in Supplementary Fig. 4.

Photoacoustic molecular dynamics simulation based on an MLFF approach

After calculating the light-induced initial strain, the real-time propagation of the strain pulse was simulated via MD simulation using the MLFF, as implemented in the on-the-fly machine learning routine in the VASP package57,58,59. Atomic trajectories were integrated with the Verlet algorithm under an energy conservation constraint with a time step of 0.5 ps and total simulation time of 50 ps (total 100,000 steps). The atomic displacements at each time step were convoluted with a Lorentzian function with a line width of 1 Å for visualization of the wave packets of the strain pulse. Vibrational power spectra of the BP and hBN layers were obtained using the DynaPhoPy package by calculating the FT spectra of velocity autocorrelation functions64. For the MD simulations, an absorbing boundary condition was used for the bottom hBN layer.