Ultrafast coherent interlayer phonon dynamics in atomically thin layers of MnBi2Te4

The atomically thin MnBi2Te4 crystal is a novel magnetic topological insulator, exhibiting exotic quantum physics. Here we report a systematic investigation of ultrafast carrier dynamics and coherent interlayer phonons in few-layer MnBi2Te4 as a function of layer number using time-resolved pump-probe reflectivity spectroscopy. Pronounced coherent phonon oscillations from the interlayer breathing mode are directly observed in the time domain. We find that the coherent oscillation frequency, the photocarrier and coherent phonon decay rates all depend sensitively on the sample thickness. The time-resolved measurements are complemented by ultralow-frequency Raman spectroscopy measurements, which both confirm the interlayer breathing mode and additionally enable observation of the interlayer shear mode. The layer dependence of these modes allows us to extract both the out-of-plane and in-plane interlayer force constants. Our studies not only reveal the interlayer van der Waals coupling strengths, but also shed light on the ultrafast optical properties of this novel two-dimensional material.


INTRODUCTION
The marriage between topology and magnetism can give birth to many exotic quantum phases and phenomena such as quantum anomalous Hall effect [1][2][3][4] and axion electrodynamics [5][6][7][8] . Recently, few-layer MnBi 2 Te 4 (MBT) crystals have emerged as a new platform for exploring these phases  . MBT is a new class of two-dimensional (2D) topological magnetic material. The crystal structure of MBT is formed by Te-Bi-Te-Mn-Te-Bi-Te septuple layers (SLs) stacked on top of each other along the c direction (Fig. 1a). The topological properties are inherited from the p orbitals of Bi and Te atoms, while the magnetism originates from the d bands of Mn atoms. The adjacent SLs are bonded by van der Waals (vdW) forces and can be cleaved to thin layers. Both electronic and magnetic properties depend sensitively on the sample thickness. Transport studies have revealed rich quantized transport phenomena in few-layer MBT (5-9 SL), suggesting odd-and evenlayer MBT samples may host the quantum anomalous Hall and axion insulator states respectively [24][25][26][27] . Magnetic circular dichroism measurements have revealed the layer-dependent magnetism with an odd-even layer-number effect, verifying the A-type antiferromagnetic structure in few-layer MBT 28,29 .
In 2D vdW materials, lattice vibrations, especially interlayer phonon modes, have not only provided unique capabilities in determining layer thickness, stacking order, and interface coupling strength [34][35][36] , but have also played an important role in engineering novel electronic, thermal and magnetic properties in 2D vdW homo/heterostructures [37][38][39][40][41][42] . For instance, in twisted bilayer graphene interlayer conductance and thermal transport are mediated by the interlayer phonons [39][40][41][42] . Several Raman studies have been performed recently, primarily on highfrequency intralayer phonons which have demonstrated spinphonon coupling [31][32][33] . In the low-frequency regime, magnon modes have been observed 43 , and a systematic study of interlayer breathing mode phonons showed a very large interfacial coupling which produced somewhat unusual results 44 . These interlayer phonons provide information about vdW coupling strength and electron-phonon scattering mechanisms, and a good understanding of their properties is useful for future work on MBT, for example in manipulating magnetic and topological states with light at a fast speed 23,45,46 , in the design and probe of heterostructure devices 37,38 , or simply as a method for determining sample thicknesses 35,36 .
Here we report ultrafast optical pump-probe reflectivity measurements in few-layer MBT, accompanied by ultralowfrequency Raman data, of carrier and coherent interlayer phonon dynamics as a function of sample thickness, with the layer number varying from 4 to 25. We observe pronounced oscillatory signals due to coherent interlayer breathing mode phonons in transient reflection measurements and find that physical properties such as the oscillation frequency, its decay time, and the carrier decay rate all depend strongly on the sample thickness. The frequency is particularly consistent and easily measured, which provides a convenient, contact-free and non-destructive tool to characterize sample thickness, especially for ultrafast studies. The increasing decay rates for thinner samples suggest that the surface and interface have a strong impact on the phonon scattering and photocarrier relaxation. The time-domain measurements are complemented by ultralow-frequency (<10 cm −1 ) Raman studies. In addition to the interlayer breathing modes, polarized Raman measurements enable the observation of the much weaker interlayer shear modes. The measured phonon frequencies all show a typical blue shift with decreasing layer number and can be fit using a linear chain model, from which the out-of-plane and in-plane interlayer force constants are calculated to be (5.7 ± 0.1) Ã 10 19 N m −3 and (2.9 ± 0.1) Ã 10 19 N m −3 respectively. The force constants allow us to calculate other mechanical parameters such as sound velocities, elastic constants and acoustic impedance. These measurements not only uncover the interlayer coupling strengths but also provide crucial information on coherent phonon and carrier lifetimes and relaxation mechanisms.

