Unveiling quasiparticle dynamics of topological insulators through Bayesian modelling

Quasiparticle - a key concept to describe interacting particles - characterizes electron-electron interaction in metals (Fermi liquid) and electron pairing in superconductors. While this concept essentially relies on the simplification of hard-to-solve many-body problem into one-particle picture and residual effects, a difficulty in disentangling many-body effects from experimental quasiparticle signature sometimes hinders unveiling intrinsic low-energy dynamics, as highlighted by the fierce controversy on the origin of Dirac-band anomaly in graphene and dispersion kink in high-temperature superconductors. Here, we propose an approach to solve this fundamental problem - the Bayesian modelling of quasiparticles. We have chosen a topological insulator TlBi(S,Se)2 as a model system to formulate an inverse problem of quasiparticle spectra with semiparametric Bayesian analysis, and successfully extracted one-particle and many-body characteristics, i.e. the intrinsic energy gap and unusual lifetime in Dirac-quasiparticle bands. Our approach is widely applicable to clarify the quasiparticle dynamics of quantum materials. Extracting quasiparticle dispersion from photoemission data is challenging and often results in confusion when investigating low-energy excitations. As a solution the authors demonstrate a technique, which applies Bayesian analysis to extract the quasiparticle dynamics of a topological insulator from angle resolved photoemission spectroscopy (ARPES) data, and could be applied to other quantum materials.

L ow-energy excitations in interacting electronic systems are known to be characterized by quasiparticles. The concept of quasiparticle was originally proposed in the Laudau's Fermiliquid theory wherein strongly interacting electrons share a similar behavior with weakly interacting counterparts. Instead of strongly interacting bare electrons (holes), one can define dressed electrons (holes) as elementary excitations (i.e., quasiparticles), which can be understood by extending the framework of singleparticle approximation. While the Fermi-liquid theory successfully captured low-energy dynamics of normal metals and 3 He, the quasiparticle concept is nowadays applied widely in solids such as electron systems interacting with lattice vibrations (phonons) and spin excitations (magnons). Angle-resolved photoemission spectroscopy (ARPES) has played a pivotal role in uncovering key quasiparticle properties by capturing the energy dispersion (E-k relation) and lifetime of, e.g., Bogoliubov quasiparticles associated with the superconducting Cooper pairing in high-temperature superconductors 1-3 and mass-renormalized quasiparticles caused by strong electron-phonon coupling on metal surfaces and quasi-two-dimensional (quasi-2D) materials [4][5][6] . As highlighted by these examples, for the understanding of the origin and mechanism of exotic physical properties of novel materials, it is crucial to experimentally establish the nature of quasiparticles.
To elucidate the quasiparticle dynamics, it is desirable to be able to unambiguously extract the original single-particle band dispersion (bare-band dispersion, E k ) and many-body effects (self-energy, Σ) from the ARPES data. Both of these physical quantities are directly linked to the ARPES spectrum through the spectral function expressed as, where ω is the energy with respect to the Fermi level (E F ). Many attempts have been hitherto made to extract the intrinsic Σ from ARPES data by assuming a reasonable shape of E k . For example, E k was referenced to the band calculation obtained with the local density approximation (e.g., refs. 7,8 ), or it was empirically approximated with a polynomial function (e.g., linear or parabola 9,10 ). While such data analysis certainly gave insights into the quasiparticle dynamics, one often faced a serious problem in clarifying the nature of many-body interactions. This is represented by the fierce debates on the absence or appearance of an intrinsic energy gap at the Dirac point in epitaxial single-layer graphene 11,12 , which is critical for feasible application of graphene as a semiconductor device. Also, the origin of dispersion kink in cuprate superconductors (phononic, magnetic, or others) is controversial for more than a decade 13 , and its relationship with the high-T c mechanism is yet to be clarified. These controversies partially originate from a few assumptions one had to make to extract E k and Σ.
To overcome these problems, we apply semiparametric Bayesian modeling to ARPES data. We start by introducing the basics of Bayesian analysis through a simple demonstration. As a prototypical example, we model the ARPES intensity of a topological insulator (TI) TlBi(S,Se) 2 14 based on the parametric form of bareband dispersion and nonparametric forms of any other elements to perform the semiparametric Bayesian analysis of spectral function. We provide a clear insight into one-particle and manybody characteristics of TlBi(S,Se) 2 by successfully extracting bareband dispersion and self-energy.

