# Fileset

[s41699-024-00475-8.pdf](https://mdr.nims.go.jp/filesets/8d9f00b9-647a-4b36-8898-72eb05a6928f/download)

## Creator

Seong-Yeon Lee, Soungmin Bae, Seonyeong Kim, Suyong Jung, [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), Hannes Raebiger, Ki-Ju Yee

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Full phonon dispersion along the stacking direction in nanoscale van der Waals materials by picosecond acoustics](https://mdr.nims.go.jp/datasets/448bdab8-b215-441a-94f0-7ab16defec68)

## Fulltext

Full phonon dispersion along the stacking direction in nanoscale van der Waals materials by picosecond acousticsnpj | 2D materials and applications ArticlePublished in partnership with FCT NOVA with the support of E-MRShttps://doi.org/10.1038/s41699-024-00475-8Full phonon dispersion along the stackingdirection in nanoscale van der Waalsmaterials by picosecond acousticsCheck for updatesSeong-Yeon Lee 1,7, Soungmin Bae 2,7 , Seonyeong Kim3, Suyong Jung 3, Kenji Watanabe 4,Takashi Taniguchi 5, Hannes Raebiger 6 & Ki-Ju Yee 1Phonon dispersion in crystals determines many important material properties, but its measurementusually requires large-scale facilities and is limited to bulk samples. Here, we demonstrate themeasurement of full phonon dispersion along the stacking direction in nanoscale systems by usingpicosecond acoustics. A heterostructure sample was prepared consisting of layers of hexagonalboron nitride (hBN) sandwiching a thin layer of black phosphorus (BP), within which a strain pulse wasgenerated by photoexcitation and observed with an optical probe in the BP layer. The strain pulsetraverses to the fewnanometer thick hBN layers,where it propagates to the edgeandechoesback, likeacoustic waves in Newton’s cradle. The echoes returning to the BP layer provide information on thefrequency-dependent time-of-flight and group velocity dispersion of the sample system. Themicroscopic origin of the photoinduced strain pulse generation and its propagation is revealed fromfirst principles. Phonon frequency combs observed in the Fourier transform spectrum confirm thestrain wave round trips and demonstrate the feasibility of determining group velocity dispersionthrough photoacoustics.The phonon dispersion of a crystal, which defines the frequency of normallattice modes, is of fundamental importance because it determines thermaland mechanical properties1,2. Phonon dispersion in bulk crystals can bemeasured by inelastic neutron or X-ray scattering, whereby phonon fre-quency and momentum information is calculated from the energy andmomentum exchange relations of neutron or X-ray beams3,4. However,these scattering approaches usually require large research facilities such asneutron sources, which are difficult to apply to nanoscale crystals due totheir small scattering cross sections.To examine the interface and fine structure hidden inside thin filmheterostructures, picosecond acoustic measurements of photoexcited strainwaves with high temporal resolution have been utilized5–7. Transverseacoustic (TA) wave propagation along the surface can be monitored bytracking the position of acoustic waves both spatially and temporally, butwith a limited frequency range below a fewGHzdue to the spatial resolutionallowed by the diffraction of light8–10. However, photoacoustic strain wavesgenerated in a predefined nanoscale region, such as in a thin metal filmtransducer, can provide precise depth information of hidden nanos-tructures. 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 andtemporal broadening when propagating in a dispersive medium. Althoughthe broadening may worsen the resolution of strain-wave nanoimaging, itmay be useful for investigating the frequency-dependent behavior of crystallattices.Among various types of nanoscale crystals, two-dimensional (2D) vanderWaals (vdW)materials have attracted considerable interest due to theirunique and versatile functional properties, including the wide-rangingmaterial selectivity of graphene, transition metal dichalcogenides, boronnitride (BN), and blackphosphorus (BP)11,12.With atomically cleanmaterialinterfaces, vdWheterostructures provide an excellent nanoscale platform tostudy strain wave behaviors at THz frequencies. For example, coherentlattice vibrations of optical and acoustic phononmodes have been reportedin 2D vdW materials13–19. Notably, if a specific layer of atomic thickness is1Department of Physics and Institute of Quantum Systems, Chungnam National University, Daejeon 34134, Republic of Korea. 2Institute for Materials Research,Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan. 3Interdisciplinary Materials Measurement Institute, Korea Research Institute of Standardsand Science, Daejeon 34113, Republic of Korea. 4Research Center for Electronic and Optical Materials, National Institute for Materials Science, 1-1 Namiki,Tsukuba 305-0044, Japan. 5Research Center for Materials Nanoarchitectonics, National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan.6Department of Physics, Yokohama National University, Yokohama 240-8501, Japan. 7These authors contributed equally: Seong-Yeon Lee, Soungmin Bae.e-mail: bae.soungmin.d6@tohoku.ac.jp; kyee@cnu.ac.krnpj 2D Materials and Applications |            (2024) 8:39 11234567890():,;1234567890():,;http://crossmark.crossref.org/dialog/?doi=10.1038/s41699-024-00475-8&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41699-024-00475-8&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41699-024-00475-8&domain=pdfhttp://orcid.org/0000-0001-5628-2914http://orcid.org/0000-0001-5628-2914http://orcid.org/0000-0001-5628-2914http://orcid.org/0000-0001-5628-2914http://orcid.org/0000-0001-5628-2914http://orcid.org/0000-0001-9922-444Xhttp://orcid.org/0000-0001-9922-444Xhttp://orcid.org/0000-0001-9922-444Xhttp://orcid.org/0000-0001-9922-444Xhttp://orcid.org/0000-0001-9922-444Xhttp://orcid.org/0000-0003-3885-9999http://orcid.org/0000-0003-3885-9999http://orcid.org/0000-0003-3885-9999http://orcid.org/0000-0003-3885-9999http://orcid.org/0000-0003-3885-9999http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0003-3969-9165http://orcid.org/0000-0003-3969-9165http://orcid.org/0000-0003-3969-9165http://orcid.org/0000-0003-3969-9165http://orcid.org/0000-0003-3969-9165http://orcid.org/0000-0002-1076-2354http://orcid.org/0000-0002-1076-2354http://orcid.org/0000-0002-1076-2354http://orcid.org/0000-0002-1076-2354http://orcid.org/0000-0002-1076-2354mailto:bae.soungmin.d6@tohoku.ac.jpmailto:kyee@cnu.ac.krphotoexcited, phonon modes with large momentum extending to theBrillouin zone edge can be accessed with the help of dimensionalmomentum in inverse proportion to the thickness20–23.In this paper, we report a concept formeasuring the phonondispersionrelation in nanoscale 2D vdWmaterials. This is demonstrated through thegeneration of broadband coherent acoustic phonons extending over 3 THzin few-layerBPcrystals and their propagationalong the stackingdirection toadjacent hexagonal boron nitride (hBN) regions heterostructured with theBP. In a time-domain study utilizing a probe energy that the BP material issensitive to, a strain wave is reflected at the air/hBN interface, generatingecho signals that canbedetectedwhen they return to theBP layer. The strainwave group velocity dispersion (GVD) is obtained from spectral-temporalanalysis of the photoacoustic strain wave. The results are readily explainedby a linear chain model (LCM) calculation simplifying the interlayer cou-pling of BN as a spring constant connecting one-dimensional (1D) particlechains; i.e., the system is visualized as a 1D system of beads connected bysprings (cf. Newton’s cradle). Bymeasuring the frequency-dependent time-of-flight (TOF) of the strain waves, we demonstrate that the proposedmethod can be a versatile tool for investigating the phonon dispersion andnanoscale structure of 2D vdW materials.ResultsOptical pump–probe experimentFigure 1 shows a schematic of the vdWheterostructure and the principles oftime-resolved acoustic strain propagation induced by optical pump exci-tation employed in the experiment. The vdW heterostructure comprises afew layers of BP encapsulated between hBN layers (top and bottom hBN,abbreviated t-hBNandb-hBN) loadedonto aquartz substrate (Fig. 1a).Twosamples were prepared with different t-hBN layer thicknesses of 12 nm and29 nm (see the Methods section). Figure 1b illustrates the time-resolvedevolution of a photoacoustic wave induced by ultrafast photoexcitation,which proceeds as follows. An optical pump pulse with a wavelength of760 nm ( = 1.63 eV) excites carriers in the BP layer (step 1). A photoacousticstrain wave is generated due to inhomogeneous photoexcitation of the760 nm pump pulse in the opaque BP layer24,25 (step 2). The photoacousticstrain pulse propagates into the adjacent hBN layers and is reflected back atthe t-hBN/air interface, generating echoes that can be detected in the BPlayer (steps 3&4). It should benoted that thewaves traveling intob-hBNareeither transmitted into the quartz substrate or scattered due to the surfaceroughness of SiO2 with minor contributions, as opposed to those reflectingfrom the t-hBN/air interface with a large impedance mismatch and a goodflatness26. The dynamics of the photoacoustic strain pulse relating its roundtrips crossing the b-hBN/BP/t-hBNheterostructures is investigated byusingan ultrafast time-resolved pump–probe method.A schematic of the pump–probe experiment used to measure thephotoacoustic strain pulse in the vdW heterostructure is shown in Fig. 2a.Femtosecond pump pulse excitation induces a photoacoustic strain pulse inthe BP layer. The optical transmission change (ΔT) of the probe pulse ismeasured with variation of the time delay with respect to the pumppulse24,25,27. The probe pulse selectively interacts with the BP layer becausethe probe energy is resonantwith optical transitions inBPbut transparent inthe hBN layers. The time-resolved transient differential transmission (ΔT/T) obtained in the hBN/BP/hBN heterostructure with t-hBN thickness of12 nm is shown in Fig. 2b. The positive value of ΔT near zero delay resultsfrom a reduced probe absorption due to the photobleaching effect28. Thezoomed-in plot ofΔT/T in the inset of Fig. 2b exhibits a transmission changeoriginated from the lattice vibrations in the BP layer29–32. When the expo-nential and slowly varying components are removed from the pump–probesignal (Fig. 2c), the signal exhibits strong strain pulse excitation at t = 0, andsubsequent first and second echoes are observed within a 50 ps window.The group velocity of the strain pulse can be obtained from the TOFmeasurement of round-trip echoes. Once the strain pulse TOF is spectrallyresolved, it is possible to determine the acoustic phonon dispersion of thetravelingmedium from the group velocity. For this purpose, slidingwindowFourier transform (SWFT) analysis of the time-domain strain wave wasconducted with a Gaussian-shaped temporal window of 1.8 ps33. Theresultant 2D contour plot shows the temporal evolution of the strain pulsespectrum (Fig. 3a). The strain signal is distributed over awide spectral rangefrom0 to 3 THz, enabling experimental determination of the group velocitydispersionof hBNphonons.A large frequencydependence is apparent fromthe 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 psshorter than that of the 2.2 THz strainwave.Asdepicted inFig. 1b, the strainwave is reflected at the air interface and travels through the t-hBN layerbeforemaking echoes. Because the optical probe at 760 nm is transparent inthe hBN layers, the echo signals appear only when the strain wave is passingthrough the BP region. Thus, the strain wave group velocity can be assessedbased on the echo times, and its dispersive characteristics in the frequencydomain lead to the GVD of hBN.To verify the experimental observations inmore detail, a 1D LCMwasdeveloped for the out-of-plane longitudinal acoustic (LA)mode of the hBNlayers (see theMethods section for details). Thenumerically obtainedarrivaltimes of the strain pulse after one, two, and three round trips are delineatedby the solidwhite lines in Fig. 3a. Excellent agreementwith the experimentaldata indicates that the LCMmodel accurately describes the out-of-plane LAphonon dispersion of bulk hBN. Figure 3c shows the lattice structure andphonon dispersion of hBN along the c axis (the stacking direction of hBNlayers), obtained from the LCM and density functional theory (DFT)SubstrateBlack phosphorus(4L-BP)bottom-hBN(b-hBN)top-hBN(t-hBN)SiO2Light excitationa b Round trip & echoTime evolution1. Carrierphotoexciation3. Strain pulse propagation2. Photoacousticlattice strain~sub fs ~sub ps ~ psElectronHoleStrain pulsepropagationLattice strain4. Strain pulse echoPumpSubstrateAirvdW heterostructureFig. 1 | Schematic of the van derWaals (vdW) heterostructure and time-resolvedacoustic pulse propagation induced by optical pump excitation. a Schematic ofthe vdW heterostructure. Few-layer black phosphorus (BP) was encapsulatedbetween bottom and top hexagonal boron nitride (hBN) layers. b Schematicdiagrams of photoexcited lattice strain generation inBP and strain pulse propagationinto hBN regions over time. Red arrows depict displacements of BP and hBN atomsinduced by photoexcited lattice strain.https://doi.org/10.1038/s41699-024-00475-8 Articlenpj 2D Materials and Applications |            (2024) 8:39 2simulations with the local density approximation (LDA) as the exchange-correlation functional. In Fig. 3d, the GVD values obtained from the LCMare compared with the DFT calculations employing the LDA or optB86b-vdW functional. Among the approximate exchange-correlation functionalsof theDFT calculation, LDA exhibits the bestmatchwith the LCM results34,while optB86b-vdW shows a relatively good agreement with the LCM (seeSupplementary Fig. 2 for benchmark calculations of exchange-correlationfunctionals for GVDpredictability). According to a simple relation of λ = v/f, thewavelengthλof the strainwave at f = 3 THz is as small as 0.7 nm for thegroup velocity of v = 2 km/s in hBN. Thus, it is in principle possible tomeasure the TOF for the case of few nm thin hBN layers.Analogous to the case of sampleA inFig. 3a, the SWFTplot of the time-domain strain wave for sample B (Fig. 4a) with t-hBN thickness of 29 nm isconsistent with the GVD relation predicted by the LCM. As a light pulseconfined in a laser cavity exhibits a frequency comb in the laser spectrum,the multiple echoes due to the phonon cavity in our 2D heterostructure areexpected tomodify the spectral features of the strainwave. The FT spectrumof the strain wave generated in sample B is presented in Fig. 4b, in which thet-hBN region constitutes the phonon cavity. The spectrum exhibits pro-nounced frequency combs, resulting from the sequential repetition ofconstructive and destructive interference of the strain wave echoes in thefrequency domain. Because constructive interference occurs when the timeinterval between echoes coincides with a multiple of the strain wave period,its frequency (fm) can be described by the relationm/fm = 2d/vg, wherem isan integer and d is the cavity length. In this case, vg can be determined fromthe frequency difference (Δf= fm+1– fm) between adjacent combmodes, i.e.,vg ≈Δf × 2d, with the reciprocal of Δf being equal to the round-trip time int-hBN. We note that the appearance of two types of combs in Fig. 4b is theconsequence of multiple phonon echoes gradually decreasing with roundtrips, where the position of each comb peak is related to the round-tripspacing. As shown in the fm versusm plot in Fig. 4c, the slope of the fm curvedecreases as the frequency increases, resulting in smaller Δf at high fre-quencies (Fig. 4d). The group velocity deduced from Δf in the frequencycomb assuming d = 28.6 nm agrees well with the GVD from the LCMcalculation. This result indicates that the frequency combs based on time-resolved strain wave analysis successfully reproduce the phonon dispersionin the thin film heterostructure.First principles simulationAlthough the photoacoustic strain pulse generation and propagation haveexperimentally led to the phonon dispersion of bulk hBN, a microscopicexplanation of the strain generation and propagation in the hBN layers isnecessary. We performed ab initio simulations to investigate the quantumorigin of the photoacoustic strain pulse generation and obtain an atomisticdescription of the strain pulse propagation. The electronic band structure oft-hBN/BP/b-hBN calculated by DFT is shown in the left panel of Fig. 5a. Afinite-sized heterostructure of 6L-hBN/4L-BP/6L-hBN was considered formechanical andquantum interactions betweenBPandhBN layers to ensurea reasonable computational cost forpractical ab initio simulations.Thebandstructure was calculated with the optB86b-vdW functional, which reveals adirect band gap at the Γ point, as reported in previous studies for bulk andfew-layer BP27,35,36. The optB86b-vdW functional was utilized because it canreliably predict both the electronic structure and phonon dispersion of thehBN/BP/hBN heterostructure and bulk hBN, whereas the LDA functionalseverely underestimates the band gap of 2D vdWmaterials (SupplementaryFig. 3 for benchmark calculations of the electronic band structure andphonondispersion). The charge density profile of the band edge states at theΓ point is displayed in the top-right panel of Fig. 5a. An occupation con-straint was applied in the delta self-consistent field (ΔSCF) calculation toFig. 2 | Time-resolved measurement of strain wave generation and propagationin the vdW heterostructure. a Schematic of the pump–probemethod for detecting thephotoacoustic strain pulse and its round trips in the vdWheterostructure of b-hBN/BP/t-hBN. b Photoinduced transmission change as a function of time delay for the hetero-structurewith a t-hBN thickness of 12 nm.The inset shows a zoomed-inplot of the strainwave signal superimposed on the pump–probe signal. c Coherent acoustic waves exhi-biting multiple echoes extracted from the time-resolved transmission change. The insetshows a schematic of the strain wave propagation and round trips in the t-hBN region.Wavevector qΓ01234Frequency (THz)LDALCM(0,0,0) A(0,0,0.5)←0 10 20 30 40 500123Time delay t (ps)Frequency (THz)0maxLinear chain model (LCM) 1st echo 2nd echo 3rd echo2.2 THz1.5 THzΔt= 11.0 psΔt= 9.1 psabFT intensity (arb. uint)0 10 20 30 40 50Time delay t (ps)1.5 THz2.2 THz2.2 THzΔt = 11.0 ps, Vg= 2.54 km/s1.5 THzΔt = 9.2 ps, Vg= 3.04 km/s  LA phonon (c axis)c dLDALCMoptB86b-vdwGroup velocity vg (km/s)Frequency (THz)0 1 2 3 401234ΓKMAc*(0,0,1)←a*←b*← Γ: (0, 0, 0)A: (0, 0, 0.5)First Brillouin zone c= 3.33 Åa←b←c←Bulk hBNFig. 3 | Frequency-resolved time-of-flight and group velocity dispersion ofphotoacoustic strain waves in hBN. a Sliding window Fourier transform (SWFT)spectrum of the heterostructure with 12 nm t-hBN (sample A). The solid lines showthe frequency-dependent strain wave arrival time at the BP region predicted bylinear chain model (LCM) calculations. b Temporal evolution of the FT intensity attwo different frequencies, 1.5 and 2.2 THz, marked with colored horizontal dashedlines in Fig. 3a. c (left panel)Atomic structure of bulk hBNand itsfirst Brillouin zone,denoted by Γ (0,0,0) and A (0,0,0.5) in reciprocal space, where the strain pulsepropagates along the c axis. (right panel) Longitudinal acoustic (LA) phonon dis-persion of hBN obtained from the LCM and density functional theory (DFT)employing the local density approximation (LDA). d Group velocity dispersioncalculated from the LCM and DFT calculations (LDA and optB86b-vdWfunctionals).https://doi.org/10.1038/s41699-024-00475-8 Articlenpj 2D Materials and Applications |            (2024) 8:39 3include photoexcited hot carrier generation24 (see the Methods section). Incontrast to the conduction band minimum (CBM), the valence bandmaximum (VBM) exhibits an inhomogeneous charge distribution in theinner two BP layers. As shown in the bottom-right panel of Fig. 5a, pho-toinduced charge depletion from the inner BP layers results in a strain wavethrough atomic relaxations. This sudden lattice strain upon photoexcitationcan generate an acoustic pulse in the BP region on a sub-picosecondtimescale, which then propagates into adjacent hBN layers and createsmultiple echo signals during successive round trips. After confirming theelectronicoriginof thephotoacoustic strainpulse generated in theBP, a real-time atomistic scale description of the acoustic pulse propagation crossingthe BP/hBN interface was calculated by molecular dynamics (MD) simu-lation based on a machine-learned force field (MLFF) (see the Methodssection for a detailed description). Figure 5b presents the simulated atomicdisplacements at different times, demonstrating strain pulse propagationinto the hBN layers from BP. Real-time travel and reflection of the strainpulse 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 intoseveral pulses of different durations. As displayed in the phonon dispersionand GVD plots in Fig. 3c, d, the group velocity decreased as the phononfrequency increased. This dispersive behavior induced a splitting of thestrain pulse into several wider (faster) and narrower (slower) wave packetsassociated with smaller and larger phonon momentum, respectively. Real-time propagation of the strain pulse in the t-hBN layer viewed on theatomistic scale is provided as a Supplementary Movie. The phonon powerspectra of the hBN and 4L-BP are provided in Supplementary Fig. 6, whichindicates that the acoustic phonon in the photoexcited 4L-BP transfersefficiently to the hBN along the stacking direction.The real-time strain wave amplitude simulations were then comparedwith 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-BPinstead of 6L-BP was adopted for sample B to avoid spurious energy-gapclosing of the semi-local functional. The atomic displacements of BPobtained from the MLFF simulation exhibit quantitively good agreementwith the experimental differential transmission regarding the oscillatorybehavior as well as the first phonon echo for both samples, indicating thatthe differential transmission reflects the time-transient atomic displace-ments (i.e., lattice strain) of theBP layer37,38.Note that the overall decayof theacoustic phonon signal with time is caused by acoustic wave loss due toeither energy penetration into the quartz substrate or phonon scattering atthe heterointerfaces. In both simulation and experiment, the echo times ofRelative index,m b c d0 1 2 3 4Frequency (THz)FT intentisty (arb. unit)b-hBN/6L-BP/82L-hBN (29 nm)1.0 1.201 23 4 5 6 7765432102.4 2.63029 31323334Frequency (THz)Group velocity vg (km/s)Δfcomb×2d (exp.)0 1 2 3 4 5012345LCMcomb-2comb-1ΔfcombFrequency (THz)01230−20 20 400 10 20 30 40 500123Time delay t (ps)Frequency (THz)0max1st echo2nd echoLinear chain model (LCM) aFig. 4 | Frequency-resolved time-of-flight and frequency comb signal for sampleB with a t-hBN thickness of 29 nm. 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 chainmodel (LCM). b FT spectrum of the time-domain coherent acoustic phonon signalof sample B. The insets show expanded views of the shaded regions highlighting thefrequency combs. The numbers in the insets are the relative indices of the combs, andthe 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. 4bare plotted as red squares, while those of the black indices are plotted as blacksquares. dGroup velocity dispersion calculated from the frequency difference (Δf) ofadjacent combs. The solid line denotes the LCM calculation.a Photoexcitation calculationb Pulse propagation and round trip00.10.200.10.2–20 –10 0 10 20Atom position (Å)CBM (Γ-point)VBM (Γ-point)Plane-averged charge densityChargen density (Å–1)-0.4-0.2 0 0.2 0.4 0.6Γ XEnergy (eV)OccupiedUnoccupiedPhotoexcitation(ΔSCF)Photoinduced strain4BPb-BNReflectionPropagationt = 0.0 pst = 1.0 pst = 2.0 pst = 3.0 pst = 4.0 pst = 5.0 pst = 6.0 pst = 7.0 psTime-resolved strain (arb. unit)ΔT/T  & lattice stra (arb.unit)0 10 20 30 40 50b-hBN/4L-BP/82L-hBN (29 nm), simulationb-hBN/4L-BP/36L-hBN (12 nm), simulationb-hBN/6L-BP/82L-hBN (29 nm), experimentb-hBN/4L-BP/36L-hBN (12 nm), experimentTime (ps)16.6 ps7.8 pscSample ASample BFig. 5 | Density functional simulation of photoacoustic strain generation andreal-time travel in the hBN/BP/hBNheterostructure. a (Left panel) Band structureof 6L-hBN/4L-BP/6L-hBN and the delta self-consistent field (ΔSCF) occupationconstraint applied to the valence band maximum (VBM) and the conduction bandminimum (CBM) states at the Γ point. (Right panel) Plane-averaged charge densitiesof the VBM and CBM states and photoinduced strain resulting from latticerelaxation after photoexcited carrier generation. b Simulated atomic displacementsover time generated by the photoacoustic strain pulse propagation through the12 nm t-hBN/4BP/b-hBN heterostructure. c Experimental ΔT/T signals and simu-lated real-time phosphorus atomic displacements obtained for two heterostructureswith different t-hBN thicknesses of 12 nm and 29 nm.https://doi.org/10.1038/s41699-024-00475-8 Articlenpj 2D Materials and Applications |            (2024) 8:39 4the strainpulse linearly correlatewith the t-hBN layer thickness, e.g., thefirstecho 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 strainpulses may provide a nondestructive tool for inspecting nanosized vdWheterostructures of several atomic layer thickness.DiscussionTime-resolved measurement of photoacoustic pulse propagation in vdWb-hBN/BP/t-hBN heterostructures was conducted in tandem withfrequency-dependent time-of-flight (TOF) analysis of the propagationand echoes through the t-hBN layer to determine the full group velocitydispersion (GVD) of LA phonons along the stacking direction. Thematerial-specific optical response of the probe transmittance, being sen-sitive to the strainwaveoriginating from theBP, enabled observation of thearrival of the strain wave back at the BP after round-trip travel within theheterostructure. TheGVDof the out-of-plane LA acoustic wavemodewasobtained from the phonon transit time as well as the energy separation inthe frequency comb of the Fourier-transformed spectrum. The reliabilityof the phonon dispersion data deduced from the photoacoustics wasconfirmed by comparison with linear chain model simulations, in whichthe hBN layerswere simplified as a linearmass chain connected by springs.We also revealed the microscopic origin of photoinduced strain pulsegeneration and its propagation from first principles calculation. The time-resolved THz acoustic wave analysis demonstrated in this work mayprovide a nondestructive tool for investigating the hidden interfaces andlattice properties of nanoscale thinfilmheterostructures. It is expected thatour method utilizing photoexcitation within a several atomic layer thicksensing layer can be extended to various nanoscale systems, including vander Waals heterostructures and semiconductor quantum wells.MethodsSample preparationEncapsulated BP heterostructures were prepared in an argon-filled glovebox (oxygen and water content below 0.5 ppm) by sequentially transferringbottom hBN (b-hBN), BP, and top hBN (t-hBN) flakes onto quartz sub-strateswith awell-establishedGel-Pakdry transfermethod.Two t-hBN/BP/b-hBN samples were fabricated with thicknesses of 12 nm/1.3 nm/30 nm(sampleA) and 29 nm/2.3 nm/27 nm(sampleB). The thickness of hBNwasestimated from thin film interference in a reflectance spectrum, whereas thethickness of BP was measured from the surface morphology determined byatomic force microscopy (AFM).Pump–probe measurementTo study the generation andpropagation of coherent acoustic phonons in thehBN/BP/hBN heterostructures, a pump–probe experiment was performedusing ultrashort pulses of 60 fs duration supplied by a 1MHz repetition-ratecavity-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 pumppulse and aweak probe pulse (Fig. 1a). Because hBN layers are transparent inthe near-infrared, the pump pulse excites carriers and coherent phononsselectively 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 delayrelative to the pump pulse. To maximize the light–matter interaction inanisotropic BP, the polarization of both pump and probe beamswas set to beparallel to the armchair crystal axis of BP. Using a 20× objective lens, 1 nJpump and 0.5 nJ probe pulses overlapped in a 5 μm diameter spot in thesample. The time delay between the pump and probe pulses was scanned at13Hz, and the signal was acquired until the noise was sufficiently suppressedto obtain a clear oscillation signal of coherent acoustic phonons.Linear chain model (LCM) for the out-of-plane longitudinalacoustic (LA) modeBased on symmetry considerations, thermal expansion of the BP layer uponphotoexcitation was expected to generate a strong out-of-plane LA mode.Themacroscopic lattice behavior of the hBN layerswas simplifiedwith a 1DLCM, where each layer was modeled as a mass particle and the interlayerinteraction as a spring force39. According to the LCM, the phonon disper-sion is given by ω(q) = ω0 × sin(q/2c), where ω0, q, and c are the maximumangular frequency at the Brillouin zone edge, the wavevector, and theinterlayer spacing of hBN (c = 3.33 Å), respectively. The group velocity, vg, iscalculated as vg(q) = v0 × cos(q/2c), where v0 =ω0/2c is the sound velocity inthe long wavelength limit. The strain wave echo time follows the relationTe = 2d/ vg, where the effective t-hBN thickness d is used to include the effectof finite BP thickness. The value of v0 was assumed to be v0 = 3.44 km/sbased on literature from inelastic X-ray scattering40.Electronic band structure calculationAb initio phonon calculations were carried out for bulk hBN, few-layer BP,and their vertically stacked structures using DFT with the optB86b-vdWexchange-correlation functional as implemented in the VASP package41–44.Phonon dispersions, harmonic eigenmodes, and phonon group velocitieswere calculated using the phonopy code45. The optB86b-vdW functional46,47was used to calculate both the finite band gap and group velocities tocompare with experimental results. Benchmark calculations of exchange-correlation functionals for predictions of electronic structures and groupvelocities are provided in the Supplementary Information.Photoinduced strain calculationTo model the light-induced strain in the hBN/BP/hBN heterostructure,atomic displacements inducedby carrier excitationwere calculated using anoccupation-constrainedΔSCF approach for a slabmodel 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 ofthe valence and conduction bands, as illustrated in Fig. 4b48–52. The excitedcarrier imposed by the occupation constraint was estimated as6.36 × 1012 cm−2, which is a moderate photoexcitation density49,53,54.Machine-learned force field (MLFF) constructionAnMLFFwas trained by on-the-flymolecular dynamics (MD) simulationsof an isothermal-isobaric molecular ensemble at 300 K for 100,000 stepswith a 2 ps time step using the Langevin thermostat55,56, as implemented inthe 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 theliterature57–61. The present study employed cutoffs of 5 Å and 8 Å for theradial and angular descriptors of the Gaussian approximation potential62,63,respectively, and aweighting factor of 1000was imposedon the ionic forces.Validation of the MLFF-predicted results compared to ab initio counter-parts is presented in Supplementary Fig. 4.Photoacoustic molecular dynamics simulation based on anMLFF approachAfter calculating the light-induced initial strain, the real-timepropagationofthe strain pulse was simulated via MD simulation using the MLFF, asimplemented in the on-the-fly machine learning routine in the VASPpackage57–59. Atomic trajectories were integrated with the Verlet algorithmunder an energy conservation constraint with a time step of 0.5 ps and totalsimulation time of 50 ps (total 100,000 steps). The atomic displacements ateach time stepwere convolutedwith a Lorentzian functionwith a linewidthof 1 Å for visualization of the wave packets of the strain pulse. Vibrationalpower spectra of the BP and hBN layers were obtained using the Dyna-PhoPy package by calculating the FT spectra of velocity autocorrelationfunctions64. For theMD simulations, an absorbing boundary condition wasused for the bottom hBN layer.Data availabilityAll relevant data are available from the authors upon reasonable request.Received: 23 December 2023; Accepted: 30 May 2024;https://doi.org/10.1038/s41699-024-00475-8 Articlenpj 2D Materials and Applications |            (2024) 8:39 5References1. Slack, G. A. Nonmetallic crystals with high thermal conductivity. J.Phys. Chem. Solids 34, 321–335 (1973).2. Kittel, C. Introduction to solid state physics, 8th edition. (Wiley, 2004).3. Shaw, W. M. & Muhlestein, L. D. Investigation of the phonondispersion relations of chromium by inelastic neutron scattering.Phys. Rev. B 4, 969–973 (1971).4. Burkel, E. Phonon spectroscopy by inelastic x-ray scattering. Rep.Prog. Phys. 63, 171–232 (2000).5. Wright, O. B., Hyoguchi, T., Hyoguchi, T., Kawashima, K. &Kawashima, K. Laser picosecond acoustics in thin films: effect ofelastic boundary conditions on pulse generation. Jpn. J. Appl. Phys.30, L131 (1991).6. Thomsen, C., Grahn, H. T., Maris, H. J. & Tauc, J. Surface generationanddetectionof phononsbypicosecond light pulses.Phys.Rev.B34,4129–4138 (1986).7. Chern, G.-W., Lin, K.-H. & Sun, C.-K. Transmission of light throughquantumheterostructuresmodulated by coherent acoustic phonons.J. Appl. Phys. 95, 1114–1121 (2004).8. Hurley, D. H., Lewis, R., Wright, O. B. & Matsuda, O. Coherent controlof gigahertz surface acoustic and bulk phonons using ultrafast opticalpulses. Appl. Phys. Lett. 93, 113101 (2008).9. Profunser, D. M., Wright, O. B. & Matsuda, O. Imaging ripples onphononic crystals reveals acoustic band structure and Blochharmonics. Phys. Rev. Lett. 97, 055502 (2006).10. Otsuka, P. H. et al. Time-domain imaging of gigahertz surface waveson an acoustic metamaterial. New J. Phys. 20, 013026 (2018).11. Ajayan, P., Kim, P. & Banerjee, K. Two-dimensional van der Waalsmaterials. Phys. Today 69, 38–44 (2016).12. Geim, A. K. & Grigorieva, I. V. Van der Waals heterostructures.Nature499, 419–425 (2013).13. Mishina, T., Nitta, K. & Masumoto, Y. Coherent lattice vibration ofinterlayer shearing mode of graphite. Phys. Rev. B 62,2908–2911 (2000).14. Bartram, F. M. et al. Ultrafast coherent interlayer phonon dynamics inatomically thin layers of MnBi2Te4. npj Quantum Mater. 7, 84 (2022).15. Bae, S. et al. K-point longitudinal acoustic phonons are responsiblefor ultrafast intervalley scattering inmonolayerMoSe2.Nat. Commun.13, 4279 (2022).16. Zhang, M. Y. et al. Light-induced subpicosecond lattice symmetryswitch in MoTe2. Phys. Rev. X 9, 021036 (2019).17. Vialla, F. & Del Fatti, N. Time-domain investigations of coherentphonons in van der Waals thin films. Nanomaterials 10, 2543 (2020).18. Jeong, T. Y. et al. Coherent lattice vibrations in mono- and few-layerWSe2. ACS Nano 10, 5560–5566 (2016).19. Ge, S. et al. Coherent longitudinal acoustic phonon approaching THzfrequency in multilayer molybdenum disulphide. Sci. Rep. 4,5722 (2014).20. Balandin, A. A. Nanophononics: phonon engineering innanostructures and nanodevices. J. Nanosci. Nanotechnol. 5,1015–1022 (2005).21. Huynh, A., Perrin, B. & Lemaitre, A. Semiconductor superlattices: atool for terahertz acoustics. Ultrasonics 56, 66–79 (2015).22. Yang, C.-Y., Wu, P.-C., Chu, Y.-H. & Lin, K.-H. Generation andcoherent control of terahertz acoustic phonons in superlattices ofperovskite oxides. New J. Phys. 23, 053009 (2021).23. Pascual Winter, M. F. et al. Selective optical generation of coherentacoustic nanocavity modes. Phys. Rev. Lett. 98, 265501 (2007).24. Cai, Y., Zhang, G. & Zhang, Y. W. Layer-dependent band alignmentandwork function of few-layer phosphorene.Sci. Rep. 4, 6677 (2014).25. Cassabois, G., Valvin, P. & Gil, B. Hexagonal boron nitride is anindirect bandgap semiconductor.Nat. Photonics 10, 262–266 (2016).26. Xue, J. et al. Scanning tunnelling microscopy and spectroscopy ofultra-flat graphene on hexagonal boron nitride. Nat. Mater. 10,282–285 (2011).27. Qiu, D. Y., da Jornada, F. H. & Louie, S. G. Environmental screeningeffects in 2D materials: renormalization of the bandgap, electronicstructure, and optical spectra of few-layer black phosphorus. NanoLett. 17, 4706–4712 (2017).28. Suess, R. J., Jadidi, M. M., Murphy, T. E. & Mittendorff, M. Carrierdynamics and transient photobleaching in thin layers of blackphosphorus. Appl. Phys. Lett. 107, 081103 (2015).29. Miao, X., Zhang, G., Wang, F., Yan, H. & Ji, M. Layer-dependentultrafast carrier and coherent phonon dynamics in black phosphorus.Nano Lett. 18, 3053–3059 (2018).30. Meng, S., Shi, H., Jiang, H., Sun, X. & Gao, B. Anisotropic chargecarrier and coherent acoustic phonon dynamics of black phosphorusstudied by transient absorption microscopy. J. Phys. Chem. C. 123,20051–20058 (2019).31. Lee, S. Y. & Yee, K.-J. Anisotropic Generation and Detection ofCoherent Ag Phonons in Black Phosphorus. Nanomaterials 11,1202 (2021).32. Wu, S. et al. Dichroic photoelasticity in black phosphorus revealed byultrafast coherent phonon dynamics. J. Phys. Chem. Lett. 12,5871–5878 (2021).33. Volpato, A. & Collini, E. Time-frequency methods for coherentspectroscopy. Opt. Express 23, 20040–20050 (2015).34. Kim, Y.-H. et al. Two-dimensional limit of exchange-correlation energyfunctional approximations. Phys. Rev. B 61, 5202–5211 (2000).35. Tran, V., Soklaski, R., Liang, Y. & Yang, L. Layer-controlled band gapand anisotropic excitons in few-layer black phosphorus.Phys. Rev. B89 (2014).36. Wang, X., Gao, W. & Zhao, J. Strain modulation of the excitonanisotropy and carrier lifetime in black phosphorene. Phys. Chem.Chem. Phys. 24, 10860–10868 (2022).37. Zhang, Z. et al. Strain-modulated bandgap and piezo-resistive effectin black phosphorus field-effect transistors. Nano Lett. 17,6097–6103 (2017).38. Kim, H. et al. Actively variable-spectrum optoelectronics with blackphosphorus. Nature 596, 232–237 (2021).39. Lüerßen, D. et al. A demonstration of phonons that implements thelinear theory. Am. J. Phys. 72, 197–202 (2004).40. Bosak, A. et al. Elasticity of hexagonal boron nitride: Inelastic x-rayscattering measurements. Phys. Rev. B 73, 041402(R) (2006).41. Kresse, G. & Hafner, J. Ab initio molecular dynamics for liquid metals.Phys. Rev. B 47, 558–561 (1993).42. Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energycalculations formetals and semiconductors using a plane-wavebasisset. Comput. Mater. Sci. 6, 15–50 (1996).43. Kresse, G. & Furthmuller, J. Efficient iterative schemes for ab initiototal-energy calculations using a plane-wave basis set. Phys. Rev. B54, 11169–11186 (1996).44. Kresse, G. & Joubert, D. Fromultrasoft pseudopotentials to the projectoraugmented-wave method. Phys. Rev. B 59, 1758–1775 (1999).45. Togo, A. & Tanaka, I. First principles phonon calculations in materialsscience. Scr. Mater. 108, 1–5 (2015).46. Klimes, J., Bowler, D. R. & Michaelides, A. Chemical accuracy for thevan der Waals. density Funct. J. Phys. Condens. Matter 22,022201 (2010).47. Graziano, G., Klimes, J., Fernandez-Alonso, F. & Michaelides, A.Improved description of soft layered materials with van der Waals.density Funct. theory J. Phys. Condens. Matter 24, 424216 (2012).48. Görling, A. Density-functional theory beyond the Hohenberg-Kohntheorem. Phys. Rev. A 59, 3359–3374 (1999).49. Krishnamoorthy, A. et al. Semiconductor-metal structural phasetransformation in MoTe2 monolayers by electronic excitation.Nanoscale 10, 2742–2747 (2018).50. Lee, J., Kim, H. S. & Kim, Y. H.Multi-space excitation as an alternativeto the landauer picture for nonequilibrium quantum transport. Adv.Sci. 7, 2001038 (2020).https://doi.org/10.1038/s41699-024-00475-8 Articlenpj 2D Materials and Applications |            (2024) 8:39 651. Lee, J., Yeo, H. & Kim, Y. H. Quasi-Fermi level splitting in nanoscalejunctions from ab initio. Proc. Natl Acad. Sci. USA 117,10142–10148 (2020).52. Bae, S., Jeong, T. Y., Raebiger, H., Yee, K. J. & Kim, Y. H. Localizedcoherent phonon generation in monolayer MoSe2 from ultrafastexciton trapping at shallow traps. Nanoscale Horiz. 8,1282–1287 (2023).53. Liu, F., Ziffer, M. E., Hansen, K. R., Wang, J. & Zhu, X. Directdetermination of band-gap renormalization in the photoexcitedmonolayer MoS2. Phys. Rev. Lett. 122, 246803 (2019).54. Trovatello, C. et al. Strongly Coupled Coherent Phonons in Single-Layer MoS2. ACS Nano 14, 5700–5710 (2020).55. Parrinello, M. & Rahman, A. Crystal structure and pair potentials: amolecular-dynamics study. Phys. Rev. Lett. 45, 1196–1199 (1980).56. Parrinello, M. & Rahman, A. Polymorphic transitions in single crystals:A new molecular dynamics method. J. Appl. Phys. 52,7182–7190 (1981).57. Jinnouchi, R., Karsai, F. & Kresse, G. On-the-fly machine learningforcefieldgeneration: Application tomeltingpoints.Phys.Rev.B.100,014105 (2019).58. Jinnouchi, R., Lahnsteiner, J., Karsai, F., Kresse, G. & Bokdam, M.Phase transitions of hybrid perovskites simulated by machine-learning force fields trained on the fly with bayesian inference. Phys.Rev. Lett. 122, 225701 (2019).59. Jinnouchi, R., Karsai, F., Verdi, C., Asahi, R. & Kresse, G. Descriptorsrepresenting two- and three-body atomic distributions and theireffects on the accuracy ofmachine-learned inter-atomic potentials. J.Chem. Phys. 152, 234102 (2020).60. Thiemann, F. L., Rowe, P., Müller, E. A. & Michaelides, A. MachineLearning Potential for Hexagonal Boron Nitride Applied to Thermallyand Mechanically Induced Rippling. J. Phys. Chem. C. 124,22278–22290 (2020).61. Bokdam, M., Lahnsteiner, J. & Sarma, D. D. Exploring librationalpathways with on-the-Fly machine-learning force fields:methylammonium molecules in MAPbX3 (X = I, Br, Cl) perovskites. J.Phys. Chem. C. 125, 21077–21086 (2021).62. Bartok, A. P., Payne, M. C., Kondor, R. & Csanyi, G. Gaussianapproximation potentials: the accuracy of quantum mechanics,without the electrons. Phys. Rev. Lett. 104, 136403 (2010).63. Bartók, A. P., Kondor, R. & Csányi, G. On representing chemicalenvironments. Phys. Rev. B 87, 184115 (2013).64. Carreras, A., Togo, A. & Tanaka, I. DynaPhoPy: A code for extractingphonon quasiparticles from molecular dynamics simulations.Computer Phys. Commun. 221, 221–234 (2017).AcknowledgementsThis research was funded by the National Research Foundation of Korea(NRF-2020R1A6A1A03047771, NRF-2023R1A2C1004284). S.B. wassupported by Grant-in-Aid for JSPS Fellows (No. 23KF0030) and H.R. (No.23K04356) from the Japan Society for the Promotion of Science. Intensivecalculations were conducted on Fujitsu PRIMERGY servers CX400M1/CX2550M5 (Oakbridge-CX) at the Information Technology Center and ISSPsupercomputers at the Institute for Solid State Physics, University of Tokyo.Author contributionsS.B. andK.-J.Y. conceived and coordinated this project. S.-Y. L. andK.-J.Y.constructed the sub-ps pump–probe setup andperformed the experiments,S.B. andH.R.performedfirst-principlescalculations, andS.-Y.K., S.J., K.W.,and T.T. synthesized and characterized the van derWaals heterostructures.S.B. and K.-J.Y. wrote the manuscript with contributions from all authors.Competing interestsThe authors declare no competing interests.Additional informationSupplementary information The online version containssupplementary material available athttps://doi.org/10.1038/s41699-024-00475-8.Correspondence and requests for materials should be addressed toSoungmin Bae or Ki-Ju Yee.Reprints and permissions information is available athttp://www.nature.com/reprintsPublisher’s note Springer Nature remains neutral with regard tojurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative CommonsAttribution 4.0 International License, which permits use, sharing,adaptation, distribution and reproduction in anymedium or format, as longas you give appropriate credit to the original author(s) and the source,provide a link to the Creative Commons licence, and indicate if changeswere made. The images or other third party material in this article areincluded in the article’s Creative Commons licence, unless indicatedotherwise in a credit line to the material. If material is not included in thearticle’sCreativeCommons licence and your intended use is not permittedby statutory regulation or exceeds the permitted use, you will need toobtain permission directly from the copyright holder. To view a copy of thislicence, visit http://creativecommons.org/licenses/by/4.0/.© The Author(s) 2024https://doi.org/10.1038/s41699-024-00475-8 Articlenpj 2D Materials and Applications |            (2024) 8:39 7https://doi.org/10.1038/s41699-024-00475-8http://www.nature.com/reprintshttp://creativecommons.org/licenses/by/4.0/ Full phonon dispersion along the stacking direction in nanoscale van der Waals materials by picosecond acoustics Results Optical pump–probe experiment First principles simulation Discussion Methods Sample preparation Pump–probe measurement Linear chain model (LCM) for the out-of-plane longitudinal acoustic (LA)�mode Electronic band structure calculation Photoinduced strain calculation Machine-learned force field (MLFF) construction Photoacoustic molecular dynamics simulation based on an MLFF approach Data availability References Acknowledgements Author contributions Competing interests Additional information