Coherent energy exchange between carriers and phonons in Peierls-distorted bismuth unveiled by broadband XUV pulses

In Peierls-distorted materials, photoexcitation leads to a strongly coupled transient response between structural and electronic degrees of freedom, always measured independently of each other. Here we use transient reflectivity in the extreme ultraviolet to quantify both responses in photoexcited bismuth in a single measurement. With the help of first-principles calculations based on density-functional theory (DFT) and time-dependent DFT, the real-space atomic motion and the temperature of both electrons and holes as a function of time are captured simultaneously, retrieving an anticorrelation between the $A_{1g}$ phonon dynamics and carrier temperature. The results reveal a coherent, bi-directional energy exchange between carriers and phonons, which is a dynamical counterpart of the static Peierls-Jones distortion, providing first-time validation of previous theoretical predictions.


I. INTRODUCTION
Ultrashort optical pulses can excite materials away from thermal equilibrium, where interactions between the degrees of freedom of solids can be explored. Electronlattice interactions in particular determine numerous material functionalities: they set fundamental limits on charge-carrier mobilities in photovoltaic applications [1], lead to the formation of electron pairs responsible for superconductivity [2], and drive the appearance of chargedensity waves [3]. The semimetal bismuth is a particularly important system for which electron-phonon interactions are being explored. A key property of bismuth is its equilibrium crystal structure: the two atoms of its primitive unit cell of the A7 structure [4] are displaced towards each other compared to the simple cubic structure. This creates a dimerized structure along the [111] direction ( Fig. 1a) in which the Bi-Bi bond length is alternatively short and long -akin to the one-dimensional distortion described in the seminal work of Peierls [5]. While the structural distortion by itself has a detrimental energy cost, it allows to lift the degeneracy of valence electronic states. This results in electrons occupying states of lower energy, which more than counterbalances the cost of the structural distortion. This phenomenon, known as a Peierls-Jones distortion [6,7], is analoguous to a Jahn-Teller effect for extended systems.
Photoexcitation gives the opportunity to tip the delicate energy balance of the system: when electrons are optically promoted from the valence states, less energy is now gained by the lifting of degeneracy. This makes the structural distortion less favorable, with the quasiequilibrium atomic position shifting towards a more symmetric and metallic structure. Time-resolved techniques allow to watch the system evolve towards this new equilibrium, uncovering spectacularly large and fast atomic motion [8]. During this relaxation, there is a bidirectional energy exchange between electronic and structural degrees of freedom, which, if it can be captured, grants a unique view into the electron-phonon interactions at the heart of a Peierls-Jones distortion. The light-driven electronic properties have been thoroughly studied using optical pump-probe techniques [9][10][11] as well as angle-resolved photoemission [12], while the transient atomic motion has been directly measured using ultrafast x-ray diffraction [13,14]. However, probing the electron, hole, and phonon responses simultaneously in a single experiment is out of reach for these approaches. Since the electronic and structural dynamics are inherently strongly coupled, such a capability is highly desirable to obtain a complete and unified picture of the photoinduced dynamics.
In the investigation presented here, photoexcited bismuth is probed using broadband transient reflectivity in the extreme ultraviolet (XUV) spectral range at the Bi O 4,5 edge. With the help of first-principles simulations based on density-functional theory (DFT) [16,17] and time-dependent DFT (TDDFT) [18,19], we measure both lattice motion and carrier dynamics, garnering separate information on electron and hole temperatures. These carrier temperatures are found to have distinct dynamics and, surprisingly, to oscillate in anti-phase with the oscillations of atomic positions -an hitherto unobserved behavior. We explain this phenomenon by en- tropy conservation arguments that can be seen as a dynamic equivalent of the static Peierls-Jones distortion. Our results highlight the coherent bi-directional energy exchange between electron and phonon sub-systems on the femtosecond timescale, and provide a compelling illustration of the use of time-dependent measurements to gain insight into time-independent interactions.

