Phonon-Mediated Quasiparticle Lifetime Renormalizations in Few-Layer Hexagonal Boron Nitride

Understanding the collective behavior of the quasiparticles in solid-state systems underpins the field of nonvolatile electronics, including the opportunity to control many-body effects for well-desired physical phenomena and their applications. Hexagonal boron nitride (hBN) is a wide-energy-bandgap semiconductor, showing immense potential as a platform for low-dimensional device heterostructures. It is an inert dielectric used for gated devices, having a negligible orbital hybridization when placed in contact with other systems. Despite its inertness, we discover a large electron mass enhancement in few-layer hBN affecting the lifetime of the π-band states. We show that the renormalization is phonon-mediated and consistent with both single- and multiple-phonon scattering events. Our findings thus unveil a so-far unknown many-body state in a wide-bandgap insulator, having important implications for devices using hBN as one of their building blocks.

flake on the graphene substrate were ascertained by LEEM and spatially resolved, surface-sensitive XPEEM measurements of the B 1s and C 1s core levels (Figs. S1c and S1d, inelastic mean-free path < 0.5 nm). Additionally, dispersive plane XPS measurements of the B 1s, C 1s, and O 1s were measured within the region of the flake to ascertain the composition of the sample. The low-energy electron reflectivity (LEER) imaging measurements (discussed in the main text) were performed in the range 0-10 eV, using an incident and coherent electron beam. The µ-LEED was collected using a 1.5 µm diameter selective area aperture placed within the region of the hBN flake.
Further photoemission measurements in both real and k-space were performed using an aberration-corrected, energy-filtered photoemission electron microscope (EF-PEEM) working at an extraction voltage 12 kV (NanoESCA III, Scienta Omicron GmbH). Real-space imaging was performed with a pass energy E P = 50 eV, using a 150 µm contrast aperture, a 0.5 mm entrance slit, and a non-monocromated Hg source (hν ≈ 5.2 eV) for photoexcitation. With these settings, the microscope had the nominal real space and energy resolutions ∆x = 35 nm and ∆E = 100 meV, respectively. The hBN bandstructure was reconstructed from a series of measured constant energy surfaces (CES, E vs. k x , k y ) obtained in the diffractive plane mode of the same momentum microscope [5,6]. The bandstructure measurements were restricted to an area within the hBN flake by introducing an iris aperture of ≈ 30×15 µm 2 . Using a He I photoexcitation source (hν = 21.22 eV), CES were acquired at E P = 25 eV with a 0.5 mm entrance slit to the energy filter, yielding energy and momentum resolutions of ∆E = 50 meV and ∆k = 0.02Å −1 , respectively. Each CES was acquired 10 meV apart and with a k-space field of view of 4.3Å −1 , i.e. spanning beyond the 1st Brillouin zone of hBN. All measurements were performed at room temperature (≈ 300 K).  Fig. S2a. Notably, their intensity diminishes at binding energies above the hBN valence band maximum as shown in Fig. S2b. Also, no energy-shifted replicas of the hBN π-band from polaron formation can be distinguished from the 2D curvature analysis [8,9].

S4. ROTATIONAL AVERAGING OF THE ELECTRONIC STRUCTURE
Before extracting anything from the measured bandstructure, the three-dimensional data cone (k x , k y , E) was replicated n = 6 times. Each copy was rotated around the energy axis, i.e. theΓ point of the BZ, by φ = n(2π/6) before adding them all together. This was done to minimize any geometrically induced differences in the photoionization matrix elements [10], and furthermore to improve the signal-to-noise ratio of the measured bands.

S5. APPROXIMATING THE NON-INTERACTING π-BAND
Energy distribution curves (EDCs) of the π-band were fitted over the relevant energy range using a Lorentzian line shape, approximating the background by an error function to accommodate the visible 'step' in intensity near the VBM. The band's peak position and linewidth (in eV) were extracted and used to estimate the Re Σ and Im Σ, respectively, using an analysis procedure similar to the ones demonstrated in Refs. 11-14. As an initial guess for the unperturbed π-band, a nearestneighbor tight-binding (TB) calculation was performed with t 1 = 2.92 eV and ∆ BN = 4.3 eV [15,16]. The initial bare band curve obtained from the TB calculation was approximated by a fifth-degree polynomial over the same energy range as the fitted experimental band. Finally, its shape was iteratively adjusted to achieve consistency between the measured, real and imaginary parts of the self-energy (Σ exp ) through the Kramers-Kronig transform [11,12].