RESULTS
Atomic force microscopy measurements Few-layer MBT crystals were prepared by mechanical exfoliation as described previously 25 . Figure 1b shows an optical image of one of the sample regions used in this work. The layer thickness was determined by atomic force microscopy (AFM) measurements ( Fig. 1c-e), where clear step sizes of around 1.4 nm can be observed, consistent with previous studies [24][25][26]29 .

Time-resolved reflectivity measurements
We measured the transient reflectivity with a standard two-color pump-probe setup, the results of which can be seen in Fig. 2. In the raw data, shown in Fig. 2a, we clearly see the presence of three distinct components. First, there is a large but short-lived change which lasts only 2-3 picoseconds with negligible layer dependence before vanishing, associated with hot electrons which rapidly transfer energy to the surrounding lattice via electronelectron and electron-phonon interactions. Next, we see a simple exponential decay, with a lifetime of around 100 ps depending on the layer thickness, probably reflecting photocarrier scattering and electron-hole recombination processes. The final component is an oscillatory signal which decays over~30-50 picoseconds, and then briefly reappears at t ≈ 100 ps (see inset) before decaying away once more. Note that the oscillatory signal is very strong, about 20% of the total transient reflection. To show the oscillatory signal more clearly, the data after the first two components are subtracted is plotted again in Fig. 2b, along with the FFT of the data in Fig. 2c, where you can see a clear upward shift in frequency as the sample thickness decreases from~130 GHz in 9 SL to~180 GHz in 7 SL. The inset shows the FFT of the echoes, which have frequencies that are identical to the main signals.
In a pump-probe experiment, coherent phonons created by the pump pulse can modulate the optical dielectric function and can therefore be detected via reflectivity changes measured by the synchronized probe pulse. These phonons can be created via rapid expansions or contractions induced by, for example, temperature changes 47,48 or carrier excitation 49 , or through a stimulated Raman scattering process 50,51 . Because the observed oscillation frequencies are much smaller than the intralayer vibration frequencies measured by Raman spectroscopy 31 (also see discussions below), we attribute these results to coherent interlayer phonons, which for ultrathin samples will form standing waves with strongly layerdependent frequencies. In our measurements, the polarization of the pump beam does not affect the observed signal (see Supplementary Note 3). This suggests that the interlayer vibrations are in the out-of-plane direction (the breathing mode), as in-plane vibrations (the shear mode) would show a phase shift when the pump polarization is changed 51,52 , a fact which has been utilized to detect very weak coherent interlayer shear mode phonons in graphene 52,53 . The assignment of the coherent phonon mode is further confirmed by our polarization-resolved ultralow-frequency Raman spectroscopy measurements as discussed below.
This mechanism also explains the echo of the signal at 100 ps, as coupling between the sample and the~300 nm thick silica layer on the substrate means that acoustic phonons can be generated, travel away from the surface, reflect from the internal silica/Si interface, and return to the surface, where they then re-excite the interlayer phonon in the sample. The timing of the reappearance gives a velocity of about 6 km s −1 , which matches well with the speed of longitudinal acoustic phonons in silica. The echo signals indicate that few layer MBT can be used as an optical transducer to launch coherent acoustic phonons in SiO 2 .
In Fig. 3a, we plot the extracted oscillation frequency at a variety of sample locations for a given layer thickness as a function of the inverse layer number (1/N). The frequency decreases with increasing sample thickness from~300 GHz in 4 SL to~50 GHz in 25 SL samples.
We analyze how the frequency of these phonons vary with thickness using a basic freestanding linear chain model [34][35][36] , in which each SL is moving as a unit and only the nearest layer interaction is considered (see Supplementary Note 9). A similar model has been used to explain the interlayer vibration modes in a variety of other 2D layered materials [34][35][36][54][55][56][57][58][59] . An N-layer system hosts N-1 breathing/shear modes with frequencies f(N,n) = f 0 sin (nπ⁄(2N)), where n = 1, 2, ⋯, N−1 is the branch index and is the out-of-plane (or inplane) force constant per unit area between the adjacent layers and μ is the mass per unit area. Fitting our results to f(N, 1) gives a resulting fit parameter f 0,B = 760 ± 5 GHz for the breathing mode, which we plot as a dashed gray line in Fig. 3a, showing excellent agreement of this theoretical model with the experimental data. Note that when the number of layers N→∞, these modes correspond to longitudinal acoustic phonons along the out-of- , with d the spacing between layers. Therefore, the longitudinal sound velocity is given by v B = πdf 0 = 3.25 km s −1 . This value is comparable to v B in other 2D layered materials such as MoS 2 (see Supplementary Table 1).
One obvious advantage of time-resolved experiments is that we can directly measure the phonon and photocarrier lifetimes. We plot the decay rate of both the oscillatory component and the simple exponential decay in Fig. 3b, c respectively. These decay rates noticeably increase for the thinner regions of the sample. For a given thickness, we performed multiple measurements on different locations. Both the phonon decay and carrier relaxation are very sensitive to the local environment especially for the thinner samples shown as scattered data points for the same thickness in Fig. 3b, c, suggesting that phonon and carrier relaxation is influenced strongly by the surface and interface (e.g., surface roughness, adatoms and charges trapped at the interface, which all act as scatterers). On the other hand, the coherent phonon frequency is measured to be robust against locations (Fig. 3a).
In general, coherent phonon dephasing in solids is governed by scattering at boundaries, lattice defects, carriers or population decaying into lower energy acoustic phonons via anharmonic interactions. As shown in Fig. 3b, for thicker samples (>10 layers) the phonon decay rate is~50 GHz and is independent of the thickness and sample locations. By contrast, in the thin layer regime (<10 layers) the decay rate grows rapidly with decreasing thickness and depends strongly on the location, indicating that the surface and interface play a dominant role. As demonstrated by the echo of the signal, the interlayer phonons can couple to acoustic phonons in the substrate, which may further increase the interfacial contribution to the decay rate. Furthermore, the layer dependence of phonon decay rate can also result from stronger electron-phonon coupling in thinner samples than in thicker ones 31 , similar to topological materials such as Bi 2 Se 3 and WP 2 60,61 . Additional insight into the anharmonic contribution to the phonon decay rate can be obtained from temperaturedependent data ( Supplementary Fig. 5). If the anharmonic contribution were large, the decay rate (linewidth) would increase significantly upon increasing temperature 62 . Our measurement shows that the decay rate stays roughly constant with temperature from 3 K to room temperature, suggesting that the anharmonic interaction is not important. We have found that the transient reflectivity simply scales linearly with the pump fluence (see Supplementary Note 6), which eliminates scattering with photoinduced carriers as a prominent dephasing factor.
The photoinduced carrier relaxation is determined by carrier scatterings with defects and lattice and electron-hole recombination, etc. The increasing decay rate in thinner samples from 2 GHz in 20 SL to~18 GHz in 4 SL shown in Fig. 3c indicates accelerated scattering and recombination rates at the surface and interface compared to the inner layers. The approximate linear scaling of the decay rate with (1/N) can be understood with a simple model, in which the total decay rate is a weighted average between the inner layer decay rate and surface/interface decay rate. Similar behavior has been observed in 2D MoS 2 and Bi 2 Te 3 63,64 .