II. OVERVIEW OF REFLECTIVITY MEASUREMENTS, XUV TRANSITIONS, AND OBSERVATIONS
The overview of the basic experiment is outlined in Fig. 1. The partial density of states (PDOS) of Bi is computed using DFT (see Appendix C 1), shown in Fig. 1c. About 24 eV below the Fermi level (E F ) lie two groups of bands with 5d character, which are split by about 3 eV due to strong spin-orbit (SOC) coupling. The topmost valence bands and lowest conduction bands mainly have 6p character. XUV radiation of appropriate energy can photoexcite the 5d semicore electrons to the Fermi level and the lowest empty 6p bands, a process that corresponds to the transition at the Bi O 4,5 edge.
In the experiment, a single-crystal of bismuth with (0001) orientation is probed by XUV radiation created by high-harmonic generation. Using sub-5 fs near-infrared (NIR) driving laser pulses, a xenon high-harmonic target gas and subsequent spectral selection (see Appendix A 1), a continuous spectrum spanning 20 − 36 eV is obtained (Fig. 1b). The static reflectivity of the sample, probed at 66 • from its surface normal (Fig. 1d), shows a maximum at 0.6 eV above the Fermi level (which lies at E F = 24 eV [15]) corresponding to the maximum in the conduction-band density of states (DOS) as obtained from DFT (Fig. 1c) and in agreement with earlier work [15]. The next large peak 3.10 ± 0.01 eV higher in energy corresponds to transitions from the deeper spinorbit split 5d bands to the empty states. A detailed assignment of features is presented together with the band structure in Appendix C 2.
The sample is photoexcited in pump-probe measurements by sub-5 fs NIR pump pulses at an intensity of 0.5-1.5 mJ/cm 2 , using the same setup as in Refs. [20,21] The pump-induced changes in XUV reflectivity, measured by scanning the pump-probe delay, processed using an edgereferencing method to reduce correlated source noise (see Appendix A 2 and Ref. [22]), are shown in Fig. 1e. The transient reflectivity shows picosecond-long signals both above and below the Fermi level, among which the most notable features are oscillations of the XUV reflectivity in several energy ranges, but predominantly close to E F . These oscillations are chirped, with instantaneous frequencies increasing from 2.4 THz to 2.8 THz as the timedelay gets longer (see Appendix B for a time-frequency analysis). This is a clear signature of the A 1g coherent phonon mode (2.93 THz at equilibrium [23]) associated with the Peierls-Jones distortion, whose softening under carrier excitation is well known [14,24]. While oscillations of XUV reflectivity have already been observed and attributed to coherent phonon motion [25][26][27], a detailed understanding of their origin in terms of interactions between excited carriers (electrons and holes) and the lattice (phonons) is lacking. The oscillations could stem from energy shifts of valence bands (due to the coupling of valence electrons to coherent phonons via electron-phonon interactions [12]) but also from corelevel shifts, which can be particularly sensitive to bond lengths [28][29][30], or even from modulations of either hole or electron temperatures around the Fermi level [31]. To elucidate their origin, we perform ab initio calculations based on TDDFT with the aim of understanding the time-dependent (dynamical) effects of the coherent A 1g displacement on the XUV spectrum of the material.