Results
Basics of Bayesian analysis. First, we explain the basic concept of Bayesian analysis by showing its application to an energy distribution curve (EDC) contributed by multiple bands. A common approach in extracting the peak positions (i.e., contributing energy bands) is to find a good reproduction of the experimental EDCs by simulated EDCs using the least-square method. However, this method cannot pin down which class of model (e.g., how many bands are contributing) is the best for given data, often posing a question on the basic applicability of the model itself. To demonstrate this problem, we show in Fig. 1a a representative experimental EDC (dots) together with the result of numerical fittings (red curves) using the least-square method assuming the existence of intrinsic single, double, and triple Lorentzian peaks (blue curves; the number of peaks K = 1-3) that represent three different class of models. One can immediately recognize that the model with triple Lorentzian peaks shows the best fit to the experimental data. This is natural because the inclusion of more peaks (more parameters) always leads to a decrease in mean square error (shown by red line in Fig. 1b). However, it does not validate that the actual number of peaks are three. One needs to select the most appropriate model (in this case, the number of peaks) to suitably reproduce the EDC. Importantly, such selection should not include arbitrariness and must not rely on "human eyes." The Bayesian framework, as an extension of the least squares, enables the evaluation of the model's appropriateness itself in terms of the posterior probability derived from the chain rule of probability, called Bayes' formula (see Methods). One can see from Fig. 1b that the double-peak model (K = 2) has the highest probability (77%) among K = 1-5. Whereas the least-square method determines a unique solution of the parameter set as a global minimum of the mean square error, the Bayesian framework treats the "statistical ensemble" of numerous solutions as random variables with the Boltzmann distribution, called the posterior probability distribution of the parameter set. Namely, many solutions for peak position and peak width are plotted in the parameter space and colored with the posterior probability density proportional to the Boltzmann factor, where the point with the highest posterior probability density corresponds to the best-fit parameters (least-square solution) as highlighted in Fig. 1c [note that "marginal" posterior probability density is plotted in Fig. 1c, because the intensity of each peak (not shown) is also a model parameter]. Like the canonical ensemble, a statistical ensemble of Bayesian analysis is conditioned on three principal factors: the number of observed data points (e.g., the number of meshes in an ARPES image), the number of parameters in each model class, and the "temperature." This temperature is nothing to do with actual temperature of the sample and is related to the signal-to-noise ratio of observed data, affecting the uncertainty of model parameters. If the "temperature" is treated as a hyper-parameter, it can be estimated by maximizing the "partition function," which is connected with the posterior probability of model class (Fig. 1b) through Bayes' formula (see Methods).
Semiparametric model of ARPES intensity. Now that the basic scheme of the Bayesian analysis is demonstrated for a single spectrum, we apply this analysis to the actual 2D ARPES intensity through a semiparametric modeling of quasiparticle bands. For a testbed system, we have chosen TlBi(S 1−x Se x ) 2 (x = 0.8) 14 , where a slight intensity suppression around the Dirac point is seen as shown in Fig. 2a, but the origin of such unusual gap-like behavior has been a target of intensive debates [15][16][17][18][19][20][21][22][23][24] . Since the Dirac gap is a prerequisite for realizing some exotic topological quantum phenomena [25][26][27] , it is important to establish whether the bare band E k of the Dirac-cone state is gapless (Fig. 2b) or gapped (Fig. 2c) to pin down the mechanism of unusual Dirac-band anomaly (it is noted here, if the chemical potential is situated within the gap, physical properties are mainly governed by the "quasiparticle band gap," which is a combination of the bare-band gap and the selfenergy effect). It is worth noting that such absence or appearance of the Dirac gap is also critical for the classification of magnetic-TI and axion-insulator phases, as highlighted by a fierce debate on whether or not the surface state hosts the Dirac gap associated with the time-reversal symmetry-breaking magnetic order in MnBi 2 Te 4 and related compounds (see, e.g., refs. [28][29][30][31][32][33] ). In our Bayesian analysis, we treat such distinct gapless and gapped states as different model classes and judge the validity of the models for given ARPES data. As shown in the bottom of Fig. 2a-c, we assume that E k is represented by the function E s (k) for the band index s = ±1 (+1 for upper, −1 for lower Dirac cones), parametrized by the binding energy (E B ) at the Dirac point ω DP , the band asymmetry α (this parameter is associated with the effective-mass asymmetry of bulk conduction and valence bands 34,35 ), a band parameter γ, and the half-width of band gap Δ (Δ = 0 for gapless, Δ > 0 for gapped). Our first goal is to extract the actual single-particle spectral function A s (k, ω) from the ARPES data by estimating a parameter set of {ω DP , α, γ, Δ} ({ω DP , α, γ} for Δ = 0 case) and also obtain a concrete form of self-energy Σ(k, ω) (for simplicity, we assume that Σ is k independent because the k range of interest is sufficiently small).
As highlighted in Fig. 2d, the ARPES intensity is composed not only of intrinsic spectral function A s (k, ω) for the Dirac-cone states but also of the spectral weight from other features such as bulk states and spectral background. Here all these "background" states are represented by a single ω-dependent function B(ω). We also assumed independent matrix-elements M s (k) and M B (k) for the photoelectron intensities 13 for the Dirac cone and the background, respectively, and simulated the total ARPES intensity I(k, ω) by neglecting the instrumental resolution.
Semiparametric Bayesian analysis of TlBi(S,Se) 2 . Based on the above modeling of an ARPES image, we formulate the semiparametric Bayesian analysis of the overall spectral quantities, E s (k), Σ(ω), B(ω), M s (k), and M B (k) and implement this analysis by a basic algorithm for Bayesian analyses (see Methods). First, we have validated our methodology by a demonstration using mimic ARPES images that the electronic structures are predefined as ground truths (see Supplementary Note 1 and Supplementary Figure 1). Then we have applied our methodology to the actual ARPES image containing sufficiently fine E-k mesh (111 × 111 for Fig. 2a) and have succeeded in simultaneously estimating all the spectral quantities [specifically, 559 (558) scalar variables for the gapped (gapless) case]. Now we examine whether the Dirac-cone state is gapped or not by using gapless (Δ = 0) and gapped (Δ > 0) models that take into account all the above spectral contributions in the Bayesian analysis. Such examination is an obvious advantage of the Bayesian analysis and can hardly be carried out by the standard least-square method. One can immediately recognize in the left inset of Fig. 3a that the probability for Δ = 0 is negligibly small as opposed to the case for Δ > 0 (100% within computational uncertainty), indicating that the Dirac gap is indeed realized in TlBi(S 0.2 Se 0.8 ) 2 , consistent with the previous study 14 . Then we estimated the posterior probability distribution of Δ for the gapped models, as shown by histogram in Fig. 3a, where the vertical axis corresponds to the probability density. As can be seen, estimated Δ values are sharply distributed at 44.3 ± 0.3 meV, as better visualized in the magnified view in the right inset. This suggests that the energy gap can be estimated with higher accuracy and reliability through our Bayesian analysis compared to the EDC analysis applied thus far to TlBi(S 1−x Se x ) 2 14,16 .
The histogram for Δ shown in Fig. 3a is obtained by integrating all the other (558) parameters so that the posterior probability distribution for the other parameters cannot be seen from the plot. We show in Fig. 3b the posterior probability distribution against three essential band parameters ω DP , Δ, and γ in the three-dimensional (3D) density scatter plot (note that the other 556 parameters are integrated out). One can see that the data points are sharply focused in the narrow region of the (ω DP , γ, Δ) parameter space. This is also visualized by plotting the distribution against two parameter sets, i.e., Δ and ω DP in Fig. 3c (Δ and γ in Fig. 3d) in the 2D density scatter plots, which are obtained by integrating γ (ω DP ). From these results, we obtain (ω DP , Δ, γ) = (0.397 ± 0.001 eV, 44.3 ± 0.3 meV, 3.31 ± 0.02 eV·Å) as the mean and standard deviation of distributed parameter sets. This demonstrates that the Bayesian analysis is useful not only to estimate the intrinsic band parameters from ARPES data but also to see a correlation between different band parameters; these characteristics can hardly be obtained by the conventional data analysis.
One can further confirm the validity of the band model used in our Bayesian analysis by seeing that the experimental data are very well reproduced numerically. A side-by-side comparison of the experimental ARPES image and the numerical semiparametric regression function I(k, ω) (the mean values of band parameters are used) in Fig. 3e, f signifies the almost identical intensity distribution except for a higher noise level in the experiment. Such a good matching is highlighted by the obviously weak and featureless subtracted intensity in Fig. 3g. Because all the parameters are obtained from our Bayesian analysis, now we are able to show any of A s (k, ω), M s (k), M B (k), B(ω), and Σ(ω) by a 2D intensity image (see Supplementary Note 2 and Supplementary Figure 2). As an example, we show in Fig. 3h, i spectral functions for the upper and lower Dirac cones A +1 (k, ω) and A −1 (k, ω), independently. The result signifies that the dispersion of each Dirac-cone branch is rounded around the Dirac point due to the Dirac-gap opening.
To highlight the degree of agreement between the experiment and Bayesian modeling, we plot in Fig. 3j the experimental EDC at the Γ point (blue open circles) together with the numerically fitted EDC (red solid curve) that includes the background B(ω) besides the peaks from the upper and lower Dirac cones. One can see that the experimental EDC is well reproduced by the fitting curve. The apparent difference in the energy position between the upper and lower peaks demonstrates the existence of a finite Dirac gap, as also corroborated by the extracted bare-band dispersions E +1 (k) and E −1 (k) shown by dashed curves in Fig. 3h, i. This suggests that the energy gap opens in the original bare band and further indicates that the experimental suppression of spectral weight at the Dirac point cannot be understood by assuming the strongly ω-dependent self-energy effect for the gapless Dirac cone, distinct from the case of graphene on SiC 11 .
One might expect that the application of Bayesian analysis to a single experimental EDC would be sufficient for concluding a finite Dirac gap in the bare-band dispersion. However, this is not the case because the background shape can be arbitrarily chosen for the sake of just numerically reproducing the single EDC. The analysis of 2D ARPES image itself, in which the background (and matrix element and self-energy as well) is a continuous function of k and ω, is essential for extracting the intrinsic band parameters. Also, the contribution from the ω-dependent selfenergy that causes the peak shift and asymmetry in the spectral line shape can never be captured by the analysis of single EDC.
One-body and many-body characteristics of TlBi(S,Se) 2 . Thanks to the extraction of intrinsic Dirac gap through our Bayesian analysis, we can access the intrinsic many-body interactions, which was  Fig. 2a with theoretical ARPES intensity I(k, ω). Left inset shows the plot of posterior probability for Δ = 0 and Δ > 0, where Δ denotes half-width of band gap. Right inset shows the magnified view of posterior probability against Δ around the peak. b 3D contour plot of posterior probability as a function of Δ, the binding energy at the Dirac point ω DP , and the Dirac velocity γ. c, d Intensity maps of posterior probability in the (Δ, ω DP ) and (Δ, γ) space, respectively. e ARPES intensity plot of TlBi(S 0.2 Se 0.8 ) 2 (same as Fig. 2a). f Reproduced I(k, ω) obtained from the semiparametric fittings to the experimental data. g Subtraction of e and f. h, i Extracted single-particle spectral function A s (k, ω) for the upper and lower Dirac cones. Bare bands E ±1 (k) are indicated by dashed curves. j Experimental EDC at the Γ point (blue open circles) and corresponding simulated EDC (red solid curve), which is decomposed into the upper (green curve) and lower (orange curve) Dirac cones as well as spectral background B(ω) (black curve).
COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00673-6 ARTICLE COMMUNICATIONS PHYSICS | (2021) 4:170 | https://doi.org/10.1038/s42005-021-00673-6 | www.nature.com/commsphys not possible in the previous studies owing essentially to the uncertainty in determining the bare-band dispersion. We plot in Fig. 4a the real and imaginary parts of self-energy ReΣ(ω) and |ImΣ (ω)| simultaneously extracted with Δ = 44.3 ± 0.3 meV from our Bayesian analysis for the gapped-state model, compared with those obtained from the gapless-state model where Δ is intentionally fixed to 0 meV (the invalid case; see also Fig. 4c). One can see the overall smooth ω dependence of both ReΣ and ImΣ for Δ = 44.3 ± 0.3 meV (Fig. 4a), whereas there exists an unusual hump feature around ω DP in both ReΣ and ImΣ for Δ = 0 (Fig. 4b). Such anomaly is unphysical and associated with an artifact originating from the assumption of gapless Dirac-cone state despite a finite Dirac gap. In fact, when the Δ value is properly incorporated in Fig. 4a, such anomaly disappears. One can see from the self-energy plot in Fig. 4a that |ImΣ| which reflects the quasiparticle scattering rate (inversely proportional to the quasiparticle lifetime) has a broad maximum at around ω~0.15 eV, whereas it shows a minimum at~0.4 eV, around ω DP . The lower scattering rate on approaching ω DP is reasonable when we consider the available phase space of the Diraccone states, because the phase space should monotonically increase on moving away from ω DP due to the expansion of equi-energy contour in k space, as can be seen from Fig. 4d. It is emphasized, however, that the broad hump seen in |ImΣ| cannot be understood by this argument, requiring the presence of additional scattering channel. As a possible source of this channel, we point out the bulk conduction band that has a bottom at ω CB~0 .15 eV (see Fig. 2a). When ω is located in the energy range of bulk conduction band (i.e., ω < ω CB ), the surface-bulk inter-band scattering would take place besides the intra-surface-band scattering, leading to the nonmonotonic behavior of |ImΣ| around ω CB . As shown in Fig. 4a, one can also recognize that |ReΣ| becomes maximally 40 meV, comparable to the size of Dirac gap. This suggests that the influence of self-energy effects cannot be neglected in the band dispersion; in particular, the bare-band dispersion cannot be determined by simply tracing the peak maxima of EDCs.