Raman spectroscopy measurements
To confirm our observations, we also performed polarized Raman measurements on our sample, shown in Fig. 4a, b. The six highfrequency peaks (at~27, 48, 68, 105, 115 and 141 cm −1 ) have been observed in previous studies 31 and are intralayer phonon modes that are Raman active both in bulk and few-layer samples. We assign these modes corresponding to their symmetry class 22 based on their polarization dependence (see Supplementary Note 1 for more details). The out-of-plane vibrations (A 1g modes) can be observed only in a parallel-polarized geometry, whereas the inplane vibrations (E g modes) are detectable in both parallel-and cross-polarized geometries. At low frequencies, we also observe new peaks which we attribute to the interlayer phonons. In addition to the breathing mode seen in the time-resolved data, (labeled B (1) ), Raman measurements also enable observation of the much weaker shear mode (labeled S (1) ). Note that unlike the breathing modes, the shear modes are sensitive to applied in-plane strain 65 , which potentially can be used as a strain probe in MBT and heterostructures when engineering electronic and magnetic properties with strain. The B (1) peak intensity is remarkably strong compared to all other peaks, which is consistent with the pronounced coherent oscillations in the time-resolved measurements as the coherent phonons are likely generated via the stimulated Raman scattering process. The frequencies of the B (1) and S (1) modes are plotted in Fig. 4c alongside the time-resolved data, showing the excellent agreement of the B (1) frequencies with those observed in the transient reflectivity. In addition, we show that the S (1) modes also fit well to the linear chain model, with a slightly lower f 0,S = 545 ± 7 GHz. The linewidths of these modes (plotted in Supplementary Fig. 2b) are very noisy, making any layer dependence difficult to observe, but are reasonable compared to the decay rates in the time-resolved data. While some small additional peaks can be observed in the parallelpolarized data (e.g., around~20 cm −1 in Fig. 4a), they have been left intentionally unlabeled as their precise origin is unclear and is not the focus of this work. Note that only the B (1) mode is seen in the time-resolved measurements, but since Raman scattering intensity is correlated with the amplitude of coherent phonons that can be generated 51 , these other modes may simply be too weak to be observed.
Since both the B (1) and S (1) modes show a large blue shift from 120 GHz in 10 SL to~300 GHz in 4 SL and from~100 GHz in 9 SL to~200 GHz in 4 SL, respectively, the frequency dispersions of these interlayer modes could be useful for determining the thickness of few-layer MBT. By contrast, among the highfrequency intralayer modes, only the A 1g (1) mode has a noticeable change 56 , red-shifting from~47.5 to~44.5 cm −1 (a difference of 90 GHz) from 10 SL to 4 SL (see Supplementary Note 2), similar to previously reported observations 31 . We also point out that while in Raman spectroscopy these interlayer mode frequencies are very low and may be difficult to measure, especially for thick samples (where for our setup, above~13 SL the frequencies drop below the detection limit), time-resolved measurements can be used to quickly and accurately obtain the breathing mode frequency across a wide thickness range. Acquiring time-resolved data could therefore be a particularly effective tool for non-destructive measurements of MBT sample thickness.