III. FIRST-PRINCIPLES CALCULATIONS OF XUV ABSORPTION
The bismuth O 4,5 absorption spectrum is first calculated at equilibrium using TDDFT (see Appendix C 1). At the incidence angle used in the experiment, the XUV penetration depth is 13.6 nm at 24 eV (the Fermi level energy). This is much larger than the extent of the surface states of Bi which have two-dimensional character [32]. Thus, in what follows we focus solely on the bulk electronic and phonon properties of bismuth.
The theoretical equilibrium distance between nearest Bi atoms along the trigonal axis is 5.50Å, which is in very good agreement with the experimental value at room temperature of 5.54Å [14,33]. To allow for comparison, the experimental absorption spectrum is obtained by inverting the measured reflectivity of Fig. 1d using a Kramers-Kronig procedure previously described [21]. Fig. 2a displays the experimental and theoretical absorption spectra, with the latter being rigidly shifted down by 2.7 eV to match the experimental position of the absorption edge. Such a shift is necessary because of the overdelocalization of the d electronic states when using (semi-)local exchange-correlation (xc) functionals in (TD)DFT, which are responsible for self-interaction errors [34] and lead to the imprecise position and shape of such states. Apart from this shift, the overall shape of the computed spectrum is in good agreement with the experimental one, with discrepancies only in the relative intensity of the above-edge peaks, which correspond to the various conduction band extrema [15]. It is worth noting that more accurate absorption spectra can be obtained by solving a Bethe-Salpeter equation on top of a GW calculation [35][36][37][38]; however, for the purposes of this work the absorption spectrum computed at the TDDFT level is fully sufficient.
The calculation is now repeated while changing the Bi-Bi distance in the unit cell along the trigonal axis, which is denoted by δ and is expressed in percent of the initial Bi-Bi distance; positive δ corresponds to a repulsion of Bi atoms away from each other, as illustrated in Fig. 2b. This approach allows us to observe changes in the imaginary part of the refractive index of Bi due to the coupling of photoexcited carriers with the lattice vibrations. The computed absorption spectrum is shown in Fig. 2a for δ evolving from -5% to 5%. The computation shows that the displacement along the phonon coordinate leads to a blueshift (resp. redshift) of the edge for a contraction (resp. elongation). For small displacements, the shift is linear with δ at a rate of 9.7 meV/pm. This value can be seen as the local gradient of the potential energy surface of the semicore-excited state along the A 1g phonon coordinate. Several eVs above the absorption edge, the behavior becomes more complicated. This is because XUV absorption probes the entire Brillouin zone and thus the contributions of higher conduction bands at various momenta are summed together. Since the electron-phonon coupling varies strongly with momentum [12], this results in a more complicated trend than a rigid shift. Additionally, absorption from the 5d 3/2 core-level further complicates the spectrum above 27 eV (see Appendix C 2). We thus focus on the shift near the absorption edge. We can elucidate the origin of this shift by looking at the DOS, shown for a repulsion of δ = +1% in Figs. 2c-d. As is customary in DFT, the computed DOS at and out of equilibrium must be arbitrarily aligned to some energy; we choose here to align them at the Fermi energy. This is contrary to previous calculations interpreting photoemission experiments in bismuth where the electronic bands were aligned at the 5d level energy [12], as the removal of an electron from the system during photoemission does not keep the Fermi energy constant. The conduction band (CB) states shift towards the Fermi level by about 100 meV while the binding energy of the 5d semicore states increases by approximately 45 meV. This is consistent with the fact that Bi in the excited state is more symmetric and metallic [39,40]: the indirect negative gap is closing, while the Coulombic repulsion between ionic cores is reduced, increasing the binding energy of the 5d levels. Here, the CB shifts more than the 5d levels, resulting in a decrease of the transition energy; however, it is important to note that this will not necessarily be true for all systems. Finally, we briefly explored the impact of the A 1g phonon mode on the electronic structure near the L point, where doping-induced topological phase transitions were observed [41]. As shown in Appendix C 3, we do not find evidence of phonon-induced phase transition, but it must be stressed that our level of theory is not well suited to describing such fine features at low energies.

IV. DECOMPOSING THE TRANSIENT REFLECTIVITY TRACE: CARRIER AND PHONON DYNAMICS
Now that we linked the A 1g phonon motion with a shift of the absorption edge, it is possible to extract quantitative information from the XUV transient reflectivity measurement. The procedure follows established methods successfully applied to elementary semiconductors [20], transition metal dichalcogenides [42], and perovskites [43]. We decompose the measured dynamics into three contributions: phonon motion (edge shift), holes photoexcited in the valence band and electrons in the conduction band. The decomposition attempts to recover three variables at each time-delay: the electronic and hole temperatures, T elec (t) and T holes (t), as well as the phonon- distinct behaviors and spectral signatures, which are separable due to the broadband nature of the spectroscopic observable. We remark that due to the reflection geometry, the total transient reflectivity is not given by the simple sum of each component, and that the electrons and holes contribute signals away from their respective bands [20]. While the ability of the model to extract accurate absolute carrier temperatures deserves further exploration, a wealth of information on the separate dynamics of holes, electrons, and phonons can be extracted. Figure 4 shows the extracted variables (edge shift and two carrier temperatures) for two sets of measurements taken at different pump fluences. The phonon-induced shift (Fig. 4a) oscillates with a very visible phononsoftening and fluence-dependent motion amplitude. The shift can be converted to real space motion using the gradient determined by TDDFT of 9.7 meV/pm: we obtain maximum excursions of 1.9 ± 0.1 pm and 2.6 ± 0.1 pm for 0.9 mJ/cm 2 and 1.4 mJ/cm 2 fluences, respectively. This conforms well with the values measured using x-ray diffraction [14] of 1.7 pm and 2.9 pm at 1.2 mJ/cm 2 and 1.7 mJ/cm 2 . Furthermore, in Appendix E we directly juxtapose the extracted time-dependent atomic positions with previous ab initio calculations [31], for a variety of pump fluences. The comparison (Fig. 12) indicates a striking agreement in both the magnitude of actual atomic displacement and the frequency of the motion, for both fluences used in the experiment. The results and comparisons strongly supports that the all-optical XUV technique yields accurate and quantitative mea-surements.