Discussion
The present study sheds light on the fiercely debated origin of Dirac gap in TIs. The Dirac gap of TlBi(S 1−x Se x ) 2 (as well as those seen in magnetically doped TIs) has been interpreted in terms of many different scenarios standing either on the intrinsically massive Dirac fermions or the massless ones. The former involves the hybridization between surface and interface Dirac cones 36 , hybridization with impurity bands 19 , local symmetry breaking 21 , disorder-driven topological phase transition 20 , and chemical-inhomogeneity-induced smearing of band inversion 23 . The latter based on the massless Dirac fermions can be associated with the extremely strong coupling with collective modes (including phonon 17 and plasmaron 11 ), spin dephasing 24 , exciton pairing 22 , the final-state effect, etc. The present study that applies the semiparametric Bayesian modeling to ARPES data suggests that the latter approach is unlikely to be responsible for the Dirac gap. To be more specific, taking into account of the intrinsic gap on the bare band as well as the behavior of self-energy around ω DP in Fig. 4a that can be basically explained in terms of the phase-space argument for the ordinary Dirac electrons, it is suggested that many-body effects such as the electron-electron scattering and the electron-mode coupling, as intensively discussed in strongly correlated systems like high-temperature superconductors, are not responsible for the formation of massive Dirac fermions in TlBi(S 1−x Se x ) 2 . It is thus inferred that the observed gap in TlBi(S 1−x Se x ) 2 is different from the "gap" seen in the Dirac-cone band of graphene that was suggested to be associated with the many-body interactions 11 .
Since the above consideration supports an intrinsically massive Dirac fermion in TlBi(S 1−x Se x ) 2 , it would be useful to compare the present result with the results for magnetic TIs where the Dirac gap is expected to open due to the time-reversal symmetry breaking but not due to the exotic many-body interactions. The magnitude of experimental Dirac gap in the magnetic TIs such as MnBi 2 Te 4 , e.g., refs. [28][29][30][31][32][33] , and TIs proximitized with ferromagnets is very small or even undetectable by ARPES, in contrast to the sizable Dirac-gap magnitude of 2Δ = 88.6 meV revealed by the Bayesian analysis for TlBi(S 0.2 Se 0.8 ) 2 . This result, together with the fact that TlBi(S 1−x Se x ) 2 shows no magnetic order, suggests that a possibility of local timereversal symmetry breaking due to the local magnetic order is ruled out to account for the observed Dirac gap.
The semiparametric Bayesian modeling of ARPES data proposed in this study can be widely applicable to various Diracelectron systems where the interplay among the Dirac gap, symmetry breaking, and many-body interactions is of interest, as represented by the Dirac-band anomaly in magnetic TIs, axion insulators, and graphene. Also, when the appropriate analytical form of bare band is established, the Bayesian-based approach would work effectively in a wider variety of systems characterized by the band anomaly occurring in a small energy scale, such as the spin-orbit gap due to the band inversion, the small band splitting associated with the spin-orbit coupling, and the dispersion kink due to the electron-mode coupling.  Here ω and Δ denote the binding energy and the half-width of band gap, respectively. c Bare-band dispersion for the gapped (red solid curves) and gapless (red dotted curves) Dirac cones used in a, b. d Schematic band diagram of bulk and gapped Dirac-cone bands, together with the energy slices (ω DP and ω CB ) where the self-energy exhibits a characteristic energy dependence. Each error bar of |ImΣ(ω)| denotes the posterior standard deviation, where each profile of ReΣ(ω) calculated from each profile of |ImΣ(ω)| satisfies the Kramers-Kronig relation under assumption of the particle-hole symmetry.
intensity} for each peak (parameter set), while A corresponds to the set {EDC data, the