DISCUSSION
As mentioned previously, we have found that for the breathing mode, f 0,B = 760 ± 5 GHz from the time-resolved data and for the shear mode, f 0,S = 545 ± 7 GHz from the S (1) Raman peaks. Using the B (1) Raman peaks to determine the breathing mode frequency results in a slightly higher and more uncertain f R 0;B ¼ 780 ± 15 GHz, likely due to the large Rayleigh scattering background in parallelpolarized Raman measurements at such low frequencies. Using a mass per unit area per layer for MBT of μ = 10.0 Ã 10 −6 kg m −2 , we find the out-of-plane and in-plane force constants per unit area to be k z = (5.7 ± 0.1) Ã 10 19 N m −3 and k x = (2.9 ± 0.1) Ã 10 19 N m −3 . The out-of-plane force constant k z compares well with a recent Raman study, where they indirectly estimated k z = (4.5 ± 0.7) Ã 10 19 N m −3 from the Davydov splitting in bulk MnBi 4 Te 7 32 . Recent first-principles calculations of the interlayer phonon modes in bilayer MBT 23 give force constants of k z = 2.85 Ã10 19 N m −3 and k x = 1.66 Ã10 19 N m −3 , which like previous studies on Bi 2 Te 3 58 are much smaller than the experimental values but show a similar ratio between the out-of-plane and in-plane force constants (~1.7 from the theory, compared to~1.9 from the experiment).
As mentioned in the introduction, a very recent study also measured the interlayer breathing mode phonons using Raman spectroscopy 44 . They observed somewhat unusual results that were À Á with f 0,B = 760 GHz. The decay rates of both the oscillatory component (b) and the simple exponentially decaying component (c) also increase for thinner samples, with the latter being particularly noticeable. Note that statistical uncertainties from fitting are very small (≈0.1%, 1% and 0.5% for points in a, b and c respectively, all within marker sizes), but practical uncertainties are clearly much larger due to factors such as variation across sample locations.
F.M. Bartram et al. attributed to a strong interfacial coupling, supported by the fact that measuring samples on two different substrates gave different results. This fact necessitated the use of a more complicated linear chain model to include the effects of the substrate coupling, which they found explained their results well. Our measurements, in contrast, agree well with a simpler model where interfacial effects are neglected (see Supplementary Discussion 1 for additional discussion of this model). This is somewhat surprising, considering one of the substrates used in their study was the same as ours (SiO 2 /Si), and suggests that more work should be done to determine under what circumstances these interfacial effects will vary. However, despite this difference, their result for the out-ofplane force constant k z = (4.3-6.4) Ã 10 19 N m −3 agrees quite well with ours.
A useful number for comparison with other 2D materials is the per bottom atom rather than per unit area force constants, which are k z/atom = (9.3 ± 0.1) N m −1 and k x/atom = (4.7 ± 0.1) N m −1 . These per-atom values are smaller than those found in Bi 2 Te 3 , but much larger than graphene (see Supplementary Table 1 for details), in line with the relative difficulties of mechanical exfoliation for these different samples. In Supplementary   Table 1, we further calculate and compare our results for the interlayer phonon frequencies f 0 , acoustic velocities, elastic constants C 33 and C 44 , and acoustic impedances of MBT with a few representative vdW materials.
We also investigate the interlayer breathing phonons as a function of temperature and magnetic field (see Supplementary Note 5) in order to determine if there are any strong links between the interlayer phonons and magnetic ordering in few-layer MBT. We find that the breathing mode frequency decreases by~5% from 3 K to room temperature ( Supplementary Fig. 5). This softening results from thermal expansion, which weakens the interlayer coupling with increasing temperature. Note that both magnetic and topological electronic properties are affected by hydrostatic pressure shown by recent theoretical and experimental studies in the MBT family compounds [66][67][68] ; the interlayer phonons provide a sensitive probe of the lattice separation and external strain/pressure. We do not see a noticeable change in the breathing phonon frequency as the temperature crosses the Neel temperature at~20 K ( Supplementary Fig. 5) nor when a magnetic field of up to 6 T is applied to the sample at 3 K ( Supplementary  Fig. 6). Ref. 31 reports that there is a 0.3% frequency shift of the With parallel polarization (a), the data contains 6 peaks at high frequencies which match previously reported Raman data, as well as a very large low-frequency peak corresponding to the phonon breathing mode, labeled B (1) . The Rayleigh scattering background is subtracted. In the cross-polarized configuration (b) this large peak vanishes, and an additional peak corresponding to the shear mode, labeled S (1) , is revealed. Plotting the frequency of the B (1) and S (1) peaks as a function of thickness (c) along with the previously shown time-resolved data (purple and green diamond symbols for Samples 1 and 2 respectively) shows a good agreement of the breathing mode with the time-resolved data as well as a good fit of the shear mode to the linear chain model with f 0,S = 545 GHz. Although the statistical uncertainties remain small (≈0.1% for B peaks and ≈1% for S peaks), compared to the time-resolved data the extracted peak locations show more variation, probably due to the influence of the background, which is very large at such low wavenumbers.
F.M. Bartram et al. ~48 cm −1 A 1g (1) mode across the phase transition temperature. We note that as the frequency change due to the spin-phonon coupling is proportional to the phonon frequency itself, similar coupling between interlayer breathing phonons and the magnetic order would likely only produce changes that are too small to observe. For example, in a 6-layer sample with~200 GHz oscillations, a 0.3% change is only 0.6 GHz, which is within the detection uncertainty.
In conclusion, we have directly observed the interlayer breathing mode in few layer MBT both in the time-domain and in the frequency domain. By measuring the vibration frequency as a function of layer number, we have obtained both the out-of-plane and in-plane interlayer force constants. A good understanding of the typical vdW coupling strength and interlayer vibrational properties will set the stage for engineering electronic and magnetic properties 69 (e.g., vdW heterostructures) as a device performance may be affected through interlayer electron-phonon interactions 37,[70][71][72] . Because the energies of interlayer phonons in MBT are small (~1 meV), we expect that they play an important role in interlayer conductance even at low temperatures 39,41 . Our time-resolved measurements have also revealed that the surface and interface have a strong impact on the relaxation of both the coherent phonons and photocarriers, which may also influence carrier mobility. Therefore, a systematic study of transport properties as a function of layer number is desirable. As proposed in ref. 23 , topological and magnetic phase transitions in MBT can be manipulated via non-linear phonon dynamics using intense terahertz (THz) light pulses. The light pulses induce lattice distortions leading to the separation of the septuple layers, and thus change the magnetic coupling between layers (from AFM to FM) and topological states. This effect should be largest when driving resonantly at the interlayer breathing phonon frequencies measured in this work. A similar approach has been used to realize a topological phase transition in WTe 2 46 , in which intense THz pulses drive the interlayer shear mode and induce a structure phase transition from a non-centrosymmetric to centrosymmetric phase, which resulted in the sample changing from a Weyl semimetal to a topologically trivial metal. Our studies also pave the way for future light-driven topological and magnetic orders in few-layer topological magnetic materials 23,45,46 , and ultrafast studies of nonequilibrium magnetic dynamics 73 and nonequilibrium axion dynamics [5][6][7][8][74][75][76][77] .