V. DISCUSSION
In addition to the atomic motion, the measurements contain another piece of the puzzle -the dynamics of both electrons and holes. We stress that the accuracy of the method in recovering meaningful absolute values for electron and hole temperatures deserves further investigation before being interpreted, as was recently done for nickel [44]. However, the relative evolution of the carrier temperatures is insightful; as seen in Fig. 4b, they decay with markedly different time constants of τ elec = 1.5±0.4 ps and τ holes = 3.8 ± 1.2 ps. This suggests that each carrier couples differently to the coherent phonon motion, which might be explained by their widely different optical masses [45]. Both temperatures show a delayed maximum at ∼ 0.3 ps which is consistent with the carrier thermalization time [45]. Most intriguingly, the experiment captures pronounced oscillations of both temperatures, at frequencies closely following the atomic motion. However, these oscillations are out of phase: when the atoms are driven away from their equilibrium position (edge shift is maximum, largest elongation), the carrier temperature diminishes, and vice-versa. The decay time of the oscillations is slightly longer for lower pump fluences, following the coherent phonon behavior [31].
The excitation of coherent phonons is well described by the displacive excitation mechanism [46], which postulates that the quasiequilibrium A 1g nuclear coordinate depends linearly on electronic temperature [47]. This leads to damped oscillations of atomic position, as observed here and in many other experiments. Importantly, an angle-resolved photoemission study [48] managed to measure both photoelectron diffraction and emission signals in the same experiment, gathering simultaneous electronic and structural information. This way, the relative phase between phonon displacement and electronic shifts was found to match the DECP process [23]. However, the observation here is fundamentally different: we probe the relative phase between phonon displacement and electronic temperature. In fact, the DECP mechanism does not explain why the electronic temperature itself oscillates. Instead, we relate to a reported thermodynamical model that anticipated this result [31], which our experiment confirms for the first time. This model succeeded in reproducing x-ray diffraction measurements of bismuth [14] in a particularly remarkable way, since it only has two adjustable parameters. In addition, it predicted antiphase oscillations of the electronic temperature and rationalized them as the way by which the system counteracts the increased entropy created by atomic displacements.
Consider the static Peierls-Jones distortion, illustrated in Fig. 5a: at equilibrium, the total energy of the system (sum of an ionic and an electronic part, respectively repulsive and attractive) is reduced by going to a less metallic and less symmetric state. The dynamic behavior of the system seems to be similar, with phonons and temperatures oscillating in quadrature as the system attempts to balance structural and electronic energies. Understanding the basis of the thermodynamic model of Giret et al. [31], we propose the mechanism sketched in Fig. 5b: first, the laser pulse suddenly deposits energy into the system, which increases the carrier temperatures and suddenly changes the quasi-equilibrium position u eq along the phonon coordinate. The total energy thereby deposited is slowly relaxed by energy transfer to the lattice on time scales longer than 10 ps [14,31]. Thus, on the femtosecond timescale studied here, the total energy, the entropy and u eq are all almost constant after laser excitation. The new value of u eq set by photoexcitation then drives oscillations of atoms around the new position. When the atomic displacement u is maximum, the system is temporarily more symmetric, and in turn, more metallic (see band structures in Fig. 5b), which increases the electronic entropy [31]. Because the entropy must remain almost constant, this excess has to be balanced, and the experiment shows that this is realized by a diminution of the carrier temperature. The observation of this bi-directional energy exchange is possible only because both carrier and structural degrees of freedom are probed simultaneously. In turn, the dynamic behavior of the system can be seen as a time-dependent counterpart of the static Peierls-Jones distortion.
In summary, we investigated the photoexcited response of the semimetal bismuth using transient broadband XUV spectroscopy and ab initio calculations. In a single experiment, real-space motion and carrier dynamics are quantified. We find that the atomic displacement and carrier temperatures oscillate in quadrature, which is a direct view into the energy balance of this Peierls-distorted system and confirms previous theoretical predictions [31]. Our approach can be applied to similar systems, with the limitation that materials showing valence or core-excitonic effects would require more advanced computations based on many-body perturbation theory. Examples of materials hosting Peierls distortions are the other group V elements (As, Sb, black phosphorus), known to also show dimerized structures [49]. Likewise, the structures of Te, Se, and many compounds they form with group V elements can be seen as n-merizations (n > 2) and might show comparable electron-phonon dynamics [50,51]. For instance, a similar behavior was recently observed in the 2D charge density wave material TaSe 2 [52] using angle-resolved photoemission. In comparison, the XUV transient technique is bulk-sensitive, has 5-fs time resolution and can probe all bands with allowed orbital character. As such, this work offers a complementary perspective and is particularly well suited to future studies of the highly coupled electron-phonon dynamics accompanying ultrafast lightinduced phase transitions [53]. More generally, the results provide an illustration of how time-resolved studies capturing the dynamics of several degrees of freedom are able to yield insight into fundamental, time-independent interactions between them.
Computational and experimental data that support the findings of this study are available on the Materials Cloud Archive [54].
The experimental apparatus uses a carrier-envelope phase stable titanium-sapphire laser system (Rainbow CEP4 oscillator + Femtopower amplifier from Femtolasers, pumped by a DM30 laser from Photonics Industries), which delivers 1.8 mJ, 25 fs pulses at 1kHz. They are broadened in a neon-filled hollow-core fiber (1.5 bar static pressure). A set of 12 double-angle chirped mirrors (PC70, Ultrafast Innovations) and a 2 mm-thick ADP crystal control the spectral phase of the pulse, yielding a duration of about 4 fs. Further details can be found in Refs. [20,21].
A 60:40 beamsplitter separates the incoming pulse. The most intense part (480 µJ) is focused by a 45 cm spherical mirror into a xenon-filled gas cell (∼ 20 Torr) to produce broadband high-harmonic radiation. A 50 nm titanium filter removes the generating infrared light and cuts XUV radiation above 40 eV. This is important for the spectrometer grating to function properly: it prevents second order diffraction of high energy components to overlap with the Bi edge region where the experiment is performed. The XUV is then focused onto the sample by a toroidal mirror. The remaining 40% of the initial short infrared pulses are recombined collinearly with the XUV light using a hole mirror and used as a pump. The pump and probe beams (p-and s-polarized, respectively) impinge on the sample surface with a 66 • angle from the sample normal. The pump light is then removed by an 150 nm-thick aluminum filter, and the spectrum of the reflected probe beam is measured using an aberrationcorrected diffraction grating (Shimadzu 30-006) and an XUV CCD camera (Princeton Instruments). The spectral resolution at 25 eV is σ = 13 meV as determined using atomic absorption lines.
The sample is a single-crystal of bismuth with (0001) orientation, polished with roughness <10 nm (purchased from Princeton Scientic).

