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 $\mathrm{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.

Low-energy excitations in interacting electronic systems are known to be characterized by quasiparticles. The concept of quasiparticle was originally proposed in the Laudau's Fermi-liquid 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 single-particle 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 quasitwo-dimensional 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- * s.tokuda.a96@m.kyushu-u.ac.jp † t-sato@arpes.phys.tohoku.ac.jp 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, A(k, ω) = 1 π − Im (k, ω) [ω − E k − Re (k, ω)] 2 + [Im (k, ω)] 2 (1) 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., [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 modelling 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 bare band 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 many-body characteristics of TlBi(S, Se) 2 by successfully extracting bare-band 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 (leastsquare solution) as highlighted in Fig. 1c [note that 'marginal' posterior probability density is plotted in Fig.   FIG. 1. Application of Bayesian analysis to the energy distribution curve. a, Numerical fittings (red curves) to a typical energy distribution curve (EDC; black dots) assuming the existence of single, double, or triple intrinsic Lorentzian peaks (blue curves). A hyper-parameter K represents the number of peaks. b, Plot of mean square error (red open circles) and posterior probability (blue rectangles) against K that supports the existence of two peaks in the EDC. c, Density scatter plot of the posterior probability distribution in the Bayesian analysis (color dots), compared with the leastsquares solution (black crosses). Inset shows the posterior probability distribution as a projection in the direction of mean square error.
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 two-dimensional (2D) ARPES intensity through a semiparametric modelling 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 self-energy 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. [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 Figs. 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 Diraccone and the background, respectively, and simulated the total ARPES intensity I(k, ω) by neglecting the instrumental resolution. Noticeably, it is not necessary to assume any particular analytical forms of Σ(ω), B(ω), M s (k), and M B (k) (note that Σ is chosen to satisfy the Kramers-Kronig relation under assumption of the particle-hole symmetry), and we could approximately treat them as vectors whose elements are values of each function at every observed points (ω, k), such as {ImΣ(ω 1 ), ImΣ(ω 2 ), · · · , ImΣ(ω n )} at {ω 1 , ω 2 , · · · , ω n }. Semiparametric Bayesian analysis of TlBi(S, Se) 2 . Based on the above modelling 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  [14]. b, c, Schematic energy dispersion of gapless and gapped Dirac-cone bands, respectively. The full energy gap is 2∆. Analytical form of the bare-band dispersion Es(k) used in the model is also shown at the bottom; s, ωDP , α, γ, ∆ are branch index [upper (s = +1) or lower (s = −1) Dirac cone], Dirac-point energy, parabolic dispersion term, Dirac velocity, and Dirac gap, respectively. d, Example of simulated ARPES intensity I(k, ω) used in the semiparametric Bayesian analysis, which is composed of photoelectron matrix-element term for the surface band Ms(k), single-particle spectral function As(k, ω), photoelectron matrix-element term for the background MB(k), and angle-integrated-type background B(ω). F (ω) denotes the Fermi-Dirac distribution function.
are predefined as ground truths (see Supplementary Note 1 and Fig. S1). 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 leastsquare 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  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 As(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).
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 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 γ (∆). 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 Figs. 3e and 3f 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 Fig. S2). As an example, we show in Figs. 3h and 3i 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 modelling, 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 Figs. 3h and 3i. 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 self-energy 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 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 selfenergy, 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 Dirac-cone states, because the phase space should monotonically increase on moving away from ω DP due to the expansion of equienergy 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 which 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 surfacebulk 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 bareband 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-inhomogeneityinduced 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], and the final-state effect etc. The present study that applies the semiparametric Bayesian modelling 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 which can be basically explained in terms of the phase-space argument for the ordinary Dirac electrons, it is suggested that manybody 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 topological insulators 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 topological insulators such as MnBi 2 Te 4 (e.g., [28][29][30][31][32][33]) and topological insulators 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 time-reversal-symmetry breaking due to the local magnetic order is ruled out to account for the observed Dirac gap.
The semiparametric Bayesian modelling of ARPES data proposed in this study can be widely applicable to various Dirac-electron systems where the interplay among the Dirac gap, symmetry breaking, and manybody 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.