Sample preparation
MnBi 2 Te 4 single crystals were grown by the direct reaction of Bi 2 Te 3 and MnTe with a ratio of 1:1 in a vacuum-sealed silica ampoule. The mixture was heated to 973 K and then cooled down to 864 K slowly. The growth method and characterization are similar to those used in previous works 18,25 . The Neel temperature and the critical magnetic field for the spin flop transition of bulk crystals are T N = 25 K and H c1 = 3.6 T (at 1.5 K) respectively 18,25 . These values are shown to be correlated with the Mn content and distribution 78 . MnBi 2 Te 4 flakes were exfoliated onto 285 nm thick SiO 2 /Si substrates treated by air plasma. The fabrication processes were carried out in an argon-filled glove box. The thickness is determined by the optical contrast and AFM measurement. For optical measurements, we mounted the sample in a vacuum chamber and immediately pumped down to minimize exposure to ambient environment.

Time-resolved reflectivity spectroscopy
Time-resolved measurements were performed using a typical two-color pump-probe setup, where one beam is the direct output of a Ti:Sapphire oscillator (repetition rate: 80 MHz), and the second beam is split off and sent through an OPO to modify the wavelength. After passing through optics and a delay stage, the beams are merged with a dichroic mirror and focused onto the sample using an NA = 0.5 aspheric lens, giving a 1-2 μm spot size for the probe beam and a~5 μm spot for the pump, which is intentionally defocused slightly to allow better overlap with the probe.
The reflected light is passed through a filter to eliminate any pump scattering and measured with a balanced photodiode. We modulate the pump beam at 377 kHz with an electro-optical modulator, which allows sensitive lock-in detection of the differential reflectivity. The data in this work was primarily taken using either a 633 nm probe and an 805 nm pump or an 805 nm probe and a 1050 nm pump (see Supplementary Note 4 for discussion of other combinations). The probe and pump powers were set to~0.1 and 1 mW (pump fluence~50 μJ/cm 2 ) average powers respectively to avoid sample damage and laser heating effects. Most measurements were performed at room temperature unless otherwise stated. Temperature and magnetic field dependent measurements were carried out in an optical superconducting magnet system with the magnetic field applied perpendicular to the sample plane.

Raman spectroscopy
The Raman spectroscopy measurements were carried out in backscattering geometry at room temperature using a micro-Raman spectrometer with a charge-coupled device. The laser wavelength was 488 nm and the laser plasma lines were filtered by a reflecting Bragg grating, while the Rayleigh line is suppressed using BragGrate notch filters, achieving a detection limit of~3 cm −1 . The sample was mounted in a vacuum chamber. The laser was focused onto the sample via a 50× microscope objective lens to a spot size of~2 μm. The backscattered signal was collected through the same objective lens and dispersed by a 2400 g/mm grating with a spectral resolution of~0.4 cm −1 . The laser power was kept below 0.25 mW to avoid sample damage and laser heating.

DATA AVAILABILITY
All data needed to evaluate the conclusions in the paper are presented in the paper and the Supplementary Information. Additional data are available from the corresponding author on reasonable request.