Edge-referencing of transient reflectivity data
At each time delay between pump and probe, we measure two spectra I on and I off , with the pump on and off, respectively. The transient reflectivity is then defined as Ion−I off I off . The dominant noise source arises due to fluctuations of the XUV flux. They are mitigated using an edge-referencing procedure that is described in detail in Ref [22]. The edge-pixels were taken as regions without transient signal; namely, photon energies below 21 eV and between 30 and 35 eV. We verified that a slight change of these regions did not influence the final transient reflectivity. Figure 6  Under intense photoexcitation, bismuth is known to undergo bond softening [55], which makes the oscillation frequency depend on the motion amplitude. This was shown to be due to carrier-induced bond softening, rather than a consequence of lattice anharmonicity [33]. Appendix C: Computational methods
Bismuth crystallizes in the A7 rhombohedral structure, and the resulting lattice has the trigonal symmetry with two atoms in the primitive unit cell [45]. The rhombohedral cell is described by two parameters: the lattice parameter a, and the rhombohedral angle α. The position of two Bi atoms in the unit cell is described by one parameter, the so-called internal parameter u which describes the atomic positions in the crystallographic axes: one Bi atom has coordinates (u, u, u) and the other one (−u, −u, −u). We optimized these parameters using DFT with the setup described below and obtained a = 4.70Å, α = 58.00 • , and u = 0.23559, which is in very good agreement with the values obtained from xray diffraction measurements [58].
The ground-state properties of Bi (equilibrium structure and projected density of states) were computed using DFT with the following computational setup. We have used a fully-relativistic norm-conserving PP, by including 5d, 6s, and 6p electrons in the valence region, while keeping the remaining electrons frozen in the core region. The details about the PP can be found in Ref. [45]. The spin-orbit coupling effect was included in the calculations [45], and for the xc energy we have used (c) Absorption spectrum of bismuth as computed using TDDFT (shaded area), and decomposed into 5d 3/2 and 5d 5/2 (red and purple dashed line, respectively). The spectrum was rigidly shifted to higher energies by 2.7 eV.
local density approximation (LDA) with the Perdew-Zunger parametrization of the Ceperley-Alder functional [59]. The Kohn-Sham wavefunctions were expanded in PWs with a kinetic-energy cutoff of 70 Ry, while the electronic charge-density and potentials were expanded in PWs with a cutoff of 280 Ry. To sample the Brillouin zone we have used a uniform 26 × 26 × 26 k point mesh centered at the Γ point, and we have used the Marzari-Vanderbilt smearing method [60] with a broadening parameter of 5 × 10 −3 Ry. The Fermi level was determined using a more refined computational setup as explained in Ref. [45], namely by using a 50 × 50 × 50 k point mesh and the tetrahedron method [61]. This latter refinement is needed in order to describe accurately the electron and hole pockets in the vicinity of the Fermi level. The excited-state properties of Bi (absorption spectrum) were computed using TDDFT [18,19], with adiabatic local density approximation. The TDDFT equations were solved in the linear-response regime in the frequency domain using the recursive Liouville-Lanczos approach [56]. Local-field effects were taken into account. The absorption spectra were obtained by computing the macroscopic dielectric function ε(ω) which then gives access to imaginary part of the refractive index, i.e. Im[ñ](E), where E is the excitation energy. Spectra have been convoluted with a Lorentzian function with a broadening parameter of 2 × 10 −2 Ry. We have computed absorption spectra for various distortions of the interatomic Bi-Bi distances to model the changes in the refractive indexñ due to the coupling of carriers with the lattice vibrations after the perturbation by the external electric field (laser pulse). The changes in the Bi-Bi distance are expressed as δ = 100% × (d new − d eq )/d eq , where d eq is the equilibrium distance between Bi atoms, and d new is the new (increased or decreased) distance.

