Ultrafast broadband optical spectroscopy for quantifying subpicometric coherent atomic displacements in $WTe_2$

Here we show how time-resolved broadband optical spectroscopy can be used to quantify, with femtometer resolution, the oscillation amplitudes of coherent phonons through a displacive model without free tuning parameters, except an overall scaling factor determined by comparison between experimental data and density functional theory calculations. $WTe_2$ is used to benchmark this approach. In this semimetal, the response is anisotropic and provides the spectral fingerprints of two $A_1$ optical phonons at $\sim$8 and $\sim$80 $cm^-$$^1$. In principle, this methodology can be extended to any material in which an ultrafast excitation triggers coherent lattice modes modulating the high-energy optical properties.

Several fundamental properties of materials, including the electrical and thermal conductivities, are influenced by lattice vibrations [1]. Recently, the possibility to control such behaviors by resonant coupling of ultrashort light pulses to specific lattice modes has been proved [2] [3]. Conversely, coherent phonon spectroscopy has emerged as a powerful method to directly observe, in the time domain, coherent lattice vibrations [4][5] [6]. Despite the many studies, coherent phonons have been exploited primarily to characterize ground state properties [7][8] [9] [10], rather than controlling the material properties. This observation motivates the quest to learn how specific lattice vibrations affect the electronic behavior in the vicinity of the Fermi energy. Structural dynamics experiments, such as time-resolved x-ray diffraction [11] and time-resolved electron diffraction [12], are currently used to measure the amplitude of synchronized collective excitations of the atoms in a solid, i.e. coherent phonon modes. Yet, subpicometer displacements in complex materials are very challenging to be resolved.
In time-resolved reflectivity experiments, such modes appear as ultrafast oscillations of the probe signal [13] [14]. Their amplitude has been related to the atomic shifts for single-element materials [15] [16] by using singlefrequency measurements. In materials with complex unit cells, a larger number of phonon modes is present and a method that avoids correlations in the estimation of individual amplitudes is required. Here, we report on a novel approach to estimate the non-equilibrium atomic displacements of coherent optical phonon modes by detecting the modulation induced in a broadband reflectivity probe signal, which is central to obtain reliable signatures of the modes. Density functional theory (DFT) calculations are applied to a displacive model to simulate the time-resolved optical reflectivity signals and to estimate the coherent phonon amplitudes in the tens of femtometers regime.
At present, our method is applied to the orthorhom-bic semimetallic (Td) phase of tungsten ditelluride (WTe 2 ), a transition metal dichalcogenide that has recently gained interest for showing an unusually high and non-saturating magnetoresistance [17] along with a possible type-II Weyl semimetal character [18]. WTe 2 has a layered structure belonging to the space group Pmn2 1 [19] (see Fig. S1 in [20]) with in-plane covalent bonding and principally van der Waals interactions holding together the individual layers. Furthermore, it has been shown that its electronic, optical and topological properties are influenced by strain forces [18] [21], which can be induced through non-equilibrium perturbations [22] allowing an ultrafast control of the functionalities of this material. Expressly, we focus on two A 1 coherent optical phonon modes at ∼8 cm -1 and ∼80 cm -1 . They are comprised of non-equilibrium displacements of the atoms along the y and z crystallographic directions. The ∼8 cm -1 optical phonon is a uniform in-plane shift of the atoms, estimated to be ∼350 fm using a ∼230 µJ/cm 2 absorbed pump fluence, while the ∼80 cm -1 mode corresponds to atomic displacements of few tens of femtometers which depend on the specific atom. A more detailed description will be given further on in this work.
High-quality tungsten ditelluride samples were grown as reported in [17]. The presence of defects was previously studied [17] [23] and has negligible impact on our estimate of the average displacements. In order to identify the in-plane crystallographic axes, LEED images were acquired under ultra-high vacuum conditions. Timeresolved reflectivity experiments (sketch of the set-up in Fig. 1(a)) were performed using a Ti:sapphire femtosecond (fs) laser system, delivering, at a repetition rate of 250 kHz, ∼50 fs light pulses at a wavelength of 800 nm (1.55 eV). The broadband (0.8-2.3 eV) supercontinuum probe beam was generated using a sapphire window.
DFT simulations were carried out using normconserving (NC) [24] scalar relativistic [25] pseudopoten- tials with the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) parametrization for the exchange-correlation functional [26] chosen from the PseudoDojo database [27] [28]. Structural optimizations and zone-center optical phonon calculations in the framework of density functional perturbation theory (DFPT) were performed using the Quantum Espresso (QE) [29] suite of codes. We calculated the diagonal macroscopic dielectric tensor components through the Yambo code [30] at the independent-particle (IP) level, starting from the calculated wavefunctions and eigenvalues obtained through QE [20].
The normalized time-resolved differential reflectivity (DR/R), measured in a near-normal incidence geometry, revealed a large anisotropy consistent with the two-fold in-plane symmetry of WTe 2 ( Fig. 1(b)) also reported for equilibrium electrodynamics [31].
In Figs. 1(c),(d), the broadband probe beam polarization was set parallel to the x (η||x) and y (η||y) crystallographic axes respectively with pump beam polarization kept perpendicular to the probe one. Measurements were acquired with a ∼710 µJ/cm 2 absorbed pump fluence and 1.55 eV pump photon energy at T=295 K. A region of the spectrum was disturbed by the scattering of pump beam photons from the sample, so it was removed.
The temporal evolution of the reflectivity is directly related to the electronic and ionic degrees of freedom [4][13] [14]. In general, the DR/R is a superposition of signals due to transient variations of electronic density of the states and population. We describe the temporal evolution after the perturbation (time-zero, t=0) as where G(t) represents the pump and probe pulses cross-correlation and A i (hν), B(hν) and C j (hν) denote the amplitudes of three different phenomena: i) electronic relaxation phenomena with time constants τ Ri , ii) heating contribution with a characteristic time τ heat and iii) oscillations due to coherent phonons with angular frequency ω j , initial phase φ j and decay time τ P j .
After the pump pulse excitation, electron-electron and electron-phonon scattering processes constitute the main incoherent relaxation phenomena in the first hundreds of femtoseconds and are responsible for the initial exponential decay of the DR/R signal. After a few picoseconds, the differential reflectivity DR/R reached a plateau, showcasing an offset with respect to the equilibrium value. It vanished in almost one nanosecond [20]. Due to its timescale, this effect is likely to be associated with local lattice heating, as described in [32] for transition metals. For our purposes, we can consider the offset as a constant over the temporal range explored in Figs. 1(c),(d).
The temporal profiles are very well fitted by using Eq. (1) [20]. Their sign and magnitude vary with the photon energy, while the main relaxation presents a time constant which is almost the same and equal to ∼1 ps with only small differences, as previously found in [33] at hν=1.55 eV.
Subtracting an exponential fit function from the experimental signals allows the extraction of the coherent component of the DR/R signal. The resulting signal displays marked periodic modulations in the time domain, arising from the excitation of coherent phonon modes. The most prominent oscillating features, as detected by Fast Fourier Transform (FFT) and in [33], have frequency of 7.9±0.4 cm -1 and 79.7±0.4 cm -1 . For notation purposes, we refer to these modes with the labels 8 cm -1 and 80 cm -1 . The associated time-decay constants, derived by fitting the modulations with two exponentially-damped cosine waves, are τ P 1 =77±4 ps for the 8 cm -1 mode and τ P 2 =12.1±0.2 ps for the 80 cm -1 mode. DFPT results indicate that these frequencies are linked to A 1 zone-center optical phonon modes, represented in Fig. 2 the experimental frequency for the 8 cm -1 mode with respect to the DFPT result (10.30 cm -1 ) can be attributed to a redshift as temperature or fluence are increased [34]. The marked difference between the time constants of the two modes could be linked to the different type of perturbation of the interatomic bonds induced by the associated displacements. The 8 cm -1 mode ( Fig. 2(a)) only alters the interplanar, mainly van der Waals, interactions. Differently, the 80 cm -1 mode ( Fig. 2(b)) induces a modification of the in-plane covalent bonds. A more detailed DFPT analysis could provide a rigorous basis for this intuitive argument, which is beyond the scope of the present work.
Weaker additional A 1 contributions at 116.5±0.4 cm -1 , 132.2±0.4 cm -1 and 210.2±0.4 cm -1 frequency appear as beats in the DR/R signal in the first hundreds of femtoseconds.
When switching the probed direction from η||x to η||y, a π phase change for the 8 cm -1 mode (Fig. 2(c)) was clearly registered for most of the probe photon energies. Analogous π phase changes were found by comparing the temporal profiles taken at different probe photon energies using the same polarization for the same mode ( Fig.  2(d)). These phase differences can be explained in terms of the peculiar anisotropy of the dielectric function and were reproduced through our numerical simulations (see discussion for Fig. 4).
In order to model the effects of the coherent phonons on the reflectivity, we first calculated the wavefunctions and eigenvalues for the equilibrium configuration using QE and the macroscopic dielectric tensor diagonal components and the reflectivity curves through Yambo. Then, we considered four out-of-equilibrium configurations labeled as 0 and π phases, corresponding to displacements in opposite directions, for each of the coherent phonon modes, using the eigendisplacements obtained through DFPT with respect to the equilibrium positions. The optical phonons led to modifications of the electronic band structure that, although small, were beyond the numerical accuracy and were confirmed by linearly rescaling the effects of larger perturbations [20]. For each configuration we calculated the associated reflectivity and DR/R with respect to the equilibrium configuration. These curves describe the effect of the optical phonons on a quasi-equilibrium adiabatic system. Indeed, a few picoseconds after time-zero, the system has relaxed through electron-phonon scattering processes and the incoherent part of the DR/R signal reaches a plateau having a decay time much larger than the coherent phonon period and damping timescale.
Experimentally, the effects of each single optical phonon mode on the DR/R were isolated with the following procedure. For the 8 cm -1 phonon, we took the mean value over one period of the 80 cm -1 modulation around two consecutive 8 cm -1 extrema (Fig. 3). For the 80 cm -1 oscillation, we directly used the spectral profiles at the extrema, considering that the period of the slower phonon is about ten times larger. The spectral profiles were collected after the main relaxation peak, as precisely indicated in the insets. The difference between the spectral profiles at the maximum (0 phase) and minimum (π phase) of the coherent phonons, named with respect to the η||x DR/R, gives an unique signature of the phonon mode.
Analogously, in our calculations, we considered the difference between the calculated DR/R curves at 0 and π, multiplied by an exponential factor derived from the experimental data, to account for the damping of the optical modulation. The DFT predictions closely match the experimental results for both the 8 cm -1 and 80 cm -1 modes by rescaling the atomic displacements, whose rel- ative amplitude is given by DFPT, with a global multiplicative factor, common to both polarizations (Fig.  4). The agreement is particularly good for η||x where the phonon effects are most prominent. Peculiar features, such as their trend and change of sign, are maintained throughout the relaxation process as previously illustrated in Fig. 2 From this comparison we can quantify the magnitude of the ionic displacements right after the perturbation, when the atoms of different unit cells move uniformly in the probed region. The estimated displacements have to be regarded as average displacements over the probed volume and photon energies. The penetration depth, derived from the optical data in [31], is almost constant and similar for the two perpendicular polarizations in the probed spectral region, with minor deviations towards the infrared for y direction.
In order to evaluate very small atomic displacements, we repeated the experiments at lower fluence (∼230 FIG. 4. Comparisons between the experimental and DFT calculated DR/R difference between the 0 and π phases, showcasing the effects of the 8 cm -1 and 80 cm -1 optical phonons for (a) η||x and (b) η||y using a ∼710 µJ/cm 2 absorbed pump fluence and T=295 K.
µJ/cm 2 ). We checked that the DR/R spectral shape does not change at lower temperature or lower fluence, where we showed that the phonon effects can be considered linear with the deposited energy [20]. This is not the case for ∼710 µJ/cm 2 , since the phonon effects measured at that fluence are only twice those at ∼230 µJ/cm 2 . At ∼230 µJ/cm 2 , for the 8 cm -1 mode, the initial positions shifts are found to be ∼350 fm along the y direction. For the 80 cm -1 mode, the displacements are quite distinct for the twelve atoms in the basis cell. The four tungsten atoms remain almost fixed (displacements smaller than 10 fm) along the y direction whereas they shift by ∼40 fm along the interlayer z direction. Four tellurium atoms are displaced by ∼90 fm along the y direction and by ∼25 fm along z. The remaining four tellurium atoms are displaced by ∼15 fm along the y direction and by ∼30 fm along z. None of these modes involves movements along the x crystallographic direction.
We highlight that with the present approach the atomic displacements can be evaluated with a precision of a few femtometers without free tuning parameters, except a scaling factor determined by an overall comparison between experimental data and numerical simulations of the broadband DR/R signal. Concerning WTe 2 , a good agreement between the measured and calculated phonon effects on the optical properties for both in-plane crystal axes has been obtained. This finding confirms that the macroscopic lattice distortions excited in WTe 2 at 8 cm -1 and 80 cm -1 can be entirely mimicked by the coherent population of selected zone-center A 1 lattice modes. We stress that the possibility to measure the optical properties on a wide spectral range is fundamental to extract reliable quantitative results about the magnitude of eigendisplacements of multiple and intertwined phonon modes. Finally, we mention that the method described here is not system-specific and in principle can be extended to any crystalline material, provided that its high-energy optical properties are affected by the coherent motion of atoms. In perspective, our findings can pave the way to the design of tailored devices in which the coherent lattice motion (at selected frequencies and up to large displacements) is exploited to finely tune the functional properties of semiconducting and metallic systems.  Fig. S1, ball-and-stick projections of the crystal structure of the layered orthorhombic phase of tungsten ditelluride (WTe 2 ) are shown. The crystallographic data were taken from [1]. All the crystallographic representations were generated using the XCrySDen [2] software.   Error intervals are expressed as ±σ, one standard deviation.