METHODS
Bayes' formula. Throughout our analyses, the chain rule of probability Pr(B|A) = Pr(A|B)Pr(B)/Pr(A) for random variables A and B, called Bayes' formula, was utilized. For the parameter estimation in the EDC analysis, B corresponds to the set width, position, intensity for each peak (parameter set), while A corresponds to the set {EDC data, the model class (peak number K), "temperature"} (note that the inset to Fig. 1c represents this Pr(B|A), referred to as posterior probability distribution of the parameter set). Under the Bayes' formula, what one needs to carry out is the modelling of Pr(A|B) and Pr(A), called here the likelihood function and prior probability distribution, respectively. Once these models are formulated, their appropriateness can also be evaluated by the Bayes' formula with the relation Pr(A) = Pr(A|B)Pr(B). For the model selection in the EDC analysis, the relation Pr(K|EDCdata) = Pr(EDCdata|K)Pr(K)/Pr(EDCdata) was utilized, where Pr(K|EDCdata) is referred to as posterior probability distribution of parameter set (Fig. 1b). Note that Pr(K|EDCdata) is derived by integrating out "temperature" in Pr(A) since Pr(A) = Pr(EDCdata, "temperature"|K)Pr(K) holds.
Bayesian analysis of EDC. In the EDC analysis (Fig.  1), the posterior probability distribution for the parameter set was formulated by p ∝ exp(−nβMSE)φ, where n, β, φ, and MSE are number of data points constructing EDC, a hyper-parameter ("inverse temperature"), the prior probability distribution, and the mean square error, respectively. The MSE for each K is defined by a difference between EDC data and sum of all Lorentzian functions. The function φ was set as the continuous uniform distribution whose support is [0, 0.1] (eV) for the peak width, [−0.5, −0.3] (eV) for the peak position, and [0, 1] (a.u.) for the peak intensity.

Formulation of semiparametric Bayesian analysis.
In the analysis of 2D ARPES image (Figs. 3 and 4), we assumed that the intensity Y ij of ARPES image at each pixel (k i , ω j ) for i = 1, · · · , m and j = 1, · · · , n is given by Y ij = I(k i , ω j ; w 0 , E k ) + ξ ij , where I is the intensity function defined by an equation in Fig. 2d, E k analytical form of bare-band dispersion, and w 0 the set of other elements including parameters of bare-band dispersion, the self-energy, the matrix elements, and the "background". The random variable ξ ij is observation noise subject to the Gaussian distribution whose mean and variance are 0 and β −1 0 > 0, respectively. In other words, Y ij is assumed to be subject to the conditional probability density function p (y ij | k i , ω j ; w, E s , β) with w = w 0 , E s = E k , and β = β 0 . Since w 0 , E k , and β 0 are unknown in practice and should be estimated, we treat them as random elements w, E s , and β subject to the posterior probability density function where D mn = {Y ij , k i , ω j } is a data set of ARPES image, φ(w) an arbitrary prior probability density function, and Z(E s , β) := p({y ij }|{k i }, {ω j }, E s , β) the partition function. Note that w consists of 2n + 3m + 4 (or 2n + 3m + 3) scalar values for the gapped (or gapless) state: the binding energy at the Dirac point ω DP , the band asymmetry α, a band parameter γ, the halfwidth of band gap ∆ (∆ = 0 for the gapless state), the imaginary part of self-energy ImΣ(ω j )}, the matrix elements {M +1 (k i ), M −1 (k i ), M B (k i )}, and the "background" {B(ω j )}. The function φ was set as follows: the exponential distribution whose mean is 10 (eV·Å 2 ) for α, 4 (eV·Å) for γ, 0.1 (eV) for |ImΣ|, and 0.2 (a.u.) for M s (k), the continuous uniform distribution whose support is [0, 0 The Bayesian analysis treats the statistical ensemble of w subject to p(w|D mn , E s , β) as an extension of the leastsquares method. The mean and standard deviation of p(w|D mn , E s , β) is respectively adopted as estimator and its error bar. We also estimated E s and β by treating them as random elements subject to the conditional probability distribution function where {E s } is a collection of candidate forms of E s . Note that this equation is derived from Bayes' formula such that p(E s , β) is an uniform distribution. Especially, E s and β that maximize p(E s , β|D mn ) are adopted as estimators. This type of estimators is known as the empirical Bayes estimator [37][38][39]. We also quantify the uncertainty of each E s by the marginal probability as shown in the inset of Fig. 3a (see also Fig. 1b).
Algorithm of semiparametric Bayesian analysis. The computation of p(w|D mn , E s , β) was performed by the exchange Monte Carlo method [40,41] (see also Table S1), where β is discretized as 128 points consisting of 0 and 127 logarithmically spaced points in the interval [1.5 × 10 −10 , 1.5 × 10 2 ]. The total Monte Carlo sweeps were 10, 000 after the burn-in, where the obtained sequence {w t l } for t = 1, · · · , 10, 000 and l = 1, · · · , 128 is regarded as a statistical ensemble of w subject to p(w|D mn , E s , β l ). Figures 1c and 3a-3d are the density scatter plots of {w t l } at β that maximize Z(E s , β). We also calculate p(E s , β|D mn ) via the bridge sampling [42,43], as shown by where Q β denotes the average of an arbitrary quantity Q over p(w|D mn , E s , β) and is approximated by sample mean of obtained sequence {Q t l }.

DATA AVAILABILITY
The data and information within this paper are available from the corresponding authors upon request.

CODE AVAILABILITY
The computer code to generate the results are available from the corresponding authors upon request.