Origin of peaks in the absorption spectrum
Here we discuss briefly the origin of peaks in the absorption spectrum of bismuth based on the band structure and density of states (DOS). Figure 8a shows the band structure of bismuth computed using DFT. The overall band structure is in very good agreement with previous studies [45,62]; in particular, the electron pockets at the L point and hole pockets at the T point are well reproduced. Due to the spin-orbit coupling, the 5d shell is split on two subshells, O 5 and O 4 . The computed 5d bands were rigidly shifted to lower energies by 2.7 eV in order to better match (on average) the experimental binding energies of the O 5 and O 4 subshells relative to the Fermi level, which are at −24.4 ± 0.6 and −26.5 ± 0.5 eV, respectively [63]. Figure 8b shows the total DOS and projected DOS (PDOS) for the 6p bands in the energy range from 0 to 7 eV above the Fermi level. The peaks in the DOS and PDOS correspond to the extrema in the conduction bands manifold.
In the photoexcitation process, electrons are promoted from the O 5 and O 4 subshells to the lowest conduction bands due to the absorption of photons of different energy. The resulting absorption spectrum is shown in Fig. 8c, which was computed using TDDFT and split into 5d 3/2 and 5d 5/2 assuming a spin-orbit splitting of 3.1 eV and respective degeneracies of 4 and 6. It is instructive to discuss the origin of peaks in the absorption spectrum. Peaks labeled from 1 to 3 in Fig. 8c are exclusively due to transitions from the O 5 subshell to the lowest conductions bands, while peaks 4 and 5 have contributions both from the excitations from the O 5 and O 4 subshells [15].