S6. ESTIMATING THE INELASTIC BACKGROUND OF Re Σ AND Im Σ
To estimate the measured inelastic energy losses due to non-bosonic interactions, a simplistic model Σ was constructed and fitted to the measured Σ exp . The model consisted of a linear combination of contributions describing the electron-impurity (Σ el-imp ), electron-electron (Σ el-el ), and electron-phonon (Σ el-ph ) scattering. The free parameters of the model were in turn adjusted to fit the measured Im Σ exp by minimizing their root mean square (RMS) difference, and subsequently, Kramers-Kronig transformed to find the model Re Σ. Each fitting parameter uncertainty was estimated from a relative RMS increase of ±5%.
The term Im Σ el-imp was modeled as a constant broadening of 380 meV (∝ T ), estimated from the full width at half maximum (FWHM) of the fitted EDC at theK point. We note here that using the FWHM directly will overestimate Im Σ imp , as it also contains additional, constant energy broadening from the instrument's finite energy resolution. The term Im Σ el-el was modeled as a S5 logarithmically corrected, quadratic expression for a two-dimensional Fermi liquid [17,18]: Finally, Im Σ el-ph was modeled by a sum of contributions Im Σ (i) el-ph (i = 1, 2) from two separate phonon modes with different maximum energies. The density of states F (i) (ω) of each phonon mode was described using a two-dimensional and isotropic Debye model (∝ ω), having a distinct and characteristic maximum (cutoff) frequency ω i . Using Eliashberg formalism, each contribution (i) to Im Σ el-ph from the individual electron-phonon interactions was calculated as [18,19]: where n(ω, T ) and f (ω, T ) are the boson and fermion distributions, respectively. The term ph is the dimensionless mass enhancement factor for coupling to the phonon mode i. Above each cut-off frequency ω i , the corresponding α 2 F (i) (ω) = 0.

S7. ESTIMATING THE ELECTRON-PHONON COUPLING FROM Re Σ
The real counterpart to the Im Σ el-ph expression in Eq. S2 is defined as [20]: Here, α 2 F (ω) is the total phonon density of states (DOS) weighted by the effective strength α 2 of the electron-phonon coupling. Thus, it describes the phonon modes participating in the interaction and contributing to Re Σ el-ph . From α 2 F (ω), the total mass-enhancement factor λ of the electrons can be found by [21]: where ω max is the maximum observable phonon energy. If Re Σ el-ph is known, e.g., from ARPES measurements, then α 2 F (ω) can be estimated from the expression in Eq. S3 by integral inversion.
The result can then be used to predict which (combinations of) phonon modes participate, and their interaction strength λ ph with the electrons can be found from Eq. S4.
To overcome the challenge of mathematical instability posed by the direct integral inversion of Re Σ el-ph a maximum entropy method (MEM) analysis procedure was employed. A detailed outline S6 of the method and its theoretical foundations can be found in Refs. [22][23][24]. Therefore, only a brief summary is given here. The optimum approximation to the experimental Re Σ el-ph is obtained by minimizing the functional: where χ 2 is the standard deviation between the measured and calculated Re Σ el-ph . The term aS serves to prevent overfitting of the experimental data by imposing physical constraints on how the α 2 F (ω) term should look. S is the Shannon-Jaynes entropy [23]: with the constraint model m(ω) reflecting our a priori knowledge of α 2 F (ω). The term a is a coefficient in the conditional probability of α 2 F (ω) under the constraints set by m(ω) [25]. In 'classical' MEM, the probability defined by a, S and m(ω) is maximized to minimize L [22].
For the analysis of the hBN π-band, the constraint model was chosen to be: This model reflects the expected quadratic increase in DOS up to the maximum energy ℏω A of the first distinguishable acoustic mode [26]. The constant m 0 is positive, thus making m(ω) positive definite over the relevant energy range. The energy ℏω max defines the cutoff energy for the observed electron-phonon interactions.
Initially, m 0 was set to approximately the average value expected for α 2 F (ω). The energy parameter ℏω A was set to approximately the E − E VBM value where the first peak in the experimental Re Σ ph could be observed (see Fig. 4a in the main text). Finally, ℏω max was fixed at E − E VBM = 400 meV at the step-like cut-off in the experimental Re Σ ph and Im Σ ph (Fig. 3b, main text). In minimizing L, the parameters m 0 and ℏω A were iteratively adjusted, achieving the final values ℏω A = 60 meV and m 0 = 0.35. This resulted in a satisfying fit to Re Σ ph , with the constraint function recreating the main feats of α 2 F (ω) (i.e., positive definite, initial quadratic increase). Yet, it was still sufficiently structureless for an unbiased fitting of the data [23].