-Measurement at 1.03 eV pump photon energy
In order to check that the reflectivity variation DR/R is independent from the pump photon energy, we repeated part of the measurements by pumping at 1.03 eV photon energy ( Fig. S2), obtained thanks to an optical parametric amplifier (OPA). We verified that the DR/R is qualitatively similar to the results obtained with the 1.55 eV (800 nm) excitation.
Moreover, since in this experimental condition the pump lies at the edge of the measured spectral range, the spectral region around 1.55 eV is not covered by scattered pump photons.
This allowed us to confirm that no structured features are present in this range, that was hidden in the datasets of the main paper acquired with 1.55 eV pump photon energy (Figs. 1(c),(d), main paper).

-Time-domain analysis
In Fig. S3 we report selected profiles extracted at two photon energies from the datasets reported in Figs. 1(c),(d) in the main paper, alongside with the Fast Fourier Transform (1) in the main paper) from the temporal profiles at 1.40 eV acquired at ∼710 µJ/cm 2 absorbed pump fluence, T=80 K and T=295 K and for η||y; the beats are due to phonon modes with frequency higher than 80 cm -1 .
features are at ∼8 and 80 cm -1 . These modes give the largest contribution to the coherent component of the DR/R signal as shown in Fig. S4, where two damped cosine waves fit the oscillatory signal. Minor contributions at higher frequencies (∼117, 132 and 210 cm -1 ) are detected and appear as beats in the oscillatory signal. Their contribution becomes more relevant as the temperature is reduced (Fig. S5).
In the visible spectral region (1.70-2.10 eV), the DR/R cannot be fitted by a single exponential decay (Eq. (1) in the main paper) for probe polarization η||x. In presence of contributions with opposite sign and comparable time constants, the DR/R reaches its maximum value at longer delay, analogously to the case in Fig. S3(a). This delay can be as large as 1 ps, alongside with a slower decay of the signal. Investigating the dependence on fluence, we learned that (Fig. S6) the delay increases with the deposited pump energy. At variance, in the infrared the peak position is independent on the fluence. For this reason, for example, it is not possible to effectively describe the behavior of this material with just a single Drude-Lorentz term [3] in the dielectric tensor for the measured spectral range.
For η||y, we do not observe such phenomenon. This may be connected to different matrix elements for the involved transitions in the two distinct linear probe polarizations.