Evaluation of the possibility of a phonon-induced topological phase transition in Bi
The energy separation between the two highest occupied electronic states at the L point in the Brillouin zone (which, for the sake of simplicity, we will refer to as the "band gap at L") is in the range from 11 to 15 meV according to experiments [64,65]. Such a small gap could be closed and even inverted, for instance by doping [41], thereby triggering a topological phase transition. Here we briefly explore whether a phonon-induced shift of the valence band could close the band gap at L.
From the computational side, reproducing such a small band gap very accurately is a significant challenge: DFT with local or semi-local xc functionals [such as LDA and generalized-gradient approximation (GGA)] largely overestimates the value of the band gap at L (see e.g. Table I in [66]). More advanced computational approaches, such as quasiparticle self-consistent GW [66], are needed in order to have accurate estimates of the band gap at L. Our current study is based on DFT-LDA (which is accurate enough to describe absorption properties of Bi in the energy range of tens of eV) and the band gap at L at equilibrium turns out to be 112 meV, which is larger than the experimental band gap at L by an order of magnitude, in accordance with [66]. Hence, DFT-LDA is not an appropriate level of theory for the investigation of the intricate details of electronic states at the L point (and hence topological properties of Bi).
Nonetheless, we computed the band gap at L as the distance between the two Bi atoms in the unit cell is changed, which corresponds to the A 1g phonon mode excitation. In addition, we repeated this study using GGA with the Perdew-Burke-Ernzerhof (PBE) parametrization [67], in order to compare the obtained trends; the results are shown in Fig. 9. We can see that at the equi- librium Bi-Bi distance, GGA gives a band gap at L of 27 meV which is greatly improved with respect to LDA, but still it is a factor of ∼ 2 larger than the experimental value. Interestingly, we see that increasing or decreasing the distance between Bi atoms gives different trends in LDA and GGA, but most importantly the band gap at L is not closed (for the magnitudes of Bi displacements considered here). This means that according to this level of theory, we cannot observe the gap closing and a subsequent re-opening, which would be a signature of the transition to a topologically insulating state [41]. However, since the band gap at L is largely overestimated at LDA and GGA, these conclusions must be tempered. More advanced studies using DFT with xc functionals of higher accuracy (e.g. hybrid functionals [68,69] or meta-GGA SCAN [70]) or the quasiparticle self-consistent GW method are needed in order to provide definitive conclusions on the effect of the A 1g phonon on the band gap at L via the electron-phonon interaction.
In order to describe carrier contributions, we begin by treating k(E, t) because it allows to directly write spacefilling effects [71]. Upon excitation, the temperature of both electrons and holes changes, which modifies the occupation of valence and conduction band states according to the Fermi-Dirac (FD) distribution, f FD . Figures 10e and h display the FD function for 300 and 2000 K, on top of the density of states (DOS) for holes and electrons, respectively. The actually occupied states are obtained by multiplying the FD distribution with the DOS of valence and conduction bands (N VB and N CB , respectively). The difference between the two thereby gives the states that become empty (resp. occupied) in the valence (resp. conduction) band. These give increased and reduced absorption, respectively (Figs. 10f and 10i). These quantities must then be broadened to account for both experimental and core-hole lifetime broadening (σ and Γ, respectively), as shown by solid lines in Figs. 10f and 10i. The changes in the imaginary part of the refractive index for electrons and holes thus read: where V is the Voigt function, a elec and a holes are timeindependent scaling factors which are needed because the amount of change in absorption as a function of the change of occupation at a given energy is unknownthis quantity depends on the oscillator strength of all semicore-to-conduction bands transitions across the Brillouin zone. We then use the Kramers-Kronig transformation with these quantities to obtain ∆n elec (E, t) and ∆n holes (E, t) (Figs. 10g and 10j). Since ∆k elec and ∆k holes are non-zero only in a small energy range, common issues with performing the Kramiers-Kronig transformation (which requires evaluating integrals from 0 to ∞ with respect to E) are circumvented.
Finally, we sum the static value of k with ∆k elec , ∆k holes and ∆k phonon , and likewise for n. The obtained complex refractive index, n(E, t) + ik(E, t), is then used in the Fresnel equation to calculate the model reflectivity. Figure 11 displays the result for this example case of a 300 meV shift and a carrier temperature of 2000 K for electrons and holes. Removing the carrier contribution highlights in which spectral region do they yield a reflectivity change. With this description in hand, the experimental transient reflectivity can be fitted by a least-square procedure. The parameters at each time delay are ∆E(t), T holes (t), and T elec (t), plus the two time-independent factors a elec and a holes , resulting in a total of N t + 2 independent variables, where N t is the number of time points. We perform the least-square minimization on the 2D matrix (time × energy) using the parallel implementation of the lsqnonlin routine in MATLAB 2018a. We used many randomly selected starting points as well as checked that lower and upper bounds did not influence the result. Finally, once the optimum was found, we computed the Jacobian at the solution point which was then used to obtain 95% confidence intervals, shown together with the optimal solution in Fig. 4. Lastly, we note that in principle there is a "cross-term" between phonon and carrier contributions: as the lattice displacement shifts the band positions, the DOS used in ∆k elec,holes should be shifted as well. However, taking into account this effect adds two unknowns at each time delay (the VB and CB shifts), making the least-squares problem impossible to solve without large uncertainties for all parameters. Fortunately, this cross-term turns out to be a second-order effect: we verified that shifting the DOS in the hole and electron contribution by the maximal edge shift that we obtained (30 meV, see Fig. 4) resulted in an error below 0.3% in reflectivity. That is an order of magnitude below all contributions shown in Fig. 3. Hence, we neglect this second-order effect to retain physically insightful quantities. The energy shift of the DOS arising due to this cross-term is most likely negligible in our case because of rather high carrier temperatures, and correspondingly wide Fermi-Dirac distributions, compounded with experimental and core-hole lifetime broadening. This might not be the case for narrower distributions. Here we directly compare the atomic motion measured in our experiment (Fig. 4) with existing data in the literature. In particular, Giret et al. [31] developed a thermodynamical model wherein most of the parameters were obtained from ab initio calculations. The agreement between their result and experimental measurements of  Fig. 4 of the manuscript). The sample temperature before excitation is 0 and 300 K for theory and experiment, respectively. time-resolved x-ray diffraction [14] is remarkable -not only in the general trend, but in the predicted magnitudes of motion. It is therefore a perfect benchmark to compare our results with. Figure 12 shows the inter-atomic Bi-Bi distance as a function of the pump-probe delay from Giret et al. [31] and from our measurements. The value before excitation is different because the initial temperature was 0 and 300 K in theory and experiment, respectively. While the calculations were not performed at the exact same pump fluences as our experiments, our data intercalates well between the theoretical curves. The agreement is striking: the chirped oscillations are well reproduced with consistent fluence-dependent frequencies, and the amplitude of displacement conforms to the ones predicted by theory. The only slight difference arises at very small time delays, which is to be expected for at least two reasons: (1) the pump pulse durations are different (70 fs for theory compared to 5 fs in our case) and, (2) as stated in the main manuscript, our description of carrier dynamics at small time delays is most likely flawed before they thermalize. This could impact the decomposition procedure and in turn influence the extracted atomic motion. In conclusion, this direct comparison suggests that broadband XUV transient reflectivity provides quantitative and accurate measurements of real-space atomic motion.