-Low fluence and low temperature measurements
The DR/R broadband image measured at ∼230 µJ/cm 2 absorbed pump fluence and T=295 K (Fig. S7) presents the same characteristics of that measured at higher fluence (see Figs. 1(c), 3(a) and 4(a) in the main paper). Such features are also well captured by ab-initio results obtained with properly rescaled eigendisplacements.
At T=80 K, phonons exhibit a blueshift in their frequencies with respect to T=295 K, as previously reported in [4] [5]. Compared to T=295 K the contribution from coherent phonons with frequency at higher than 80 cm -1 becomes more prominent in the DR/R images measured at ∼710 µJ/cm 2 (Figs. S8, S9) and ∼230 µJ/cm 2 (Fig. S10) absorbed pump fluence. The phonon effects for η||x shown in Fig. S11 are analogous to the room temperature results reported in Figs. 4(a), S7(d), although it becomes harder to separate them using the same method. This becomes especially challenging for η||y, where the amplitude of phonons with frequency higher than 80 cm -1 is comparable to the lower modes (Fig. S5) throughout the investigated spectral region. In this case, a more advanced procedure involving a timeweighted average over multiple 0 or π phase profiles to determine the initial amplitudes of the atomic shifts could be employed.

PART 3 -COMPUTATIONAL RESULTS
Density functional theory (DFT) simulations were carried out using norm-conserving (NC) [7] scalar relativistic [8] pseudopotentials with the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) parametrization for the exchange-correlation functional [9] chosen from the PseudoDojo database [10] [11]. An orthorhombic simulation cell (Figs. S1(a)-(c)), with lattice constants a=3.477Å, b=6.249Å and c=14.018Å along the x, y and z directions respectively taken from the crystallographic data in [1], was used. It contains four tungsten and eight tellurium atoms, with 14 and 16 valence electrons respectively, for a total number of 184 electrons per unit cell.
Structural optimizations and phonon calculations in the framework of density functional perturbation theory (DFPT) were performed using the Quantum Espresso (QE) [12] suite of codes with a plane wave kinetic energy cutoff of 70 Rydberg and a 12×10×6 uniform kpoint mesh for integrations over the Brillouin zone. Van der Waals contributions were not included, since, although in general relevant to describe layered compounds, in this specific case they do not improve the description of structural and electronic/optical properties of the system. Td-WTe 2 presents 33 zone-center optical phonon modes with four possible symmetry representations given by the C 2v point group [1]. In the following sections, we will focus on the two lowest energy A 1 modes. We calculated the diagonal macroscopic dielectric tensor components through the Yambo code [13] starting from the wavefunctions and eigenvalues obtained through the PW package of QE using the same kinetic energy cut-off and a 16×14×10 k-grid. Convergence tests showed that these parameters ensure the numerical accuracy needed to appreciate the relative effects of subpicometer atomic displacements on the electronic states in the near-infrared and visible spectral regions. The number of conduction bands to be considered in the calculations was also determined from convergence tests. The reflectance curves were derived at the independent-particle (IP) level using a linear response approach.

-Equilibrium reflectance
We report the calculated reflectance curves along the two perpendicular polarization directions in Fig. S14, compared with the reflectance data found in [14]. We included 119 bands for the response function to reproduce the experimental data with a very good quantitative agreement. Although the experimental data was measured at 10 K and DFT results were obtained at T=0 K, they can be still considered adequate for our DR/R comparisons at T=295 K in good approximation, since the optical properties in this spectral region are only slightly affected by temperature [15].
FIG. S14. Equilibrium reflectance curves in the two perpendicular probe polarizations; the experimental curves were derived from the data present in [14] measured at 10 K.

-Coherent optical phonon effects on the energy levels
The subpicometer atomic displacements in the out-of-equilibrium structural configurations lead to modifications of the electronic band structure (Fig. S15) of the order of a few meV.
Such variations are small, but beyond our relative numerical accuracy. Furthermore, we checked that these effects are almost linear with the displacements if we increase them by a factor up to twenty. We focus on the extrema of the electron and hole bands along Γ-X. In  Fig. S17), while for the 80 cm -1 modes the difference becomes ∼1.7 (0) meV. By properly tuning the Fermi level, these variations, although small, could have an impact on the extremely high magnetoresistance in WTe 2 [16], which is ascribed to a delicate balance between electron and hole carriers. Additional contributions may come from carrier mobility and the magnetic field [17] [18]. We caution the reader that, as stated in the main paper, we used scalar relativistic pseudopotentials. The inclusion of spin-orbit coupling (SOC) using fully relativistic pseudopotentials would provide a more accurate description of these effects, since it is known that SOC plays an important role on the electronic bands around the Fermi level in WTe 2 [19]. for the 80 cm -1 phonon mode; the 0 and π phases are labeled as in the main paper, taking as reference the η||x polarization; the displacements used for the graphs are double with respect to the low fluence values reported in the main paper to highlight the phonon effects.

-Bands contribution to the response function
To accurately describe the effects induced by the optical phonons, it is necessary to include in the response function an appropriate number of bands around the Fermi level.
The "converged" curves shown in Fig. S18 were obtained including 119 bands, which is safe enough compared with the 184 electrons per unit cell and the spectral range under consideration for the transitions. In the same figure, we show the effect of varying the number of "conduction bands", i.e. the totally unoccupied bands above the Fermi level. We notice that the main features of the DR/R profile are quickly captured even using few conduction bands, while to correctly compare the DFT amplitudes with the experimental data and give a reasonable estimate of the ionic displacements, more than 15 conduction bands have to be included.
However, since major features are still retained adding only the first conduction band in the response function, we deduce that transitions having final states that belong to bands crossing the Fermi level are involved. This supports the claim that for this material the phonon effects cannot be attributed to single transition in k-space, but to multiple contributions.