Hidden self-energies as origin of cuprate superconductivity revealed by machine learning

TL;DR

This study uses machine learning to uncover hidden self-energy structures in cuprate superconductors, revealing their role in high-temperature superconductivity and demonstrating the potential of ML-based spectroscopy for fundamental insights.

Contribution

The paper introduces a machine learning approach to extract hidden self-energies from experimental spectra, revealing their critical role in cuprate superconductivity, which was previously inaccessible.

Findings

Peak structures in self-energies are linked to superconducting gaps.

Hidden self-energy peaks cancel in total self-energy, making them invisible.

Universal carrier relaxation correlates with high critical temperatures.

Abstract

Experimental data are the source of understanding matter. However, measurable quantities are limited and theoretically important quantities are sometimes hidden. Nonetheless, recent progress of machine-learning techniques opens possibilities of exposing them only from available experimental data. In this paper, after establishing the reliability of the method in various careful benchmark tests, the Boltzmann-machine method is applied to the angle-resolved photoemission spectroscopy spectra of cuprate high temperature superconductors, Bi$_2$Sr$_2$CuO$_{6+\delta}$ (Bi2201) and Bi$_2$Sr$_2$CaCuO$_{8+\delta}$ (Bi2212). We find prominent peak structures both in normal and anomalous self-energies, but they cancel in the total self-energy making the structure apparently invisible, while the peaks make universally dominant contributions to superconducting gap, hence evidencing the signal that generates the high-$T_{\rm c}$ superconductivity. The relation between superfluid density and critical temperature supports involvement of universal carrier relaxation associated with dissipative strange metals, where enhanced superconductivity is promoted by entangled quantum-soup nature of the cuprates. The present achievement opens avenues for innovative machine-learning spectroscopy method to reveal fundamental properties hidden in direct experimental accesses.

Figures

D of Ref. [13](#bib.bib13) around -50 meV is interpreted as a remnant of the original peak structure in the genuine self-energies, which is consistent with our present results. Therefore, the difference between our results and those of Ref. [13](#bib.bib13) is attributed to the artificial reduction in the amplitude of ${\rm Im}\Sigma^{\rm nor}$ due to the assumption ${\rm Re}N(\omega)=\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)$ in Ref. [13](#bib.bib13).

D of Ref. [13](#bib.bib13) around -50 meV is interpreted as a remnant of the original peak structure in the genuine self-energies, which is consistent with our present results. Therefore, the difference between our results and those of Ref. [13](#bib.bib13) is attributed to the artificial reduction in the amplitude of ${\rm Im}\Sigma^{\rm nor}$ due to the assumption ${\rm Re}N(\omega)=\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)$ in Ref. [13](#bib.bib13).

Tables

Table 1. Table 1: Doping dependence of superconducting order parameters and quasiparticle gap. Physical quantities are calculated by the self-energy inferred from ARPES data of three Bi2201 samples and Bi2212 at the optimum doping  [ 12 ] . The doping dependence of the superconducting order parameter F ​ ( k ) 𝐹 𝑘 F(k) , gap amplitude determined from the peak position in EDC, the superconducting gap amplitude Δ 0 ​ ( k ) subscript Δ 0 𝑘 \Delta_{0}(k) obtained by the present Boltzmann machine learning (BML), the quasiparticle renormalization factor z qp subscript 𝑧 qp z_{\rm qp} , c 1 ​ ( k ) subscript 𝑐 1 𝑘 c_{1}(k) , Q ¯ ¯ 𝑄 \overline{Q} , W PEAK ¯ ¯ subscript 𝑊 PEAK \overline{W_{\rm PEAK}} , ω PEAK subscript 𝜔 PEAK \omega_{\rm PEAK} and Ω 0 subscript Ω 0 \Omega_{0} are shown mainly at k = k AN 𝑘 subscript 𝑘 AN k=k_{\rm AN} , where the doping p 𝑝 p is estimated by doping dependence of T c subscript 𝑇 c T_{\rm c} in Ref.  87 . The average over the Fermi surface of F 𝐹 F , F ​ ( k F ) ¯ ¯ 𝐹 subscript 𝑘 F \overline{F(k_{\rm F})} is also estimated.

Sample		Doping	Order Parameter	Gap from EDC	Gap from BML
	$T_c$ (K)	$p$	$2F\left(k_{\mathrm{AN}}\right)\left(2\overline{F\left(k_F\right)}\right)$	(meV)	${\Delta }_0\left(k_{\mathrm{AN}}\right)$ (meV)
Bi2201 UD	23	0.12	0.085 (0.075)	60	18$\pm 1$
Bi2201 OP	35	0.15	0.152 (0.126)	30	24.1 $\pm$ 0.2
Bi2201 OD	29	0.18	0.108 (0.089)	12	11 $\pm 2$
Bi2212 OP	90	-	0.13 (-)	33	30

Equations

$$\displaystyle\hat{G}(k,\zeta)=\left[\begin{array}[]{cc}\zeta-\epsilon_{k}-\Sigma^{\rm nor}(k,\zeta)&-\Sigma^{\rm ano}(k,\zeta)\\ -\Sigma^{\rm ano}(k,\zeta)&\zeta+\epsilon_{k}+\Sigma^{\rm nor}(k,-\zeta)^{\ast}\\ \end{array}\right]^{-1},$$

$$A(k,\omega )=-\frac{1}{\pi }\mathrm{Im}{\left[\{\hat{G}(k,\zeta ){\}}_{11}\right]}_{\delta \to +0},$$

$$G^{\mathrm{nor}}(k,\omega )\equiv \hat{G}(k,\omega {)}_{11}=[\omega -{ϵ}_k-{\Sigma }^{\mathrm{tot}}(k,\omega ){]}^{-1}.$$

$${\Sigma }^{\mathrm{tot}}(k,\omega )$$

$$W(k,\omega )={\Sigma }^{\mathrm{ano}}(k,\omega {)}^2/[\omega +{ϵ}_k+{\Sigma }^{\mathrm{nor}}(k,-\omega {)}^{\ast }].$$

$$\Delta (k,\omega )=Q(k,\omega ){\Sigma }^{\mathrm{ano}}(k,\omega ),$$

$$Q(k,\omega )={\frac{1}{1-[{\Sigma }^{\mathrm{nor}}(k,\zeta )-{\Sigma }^{\mathrm{nor}}(k,-\zeta {)}^{\ast }]/(2\zeta )}|}_{\delta \to +0}.$$

$$F(k)$$

$$z_{\mathrm{qp}}^{-1}(k)$$

$${\Delta }_0(k)$$

$$\mathrm{Im}{\Sigma }^{\mathrm{nor}}$$

$$\mathrm{Im}{\Sigma }^{\mathrm{ano}}$$

$${\mbox{\boldmath$\sigma$}\equiv(\sigma_{0},\sigma_{1},\cdots,\sigma_{L-1}),}$$

$${I(\mbox{\boldmath$\sigma$})=\sum_{\ell=0}^{L-1}\sigma_{\ell}\cdot 2^{\ell}.}$$

$$\displaystyle\Theta^{L}_{\mbox{\boldmath$\sigma$}}\left(\omega\right)=\left\{\begin{array}[]{cl}1&{\rm for}\ x\in\left[I(\mbox{\boldmath$\sigma$})/2^{L},\{1+I(\mbox{\boldmath$\sigma$})\}/2^{L}\right)\\ 0&{\rm otherwise}\\ \end{array}\right..$$

$$\mathrm{Im}{\Sigma }^{\mathrm{nor}}(\omega )$$

$$\mathrm{Im}{\Sigma }^{\mathrm{ano}}(\omega )$$

$${\mathcal{B}(\mbox{\boldmath$\nu$}|\mbox{\boldmath$W$},\mbox{\boldmath$B$})=\exp\left[\sum_{\ell,m}W_{\ell m}\nu_{\ell}\nu_{m}+\sum_{\ell}B_{\ell}\nu_{\ell}\right],}$$

$${\mathcal{C}({\mbox{\boldmath$S$}})=\sum_{\mbox{\boldmath$h$}}\mathcal{B}_{C}({\mbox{\boldmath$S$},\mbox{\boldmath$h$}}|\mbox{\boldmath$W$},\mbox{\boldmath$B$}),}$$

$${\displaystyle\mathcal{C}({\mbox{\boldmath$S$}})=e^{b}\prod_{m=0}^{L_{\rm h}-1}2\cosh\left[{S_{\ell}}W_{\ell m}\right],}$$

$$\mathrm{Re}{\Sigma }^{\mathrm{ano}}(\omega )$$

$$\mathrm{Im}{\Sigma }^{\mathrm{ano}}(\omega )$$

$${{\color[rgb]{0,0,0}\mathcal{D}({\mbox{\boldmath$S$}})=\sum_{\lambda=1}^{M}w_{\lambda}\mathcal{B}_{D}(\mbox{\boldmath$S$}|\mbox{\boldmath$W$}^{\lambda},\mbox{\boldmath$B$}^{\lambda})},}$$

$${\mathcal{B}_{D}({\mbox{\boldmath$S$}}|\mbox{\boldmath$W$}^{\lambda},\mbox{\boldmath$B$}^{\lambda})=e^{\sum_{\ell,m}S_{\ell}W_{\ell m}^{\lambda}S_{m}+\sum_{\ell}S_{\ell}b_{\ell}^{\lambda}}.}$$

$$\mathrm{Re}{\Sigma }^{\mathrm{nor}}(\vec{k},\omega )$$

$$\mathrm{Re}{\Sigma }^{\mathrm{ano}}(\vec{k},\omega )$$

$$\mathcal{P}\int d{\omega }^{\prime }\frac{f({\omega }^{\prime })}{{\omega }^{\prime }-\omega }$$

$$\mathrm{Re}\int d{\omega }^{\prime }\frac{f({\omega }^{\prime })}{{\omega }^{\prime }+i\delta -\omega }.$$

$${\chi }^2=\frac{1}{2N_{\mathrm{d}}}\sum _j{\left\{A^{\exp }({\omega }_j)-f({\omega }_j)A({\omega }_j)\right\}}^2,$$

$${\sigma }_n^2=N_{\mathrm{d}}^{-1}\sum _j(A^{\mathrm{fit}}({\omega }_j)-A^{\exp }({\omega }_j){)}^2.$$

$$\overline{{\chi }^2}=\frac{1}{N_{\mathrm{d}}{N}_{\mathrm{r}}}\sum _{s=1}^{N_{\mathrm{d}}}\sum _{r=1}^{N_{\mathrm{r}}}\{A_r^{\mathrm{syn}}({\omega }_s^{(r)})-f({\omega }_s^{(r)})A({\omega }_s^{(r)}){\}}^2,$$

$$\Delta (\omega )$$

$$\left[1-Z(\omega )\right]\omega$$

$$\displaystyle K_{\pm}(\omega,\omega^{\prime})={\color[rgb]{0,0,0}g_{\rm el\mathchar 45\relax ph}^{2}}\left[\frac{1}{\omega^{\prime}+\omega+\Omega-i\eta}\pm\frac{1}{\omega^{\prime}-\omega+\Omega-i\eta}\right],$$

Full text

I Introduction

Momentum $k$ and energy $\omega$ dependent electron single-particle spectral function $A(k,\omega)$ can be measured with recent revolutionarily refined resolution of angle resolved photoemission spectroscopy (ARPES)[1]. From $A(k,\omega)$, the interaction effects crucial for unconventional superconductors can be identified in the self-energy [2, 3]. Scanning tunneling microscope (STM) and its spectra (STS) including the quasiparticle interference method [4] also give us insights into the self-energy [5, 6].

In superconductors, the self-energy consists of the normal and anomalous (superconducting) contributions, $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$, respectively. ARPES and STS provide us with only the total self-energy $\Sigma^{\rm tot}$ in a specific combination of these two [2] (see below for details). However, to understand the superconducting mechanism, it is crucially important to extract these two separately, because they represent different part of interaction effects: $\Sigma^{\rm ano}$ is proportional to the superconducting gap function, at the heart of superconducting properties, while normal-electron correlation effects, such as renormalized mass and life time, are encoded in $\Sigma^{\rm nor}$. Despite its importance, $\Sigma^{\rm ano}$ can be straightforwardly extracted separately only when $\Sigma^{\rm nor}$ is non-singular as in the BCS superconductivity of weakly correlated systems [5, 6]. Indeed, the decomposition of the self-energies of the BCS superconductors has played the role of establishing the phonon mechanism since the anomalous part contains the information of the phonon density of states, which is crucial for the identification of the glue for the superconductivity. In case of the cuprate high-$T_{\rm c}$ superconductors, because of the strong electron correlation, the subject of extracting normal and anomalous self-energy separately belongs to an open enigmatic inverse problem, which has hampered the full understanding of the superconducting mechanism for decades.

Recently machine learning and data science are developing rapidly as tools to analyze accumulated data across various research domains. The expectation is to solve complex, important problems for human being, which are hardly solved by conventional tools, as in the attempt of forecasting future and designing new functional materials by utilizing existing big data. More specifically, machine learning has potential for innovating routes of exposing quantities, which is invisible in direct measurements. Solving inverse problems to construct theories with predictive power by using existing experimental data is a typical target of the machine learning innovation.

In this paper, we develop a scheme of machine-learning technique to extract physical quantities hidden in experimental data. To demonstrate the power of our method, we apply it to the electronic structure of the cuprate high-temperature superconductors under strong correlation effects manifested by the formation of the pseudogap in the normal state. Specifically, the Boltzmann machine [7, 8] is examined to extract $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ separately from available ARPES spectra even when the normal self-energy is subject to prominent or singular correlation effects. Discovered prominent peak structure in the energy dependence of $\Sigma^{\rm ano}$ hidden in the ARPES is shown to generate most of the superconducting gap, namely more than 90% of the gap, and hence to be the driving force of the superconductivity in the cuprates. From the extracted self-energies, we elucidate the factors that determine the superconducting transition temperatures.

The organization of the present paper is as follows: In Sec. II, we summarize the essence of the regression scheme to extract the self-energies from the spectral functions. The details of the method is given in Sec. III. The readers who are not interested in the technical details of the method but are interested only in the results may skip Sec. III and directly jump to Sec. IV. Various benchmark tests of the regression by using the present method are presented in Sec. IV for simple metals and conventional Bardeen-Cooper-Schrieffer-type superconductors to show the reliability of the method. The main results of the present paper studying the regression for the cuprate superconductors are given in Sec. V. The self-energies of cuprate superconductors extracted separately for the normal and anomalous contributions are shown, which have prominent peak structures each in the normal and anomalous part, but cancels in the spectral function. We also show that the prominent peak in the anomalous self-energy gives rise to the major part of the high temperature superconductivity. Sec. VI is devoted to the discussion on the implications of the present results. Then, we summarize this paper and give our outlook in Sec. VII.

II Methods

In this section, we introduce the fundamentals and basic concept of the present method. The detailed regression scheme is given in Sec. III for readers who are interested in technical details.

II.1 Regression

The present machine learning scheme is classified into the category of a general regression task, which optimizes a function $A$ to find $B$ and/or $F$, when $A$ is a nonlinear and complex functional of another function $B$ as $A=F(B)$, where the goal of solving the inverse problem $B=F^{-1}(A)$ by the optimization of $A$ is set from a physical purpose [9]. Another example belonging to the regression task as an application of the machine learning is found in the use for quantum many-body problems or statistical physics problems (see for instance, Ref. 10).

In the present case, $A$ is the ARPES spectral function $A(k,\omega)$, $B$ is the self-energies $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ and $F$ is given by Eqs. (5) to (7) below. The training data is given by experimental $A(k,\omega)$ at a discrete and limited number of $\omega$. Then the present machine learning is a typical regression task to infer both $\Sigma^{\rm ano}(k,\omega)$ and $\Sigma^{\rm nor}(k,\omega)$ separately as continuous functions of $\omega$.

Here, we note that the present scheme is not a simple interpolation for functions of $\omega$. The regression with artificial neural network provides the self-energies $\Sigma^{\rm ano}(k,\omega)$ and $\Sigma^{\rm nor}(k,\omega)$ that simultaneously satisfy the constraints from the several rigorous and physically sound prior knowledge together with the experimental data.

In the machine learning, flexible regression models that can minimize the training error without any limitation are employed. In contrast, when we employ phenomenological regression models, the constraints arising from phenomenological functions substantially increase the cost function far beyond the noise in the experimental data as demonstrated in the following discussion (see Appendix A). If the phenomenological form of the regression model is not a priori justified, the machine learning scheme should have an advantage over standard phenomenological regression schemes.

The basic procedure is illustrated in Fig. 1. More technical details are found in the following Sec. III and the robustness, accuracy and reliability of the present machine learning are shown in detail in Appendices F and D. An important advantage of the machine learning is that the approximation converges to the correct results without overfitting and bias if one increases the data point. We do not call any other approach that does not reach this systematic improvement as a machine learning.

II.2 Green function

We propose a theoretical method to extract $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ from experimentally observed spectral functions $A(k,\omega)$ of superconductors. In the superconducting phases, the single-particle retarded Green function at a given momentum $k$ as a function of frequency $\omega$ is given by a diagonal component of $2\times 2$ matrix in Nambu representation,

[TABLE]

with $\zeta=\omega+i\delta$ ($\delta$ is a small positive real number). The bare dispersion is given by $\epsilon_{k}$. $A(k,\omega)$ measurable by ARPES is related to $\hat{G}$ as

[TABLE]

with the normal component of the Green function

[TABLE]

Here, the total self-energy $\Sigma^{\rm tot}$ is given by [11]

[TABLE]

with $W$ given as a specific combination,

[TABLE]

The gap function

[TABLE]

which is a measure of superconducting order, is proportional to $\Sigma^{\rm ano}(k,\omega)$ with the coefficient $Q(k,\omega)$ called the frequency dependent renormalization factor defined as

[TABLE]

The $\omega\rightarrow 0$ limit of $Q$ is theoretically equivalent to the quasiparticle weight (renormalization factor) defined in Eq. (12) below as we calculate in Appendix B with the help of the procedure in Appendix C. The real part of $\Delta(k,\omega)$, ${\rm Re}\Delta(k,\omega=0)$ is nothing but the superconducting gap (see the definition of $\Delta$ in Eq. (9)). In the present report, $\delta$ is chosen to be equal to the experimental resolution as $\delta=10$ meV [12], instead of taking $\delta\rightarrow 0^{+}$.

To estimate the density of the Cooper pairs, mass renormalization, and gap amplitude from the spectral function, we define $F(k)$, $z_{\rm qp}(k)$, and $\Delta_{0}(k)$ respectively as,

[TABLE]

where $z_{\rm qp}$ is called the renormalization factor (see Appendix B). As is well known, the gap function is interpreted by the product of the Cooper pair density and the effective attractive interaction to form the Cooper pair as the mean field acting on the Cooper pair formation.

II.3 Prior knowledge

When we can measure entire $\omega$ dependence of $\hat{G}(k,\omega+i\delta)$ at a fixed $k$, we can reconstruct $\Sigma^{\rm nor}(k,\omega)$ and $\Sigma^{\rm ano}(k,\omega)$ solely without any information at momenta other than $k$, because Eqs. (6) and (7) are all diagonal in the $k$ space. However, we can not measure the entire $\omega$ dependence of the complete Green function matrix $\hat{G}(k,\omega+i\delta)$, which makes the information through the Kramers-Kronig relation incomplete. In the literature [13], it has been assumed that the momentum dependence of the normal-state spectrum at a fixed $\omega$ is a single Lorentzian curve so that $\Sigma^{\rm nor}(k,\omega)$ and $\Sigma^{\rm ano}(k,\omega)$ can be extracted separately without knowing the $\omega$ dependence at large $|\omega|$. In the present article, to overcome the lack of information at large $|\omega|$, we utilize physically sound constraints and extract $\Sigma^{\rm nor}(k,\omega)$ and $\Sigma^{\rm ano}(k,\omega)$ from experimentally observed $A(k,\omega)$ at a single fixed momentum $k$, instead of assuming specific momentum dependences of the spectra. Indeed, the present Boltzmann machine learning as detailed below successfully reproduces $\Sigma^{\rm nor}(k,\omega)$ and $\Sigma^{\rm ano}(k,\omega)$ separately from benchmark spectral functions without the momentum dependence of the spectra as shown in Appendix D.

The physical constraints employed in the present article are classified into two categories. The constraints in the first category are the structure of the Green function given in Eq. (4), the Kramers-Kronig relationship between the real and imaginary part of the self-energies, negative definiteness of ${\rm Im}\Sigma^{\rm nor}$, and odd nature of ${\rm Im}\Sigma^{\rm ano}$ as ${\rm Im}\Sigma^{\rm ano}(-\omega)=-{\rm Im}\Sigma^{\rm ano}(\omega)$. The constraint in the second category is the sparse and localized nature of the ${\rm Im}\Sigma^{\rm ano}$ along the $\omega$ axis. In the present context, the sparseness is defined as the property of ${\rm Im}\Sigma^{\rm ano}$concentrated and localized in the small frequency range around the Fermi level. It should be mentioned that the optimization procedure of the present machine learning does not explicitly impose this latter constraint because of the flexible representability of the present (restricted) Boltzmann machine [14]. The sparse nature of ${\rm Im}\Sigma^{\rm ano}$ is, however, justified a posteriori in the optimized solution, which turns out to satisfy physically reasonable sparseness although the machine learning procedure does not explicitly impose this constraint.

From physical grounds, the sparse and localized nature of ${\rm Im}\Sigma^{\rm ano}$ is a natural consequence as clarified in Ref. 15, where irrespective of the mechanism and symmetry of the realistic pairing, strong and long-range nature of Coulomb repulsion in general causes severe pair breaking at large energies and suppresses $\Sigma^{\rm ano}$ at energies far away from the Fermi level.

II.4 Overview of optimization

We represent the self-energies by using artificial neural networks in the present scheme, and optimize the neural network to reproduce the experimentally observed spectral functions. The optimization of the neural network consists of inner and outer optimization loops (see Figure 1). In the inner optimization loop, starting with given initial parameters for the neural network, all the parameters are optimized to minimize the training error (see the definition given in Eq. (32) below) by following the natural gradient method (see Appendix E for details). On the other hand, in the outer optimization loop, test errors (see Eq. (34)) are minimized by the Bayesian optimization. The updated distribution delivers the initial values for the next inner loop. The next inner loop starts again with the updated neural network parameters, while, if the number of repetitions of the outer loop already reaches an upper limit (typically less than a hundred in the present paper), the optimization is completed and the current best neural network gives the optimized self-energies. Here, we combine these inner loop and outer loop optimization to avoid the overfitting, which would be inevitable if we would use only the inner loop optimization with the training data.

III Technical Details

In this section, we show the details of the present regression scheme to extract the self-energy from ARPES data. The neural network representation of the self-energy is detailed in Sec. III.1, which is supplemented by the Kramers-Kronig relation in Sec. III.2. The procedures to optimize the self-energy are given in Sec. III.3, III.4, and III.5. You may skip this section if you are not interested in the technical details of the present method.

III.1 Wavelet analysis and Boltzmann machine representation of imaginary part of self-energy

Although high-resolution ARPES data for $A(k,\omega)$ are available in experiments, $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ are not directly given separately from $A(k,\omega)$, while if $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ and $\epsilon_{k}$ are given, $A(k,\omega)$ can be determined easily by using Eqs. (4) and (5). Therefore we need to solve an underdetermined non-linear inverse problem. To overcome the underdetermined nature of the problem, we employ physically sound constraints justifiable even in strongly correlated electron systems as prior knowledge as introduced above. By incorporating these physical constraints, we try to optimize $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ so as to reproduce experimental $A(k,\omega)$. For this purpose we employ a machine-learning method by applying a Boltzmann-machine algorithm [7]. The reliability, accuracy and robustness of the present machine learning procedure are shown in several robustness test against noise in Appendix F and benchmark tests in Appendix D.

In the retarded Green function representation, ${\rm Im}\Sigma^{\rm nor}$ is negative definite and ${\rm Im}\Sigma^{\rm ano}$ is an odd function of $\omega$. The negative definiteness and odd nature are guaranteed by expanding ${\rm Im}\Sigma^{\rm nor}$ and ${\rm Im}\Sigma^{\rm ano}$ with positive definite bases $\{\theta_{j}\}$ as

[TABLE]

where $c_{j}$ and $d_{j}$ are real coefficients. Due to the positive definiteness of the normal part, $c_{j}$ is positive.

One of the simplest basis set $\{\theta_{j}\}$ is a set of rectangular functions, which gives a step-wise representation of the self-energies. To obtain a flexible and compact representation for the coefficients $c_{j}$ and $d_{j}$, here we will combine a well-established wavelet-type representation [16, 17, 18] and Boltzmann machines as follows. The high representability of the wavelet formalism with the rectangular basis is discussed in Appendix G. The piece-wise rectangular representations of ${\rm Im}\Sigma^{\rm nor}$ and ${\rm Im}\Sigma^{\rm ano}$ are shown in Fig. 2(a) and (b), respectively.

In this fitting, the frequency range of our interest $\omega\in[-\Lambda/2,+\Lambda/2]$ is first divided into $2^{L}$ grids using an integer $L$ and assign an $L$-digit binary representation as

[TABLE]

where $\sigma_{i}={\rm mod}(I/2^{i},2)$ for the decimal representation I(\mbox{\boldmath\sigma}) in the range 0\leqq I(\mbox{\boldmath\sigma})\leqq 2^{L}-1 of the grid number coordinate;

[TABLE]

Then the unit rectangular function \Theta^{L}_{\mbox{\boldmath\sigma}}\left(x\right) is defined as

[TABLE]

The correspondence between the position of the rectangular functions and the decimal representation I(\mbox{\boldmath\sigma}) is shown in Fig. 2(c).

Then, the basis set expansions of ${\rm Im}\Sigma^{\rm nor}$ and ${\rm Im}\Sigma^{\rm nor}$ are obtained as

[TABLE]

with two sets of $2^{L}$ fitting parameters, C(\mbox{\boldmath{\sigma}}) and D(\mbox{\boldmath{\sigma}}).

In the actual calculation, we take $\Lambda\sim 0.8$ eV to efficiently fit the experimentally observed $A(k,\omega)$ confined to the range $\omega>-0.4$ eV. Note that we assume ${\rm Im}\Sigma^{\rm nor}(\omega)$ is thus restricted and nonzero in the comparable range ($\omega>-\Lambda/2$ eV, $\Lambda\sim 0.8$ eV ) to the experimentally measurable one of $A(k,\omega)$ to fit them, while we do not impose any constraint on ${\rm Im}\Sigma^{\rm ano}(\omega)$, because $\Sigma^{\rm ano}$ is restricted within this energy range as we noted above. We have confirmed in Appendix F that the result is basically unchanged even when the energy cut-off $\Lambda$ imposed on ${\rm Im}\Sigma^{\rm nor}(\omega)$ is removed by adding possible high-energy tail or peak in the range $\omega<-\Lambda/2$. This means that our main finding is not altered by the experimentally unknown shape of the spectral function at $\omega<-0.4$ eV.

Irrespective of the cut-off $\Lambda$, the Kramers-Kronig relations given below (Eqs. (30) and (30)) are applied for all the self-energies to the whole energy range ($-\infty<\omega<\infty$) in the machine learning procedure so that the optimized solution satisfies the strict causality. In fact, for the converged $A(\omega)$, the sum rule $\int_{-\infty}^{\infty}d\omega A(\omega)=1$ is satisfied.

As a compressed representation for C(\mbox{\boldmath{\sigma}}) and D(\mbox{\boldmath{\sigma}}), in this paper, we employ two different types of Boltzmann machines to enhance their representational power and determine each $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$, separately. ${\rm Im}\Sigma^{\rm nor}$ is negative definite and a widely distributed function within energy scale set by the Coulomb repulsion. Thus, we use the flexible and nonnegative restricted Boltzmann machine, which will be introduced below. However, ${\rm Im}\Sigma^{\rm ano}$ has different properties: It is sparse, which is justified a posteriori as we addressed already. Therefore, we employ a mixture distribution of the Boltzmann machine without the hidden variables to accelerate the optimization. There may still remain ambiguities in determining $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ from an observed $A(k,\omega)$. As proposed below, the ambiguities are removed by imposing physical constraints of ${\rm Im}\Sigma^{\rm ano}$.

Now we introduce the Boltzmann machine to represent C(\mbox{\boldmath\sigma}) and D(\mbox{\boldmath\sigma}) in Eqs. (22) and (23), respectively.

We first change the binary variable $\sigma$ introduced in Eq. (16) to the Ising variable \mbox{\boldmathS}=(S_{0},S_{1},...S_{L-1}) using the relation $S_{\ell}=2\sigma_{\ell}-1$ for later convenience and rewrite as {\mathcal{C}}(\mbox{\boldmathS})(=C(\mbox{\boldmath\sigma})) and {\mathcal{D}}(\mbox{\boldmathS})(=D(\mbox{\boldmath\sigma})). Then by adding hidden Ising variables \mbox{\boldmathh}=(h_{1},h_{2},...), the Boltzmann machine is generally defined as a Boltzmann weight for Ising variables \mbox{\boldmath\nu}_{\ell}=\pm 1 consisting of $S$ and $h$ in the notation \mbox{\boldmath\nu}=(\mbox{\boldmathS},\mbox{\boldmathh}) as

[TABLE]

where (\mbox{\boldmathW})_{\ell m}=W_{\ell m} represents interaction among $\nu$, and (\mbox{\boldmathB})_{\ell}=B_{\ell} represents bias fields applied to $\nu$. $W_{lm}$, and $B_{\ell}$ are variational parameters to minimize the difference between the resultant $A(k,\omega)$ and the measured spectral functions. The role of the hidden variables $h$ is to enhance the representability of $\mathcal{B}$ to approximate {\mathcal{C}}(\mbox{\boldmathS}) and {\mathcal{D}}(\mbox{\boldmathS}).

Thanks to the non-negativity of \mathcal{C}(\mbox{\boldmathS}), it is efficiently represented by the restricted Boltzmann machine (RBM) [8, 19, 20], one of the most widely used one, as represented by

[TABLE]

where $\mathcal{B}_{C}$ restricts the interactions in $\mathcal{B}$ only between visible and hidden variables in the form $W_{\ell m}S_{\ell}h_{m}$.

The advantage of the RBM is that one can analytically trace out the hidden variables $h_{m}$, leading to

[TABLE]

where $L_{\rm h}$ is the number of the hidden variables. Any $\omega$-dependent line shape in the energy range $[-\Lambda/2,\Lambda/2]$ can be flexibly represented by optimized Boltzmann-machine parameters, if they are nonnegative. $L_{\rm h}$ is typically set $L_{\rm h}=2L$ to achieve a convergence with reasonable computational costs. The restricted Boltzmann machine representation is schematically illustrated in Fig. 2(d).

For $\Sigma^{\rm ano}$, to remove the ambiguities, we impose the physically required symmetries,

[TABLE]

which can be constrained by employing the odd function Eq. (23). If ${\rm Im}\Sigma^{\rm ano}$ is sparse, namely confined in a certain range of $\omega$, ${\rm Im}\Sigma^{\rm ano}$ can be better represented by a mixture distribution, namely by a linear combination of the full Boltzmann machine in the form

[TABLE]

where $M$ is the number of the Boltzmann machines in the linear combination. Here, $\mathcal{B}_{D}$ allows only the physical variable $S$ without the hidden one $h$, but allows the interaction between $S$ as

[TABLE]

The Boltzmann machine representation is schematically illustrated in Fig. 2(e).

Note that the linear combination of the Gaussian distributions is one of the standard procedure to approximate a smooth function[21] and can be used as an initial guess of $\mathcal{D}$ (for the detailed procedure, see in Appendix E.4). Of course, the Boltzmann machine has representability far beyond the Gaussian distribution after the optimization.

Since it is sufficient to take the number of variables $S_{\ell}$ at most ${L=9}$ to fit the experiment data containing the resolution limitation, we can explicitly take the trace summation over $S_{\ell}$ for all $\ell$ (with $2^{L}$ terms) at each iteration step. Therefore the drawback of the mixture distribution of BM beyond the RBM (namely, the complexity arising from containing the interaction term between two physical variables proportional to $W^{\lambda}_{\ell m}S_{\ell}S_{m}$) is not a serious problem. We set the number of the Boltzmann machines in Eq. (28) up to 3 ($M\leq 3$) to promote the faster optimization of the imaginary part of the anomalous self-energy.

When the experimental spectral data contain small noise, the step-wise representation for the imaginary part of the self-energies may introduce a systematic increase in the test errors. To reduce the possible increased error, we introduce a piecewise-linear representation instead of the step-wise representation ${\rm Im}\Sigma^{\rm nor/ano}$ (see Appendix G). Since the estimated noise is very small for the experimental ARPES spectra of the optimally doped Bi2212 ($T_{\rm c}=90$K) at 12K, the piecewise linear representation is helpful to achieve the comparable size of test errors with the noise in the experimental data (see Appendix F for quantitative discussions).

III.2 Real part of self-energy

The real part of the retarded self-energy is obtained through the Kramers-Kronig relation as,

[TABLE]

where a broadening factor $\delta=10$meV is introduced to represent a principal value,

[TABLE]

by

[TABLE]

See Appendix H for details.

III.3 Numerical procedure for optimization

We optimize the Boltzmann machines to reproduce experimentally observed spectral functions. The optimization of the Boltzmann machines consists of an inner and outer optimization loops (see Figure 3). In the inner optimization loop, starting with given initial parameters for the Boltzmann machines, all the parameters of the Boltzmann machine are optimized to minimize the training error defined in Eq. (32) by following the natural gradient method detailed in Appendix E. On the other hand, in the outer optimization loop, a test error defined below (Eq. (34)) is minimized by the Bayesian optimization that only updates the centers of mass of the distributions in \mathcal{D}(\mbox{\boldmathS}) by fixing other BM parameters to the values obtained in the inner loop. The updated distribution delivers the initial values for the next inner loop.

To find the optimized self-energy, these inner and outer loop optimization processes are combined as follows. First, Boltzmann machines are initialized to follow the prior knowledge explained above. Then, these Boltzmann machines are optimized to minimize the training error in the inner loop. Once the inner-loop optimization converges, the test error is evaluated by the optimized Boltzmann machines, and the self-energies given by the present Boltzmann machines are stored as the current best ones. When the outer loop is repeated, only if the test error becomes minimum in the whole optimization history, the present self-energies are stored as the current best. Then, the centers of mass of the $M$ Boltzmann-machine distributions are updated in the outer loop by using the Bayesian optimization scheme. The next inner loop starts again with the updated Boltzmann machines parameters, while, if the number of repetitions of the outer loop already reaches an upper limit (typically less than a hundred in the present paper), the optimization is finalized and the current best Boltzmann machines give the optimized self-energies. For an efficient optimization, the initial condition at the largest $L$ ($L=9$) is prepared from the optimization at $L=8$ (see Appendix E.)

III.4 Minimization of training error

For given initial parameters of the Boltzmann machines, the least square error of the training defined by

[TABLE]

is minimized, where $N_{\rm d}$ is the number of the experimental data points, $A^{\rm exp}(\omega)$ is an experimentally observed $A(k,\omega)$, $\{\omega_{j}\}$ ($j=1,2,\dots,N_{\rm d}$) is the set of frequency where $A(k,\omega)$ is observed in the experiment, and $f(\omega_{j})$ is a convolution of the Fermi-Dirac distribution and a Gaussian distribution. The experimental data $A^{\rm exp}(\omega_{j})$ involves the Fermi-Dirac distribution broadened by the resolution of the experiments. Therefore, we introduce the convolution $f(\omega)$ of the Fermi-Dirac distribution at 12 K for Bi2212 and 11 K for Bi2201 and the Gaussian distribution with standard deviation $\sqrt{\sigma^{2}}=5$ meV. Here, we normalize the experimental data $A^{\rm exp}(\omega)$ by assuming $(1/N_{\rm d})\sum_{j}A^{\rm exp}(\omega_{j})=n_{0}$, where $0<n_{0}<1$. In this paper, we infer $n_{0}=0.3$ per spin, which means that 60% of an electron is assumed to be distributed in the experimentally observed range ($\omega\gtrsim-0.4$ eV). We show that $n_{0}=0.3$ is indeed the optimized value of the least square fit in Appendix I while the result of the self-energies does not sensitively depend on the choice of $n_{0}$ around 0.3.

III.5 Minimization of test error

To further explore the multi-dimensional parameter space of the Boltzmann machine and find an optimized solution, we employ the Bayesian optimization scheme in the outer loop. Before performing the Bayesian optimization, as explained above, we perform sufficiently large number of optimization steps, which is typically $4\times 10^{3}$, in the inner optimization loop to minimize the training error.

Then, we update the center of mass of each component of mixture distribution represented by Boltzmann machine in ${\rm Im}\Sigma^{\rm ano}$ defined in Eq. (23), in the following procedure. First, we extract the weight, center of mass, and variance of each Boltzmann-machine distribution by zeroth, first, and second moments as the function of $\omega$. To update the center of mass, we use the Bayesian optimization scheme [21] depending on the history of the optimization process for the center of mass, where the test error $\overline{\chi^{2}}$ defined below is the cost function to be minimized instead of $\chi^{2}$, to avoid overfitting. Then, we construct the initial values for \mathcal{D}(\mbox{\boldmathS}) that defines ${\rm Im}\Sigma^{\rm ano}$ for the next inner loop by using the updated center of mass with the weight and variance obtained above. Each Boltzmann machine in the mixture distribution \mathcal{D}(\mbox{\boldmathS}) is initialized as the Gaussian distribution with the weight, the updated center of mass, and variance. One may wonder why we initialize the Boltzmann machine again with the Gaussian distribution: Since the next inner loop optimizes the Boltzmann machine again, the initialization by the Gaussian does not alter the final results, where the final convergence to the optimized self-energies is reached after the inner loop. The reason to reduce the distribution temporarily to the Gaussian is that the outer loop optimization can be handled easily since only the three lowest moments are needed.

To define the test error, first, we generate synthetic experimental data from the original data. Because the overfitting originates from reproducing detailed noisy behaviors in the experimental data finer than the experimental resolution [22], to eliminate the noise, the experimental data $A^{\rm exp}$ is fitted by a smooth function $A^{\rm fit}$ defined as a linear combination of the Gaussian distributions [21] with standard deviation $\sqrt{\sigma^{2}}=10$ meV equal to the experimental resolution. Then, we can estimate amplitude of noise in the experimental data as

[TABLE]

This error estimate is a standard procedure in the linear regression problem [23]. Inferring error or noise of the experimental data from the smoothed curve represented by an interpolation of the experimental data based on the physical assumption of smooth and continuous behavior in nature is a general established procedure in the linear regression problem. (See for instance, Ref. 21). This is a natural regularization to infer the generalization error reliably. There, it is important to assume the smoothness within the scale of the experimental resolution and with a frequency scale sufficiently longer than the interval of the experimental discrete points to exclude overfittings. This is achieved by the superposition of the Gaussian with the 10meV width to meet the experimental resolution, 10meV (see Ref. 9 as well).

By using the probability distribution $p(A^{\rm syn}|A^{\rm fit},\omega)\propto\exp\left[-\left(A^{\rm syn}(\omega)-A^{\rm fit}(\omega)\right)^{2}/2\sigma_{n}^{2}\right]$, we can generate synthetic experimental data $A^{\rm syn}$. By assuming that $p(A^{\rm syn}|A^{\rm fit},\omega)$ well reproduces real experimental data, the cost function to avoid the overfitting is defined by,

[TABLE]

where $A^{\rm syn}_{r}$ is the $r$th synthetic experimental data independently generated by the probability distribution $p(A^{\rm syn}|A^{\rm fit},\omega)$ and $\{\omega_{s}^{(r)}\}$ is a set of randomly chosen frequency points for each synthetic data $A^{\rm syn}_{r}$. See also Appendix F.

The optimization of the internal parameters of the Boltzmann machines in the inner loop and the optimization of the center of mass of each mixture distribution in the outer loop is repeated together several tens. The self-energies that give the minimum of $\overline{\chi^{2}}$ for the test data are called the optimized self-energies.

IV Benchmark Tests

In this section, we benchmark the performance of the present self-energy inference by utilizing the artificial neural network representation. First, we reproduce the normal self-energy of metals. Since there is no anomalous component, the regression becomes easier and less underdetermined. As typical examples, the surface state of Be and normal state of (La,Sr)2CuO4 are analyzed.

Next, we perform the regression to reproduce the known normal and anomalous self-energies of superconductors. The present method indeed reproduces these self-energies for a model of strong coupling Bardeen-Cooper-Schrieffer (BCS) superconductors. The readers who are interested in the main results on the superconducting cuprates may skip this section and go to Sec. V.

IV.1 Normal metals

To benchmark the capability of the present scheme, we analyze non-superconducting metals as trivial examples.

IV.1.1 Surface state of Be

The ARPES spectrum of Be (0001) surface state [24] is analyzed as a typical metal. In Fig. 4(a), the spectral function at the Fermi momentum $k_{\rm F}$ obtained by our regression task reproduces the experimental data. The optimized normal self-energy is consistent with the self-energy given by the combination of Migdal-Eliashberg (ME) theory and data from experiments by considering a constant shift of the imaginary part due to elastic scatterings, as shown in Fig. 4(b). Note that our machine learning procedure does not assume the ME theory and shows the ability of reliable regression later even when the ME theory is not a priori justified as in the case of strongly correlated electron systems.

IV.1.2 Normal state of (La,Sr)2CuO4

We have also analyzed the ARPES data of Zhou, et al. [25], to see the normal-state self-energy of (La,Sr)2CuO4. After subtracting the background proposed in this paper, we have analyzed intrinsic spectral function by the machine learning and succeeded in reproducing the spectral function from the optimized self-energy as shown in Fig. 5.

We also compare thus obtained self-energy with the results of the standard Kramers-Kronig (KK) transformation scheme and found that they are essentially consistent each other. Here, the standard KK transformation scheme employed in Ref. 25 consists of three steps: First, a high-energy tail is attached to the spectral function to obtain the full $\omega$ dependence of $-(1/\pi){\rm Im}G(\omega)$. Second, by assuming the particle-hole symmetry, the KK transformation is used to obtain the real part of $G$. Then, from the Dyson equation, the self-energy $\Sigma=\omega-G^{-1}$ is obtained. Without the assumption on the high-energy tail and the particle-hole symmetry, our regression scheme indeed reproduces the results of the standard KK transformation.

IV.2 BCS superconductors

Here, we show that the present scheme reproduces the self-energies of the strong coupling BCS superconductors. In the BCS superconductors, the amplitude of the self-energies is smaller than those in strongly correlated electron systems such as cuprate superconductors. The weaker self-energy effects are suitable for perturbative treatments for forward problems such as the ME theory if the ME theory is justified. In contrast, the present self-energy regression from the spectral functions works well for the strongly correlated electron systems because the spectral functions provide more information necessary for solving the inverse problem to obtain the self-energies, such as large superconducting gaps and broader quasiparticle peaks. Then, the regression of the BCS superconductors are difficult tasks to perform by using the present scheme when the amplitude of the anomalous self-energy and its influence on the spectral function are small. Nevertheless, in this section, we demonstrate the successful regression of the self-energies for a BCS superconductor when it is close to the strong-coupling limit, in which we have substantial amplitude of the anomalous self-energy.

As a typical model, the superconducting state described by the following Eliashberg equations is examined. By following the standard strong coupling theory for boson-mediated superconductors [26, 27, 5], the superconducting gap $\Delta(\omega)$ and particle-hole symmetric component of the normal self-energy $\omega(1-Z(\omega))$ at zero temperature are given by the Eliashberg equations,

[TABLE]

where $\Delta_{0}={\rm Re}\Delta(\Delta_{0})$, $\omega_{\rm c}$ is the cutoff frequency, and $\eta$ is a positive broadening factor. Again, we assume that the noninteracting density of states is given by a momentum and energy independent constant $N_{F}$ for simplicity and the superconducting symmetry is momentum independent $s$-wave. In the following, we use the name “phonon” for the boson, although Eqs. (35), (36), and the following Eq. (37) can be used to describe electrons coupled to any localized optical boson mode, irrespective of the origin of the boson mode. Here, we assume that the kernel functions $K_{\pm}$ originate from the Einstein phonon as

[TABLE]

where $\Omega$ is the Einstein phonon frequency and $g_{\rm el\mathchar 45\relax ph}$ is the electron-phonon coupling constant. With the assumption that the density of states is a constant around the Fermi level, the self-energies obtained by the Eliashberg equations are independent of the electron density and dimension of the system. Then, the normal component of the Green function is defined as

[TABLE]

where $\phi(\omega)=Z(\omega)\Delta(\omega)$ and $\epsilon$ is the energy measured from the Fermi energy. The spectral function is defined as $A_{\rm b}(\omega)=-(1/\pi)f_{\rm FD}(\omega){\rm Im}G^{\rm nor}_{\rm b}(\epsilon=0,\omega)$, where $f_{\rm FD}(\omega)$ is the Fermi-Dirac distribution, and the superconducting density of states is given by

[TABLE]

The notation of the self-energies used in the literature on the phonon-mediated superconductors [26, 27], $\omega[1-Z(\omega)]$ and $\phi(\omega)$, is different from $\Sigma^{\rm nor}(\omega)$ and $\Sigma^{\rm ano}(\omega)$ in the present paper. The normal and anomalous components of the self-energies are obtained as

[TABLE]

In Fig. 6, an example of the self-energy inference for the phonon-mediated superconductors is shown. Here, we set the coupling constant {\color[rgb]{0,0,0}g_{\rm el\mathchar 45\relax ph}^{2}}N_{\rm F}=0.275 eV, the Einstein phonon frequency $\Omega=0.075$ eV, the cutoff frequency $\omega_{\rm c}=4$ eV, and the broadening factor $\eta=0.0075$ eV2. The parameters in $\Sigma^{(0)}(\omega)$ are chosen as $a=0.008$ eV2, $b=0.016$ eV2, and $\alpha=0.005$ eV3. For normalization of the spectral function, we choose $n_{0}=0.3$ without tuning. Although its effect is negligibly small in the inferred self-energies, the Fermi-Dirac distribution with $T=40$ K is introduced in the spectral function used in the machine learning just by following the scheme with finite-temperature experimental data. The machine learning results capture essential features of the original normal and anomalous self-energies that generate the target spectrum $A_{\rm b}$, where the anomalous self-energy has a dip around $-(\Delta_{0}+\Omega)$ and $-(\Delta_{0}+2\Omega)$ and the normal self-energy shows a sharp drop. The dip in the anomalous self-energy, which arises from the electron-phonon coupling and gives the superconducting gap through the Kramers-Kronig relation, is responsible for the $s$-wave superconductivity. These features are characteristic of the strong coupling BCS (phonon-mediated) superconductors and the anomaly of $N_{\rm b}(\omega)$ at $-(\Delta_{0}+\Omega)$ shown in the inset of (a) is regarded as the evidence of the phonon mechanism. Although the slight deviation between ${\rm Im}\Sigma^{\rm nor}_{\rm b}$ and the inferred self-energy ${\rm Im}\Sigma^{\rm nor}$ for $\omega\lesssim-0.3$eV due to the finite high-energy cutoff in ${\rm Im}\Sigma^{\rm nor}$, in contrast to ${\rm Im}\Sigma^{\rm nor}_{\rm b}$ that stays constant even for $\omega\ll-0.3$eV, the machine learning well reproducing the exact results of the dip in ${\rm Im}\Sigma^{\rm ano}$ and the sharp drop in ${\rm Im}\Sigma^{\rm nor}$ indicates the reliability of the present method.

V Results: Cuprate Superconductors

In this section, we will show the results of the regression for ARPES data of Bi2212 and Bi2201. The prominent peak structures are found both in the normal and anomalous self-energies, which cancel each other.

V.1 Experimental data

We utilize high resolution ARPES data taken for two cuprate compounds, Bi2212 for optimally doped sample with critical temperature $T_{\rm c}\sim 90$K [12] and Bi2201 underdoped sample with $T_{\rm c}\sim 29$K [28]. We analyze Bi2212 data at temperature $T=12$ K and Bi2201 at $T=11$K, which are both well below $T_{\rm c}$. The machine learning enables us to obtain $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ separately, and reveals prominent peak structures in both of them, which are apparently hidden in the original ARPES data, because of the cancellation of these two contributions. We elucidate its profound consequences for the superconducting mechanism.

Although tremendous efforts have been devoted since the discovery of the cuprate superconductors with many fruitful clarifications, various puzzling issues remain open. The normal-state $A(k,\omega)$ is highly unusual including the pseudogap. Nevertheless, the superconducting phase does not look so unusual except for the $d$-wave-type nodal gap itself and somewhat inconspicuous “peak-dip-hump” structure (see red square symbols in Fig. 7(a)[12, 28]): Outside the sharp quasiparticle peak (at $\sim-40$ meV in Fig. 7(a)) expected at the superconducting gap edge, $A(k,\omega)$ (energy distribution curve (EDC)) particularly at the antinodal point $k=k_{\rm AN}$ is characterized by a deeper-energy weak dip followed by a broad hump [1, 2]. In contrast, the underdoped sample does not show the gap-edge peak (Fig. 7(b)), although comparable gaps $\sim 30$ meV open as a first look. They are in contrast with the strong-coupling BCS superconductors, where the solution of the Eliashberg equation using the phonon density of states predicts prominent peaks (or saw-tooth-like) structures outside the gap in $A(k,\omega)$ (or density of states after angle integration), which has finger-print correspondence to the actual peak measured by the tunneling spectra, while the peak is shown to be crucial in the emergence of the superconductivity, thus constitutes the decisive testimony of the electron-phonon mechanism [5, 6, 11].

There are limitation on information available from the ARPES measurements. The unoccupied states invisible in the ARPES spectra may impose such a limitation. The bilayer nature of Bi2212 may also affect our regression scheme. There are also issues such as photon energy dependence and effects of matrix elements. The influence of these uncertainties and issues do not change our results in the following sections, as examined in Appendices F and J.

V.2 Prominent peaks in self-energies revealed by Boltzmann machine

By using the Boltzmann-machine learning, a dramatic consequence is revealed for $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ by reconstructing them from the mild structure of $A(k,\omega)$ given by ARPES. The present reconstruction is a non-linear underdetermined problem as in many of machine learning problems. To obtain a reliable solution, we utilize physically sound constraints such as the rigorous causality encoded as the Kramers-Kronig relation. Sparse and localized nature of ${\rm Im}\Sigma^{\rm ano}$ has resulted a posteriori as the optimized solution under physical constraint as detailed in Sec. II.3. To represent the self-energies and incorporate the physically sound constraints, the Boltzmann machine as universal function approximators developed in machine learning is employed. See Sec. III and also Figure 3 for the flow chart.

The obtained $A(k_{\rm AN},\omega)$ (cross points in Fig. 7) perfectly reproduces the distinct behaviors of both the optimally doped and underdoped samples. These EDC curves are constructed from ${\rm Im}\Sigma^{\rm nor}(k_{\rm AN},\omega)$ and ${\rm Im}\Sigma^{\rm ano}(k_{\rm AN},\omega)$ in Figs. 8(a) and (b). Remarkably, prominent peaks are found in ${\rm Im}\Sigma^{\rm ano}(k_{\rm AN},\omega)$ at $\omega=\omega^{\rm OP}_{{\rm PEAK}}\sim\pm 70$ meV for Bi2212 and at $\omega=\omega^{\rm UD}_{{\rm PEAK}}\sim\pm{45}$ meV for Bi2201, accompanied by weaker peaks at $\pm$ 180 meV and $\pm$ 160 meV, respectively. We will show below that the discovered peaks are the main source of superconductivity. Although the peak of ${\rm Im}\Sigma^{\rm ano}(k_{\rm AN},\omega)$ had been searched for long time in analogy to the strong coupling BCS superconductors [2, 13], its clear signature was missing in the cuprates. The machine learning now has succeeded in its identification. Surprisingly, ${\rm Im}\Sigma^{\rm nor}(\omega)$ also has distinct (positive or negative) peaks at the same energy as ${\rm Im}\Sigma^{\rm ano}(\omega)$ and, as we clarify below, their contributions to the spectral function cancel out each other.

The robustness of our finding against noise, uncertainty in experimental data and the experimentally uncertain high-energy part is demonstrated in Appendix F and D. In Fig. 8(a), we also plot the error bars of the peak position and height inferred from the experimental uncertainty and the machine learning error (for detailed procedure of the error-bar estimate, see Appendix F.1, Figure 17, and Sec. III (“Minimization of training error” and “Minimization of test error” as well as Appendix G). The small error bars indicate that the existence of prominent peak is reliable. For this analysis, we have utilized the inferred experimental noise estimated from the interpolated smooth $A(k,\omega)$ (namely, the inferred true smooth $A(k,\omega)$) obtained in the standard linear regression analysis . Then the error bars of our peak estimate here is given from the optimization to hypothetical experimental data points generated with the same level of noise to the inferred true $A(k,\omega)$. See Sec. III.5 for details. We note that in the case of underdoped Bi2201 sample below $T_{\rm c}$, there exists subtlety in the machine learning solution. Although the superconducting solution presented here is the optimized solution that gives the smallest mean-square error in the fitting of $A(k,\omega)$, an insulating solution is also found with larger error. This may be related to severe competition of insulating and superconducting behaviors in the real sample. We show the properties of the superconducting solution because it gives the best optimized solution and the sample is indeed superconducting.

Despite the peaks in ${\rm Im}\Sigma^{\rm nor}$, and ${\rm Im}\Sigma^{\rm ano}$, prominent peaks are missing in ${\rm Im}\Sigma^{\rm tot}(\omega)$ as shown in Figs. 8(c) and (d) (black symbols). We discuss below why the peaks in ${\rm Im}\Sigma^{\rm ano}$ at $\omega\neq 0$ necessarily show up and their contribution cancels with ${\rm Im}\Sigma^{\rm nor}$, when we impose physically justifiable constraints such as the Kramers-Kronig transformation.

Instead of the peak in ${\rm Im}\Sigma^{\rm nor}(\omega)$ and ${\rm Im}\Sigma^{\rm ano}(\omega)$, a negative prominent peak [29] generating the superconducting gap is found centered at $\omega\sim 0$ in ${\rm Im}\Sigma^{\rm tot}$, which arises from the zero of the denominator in Eq. (8), commonly to the conventional BCS superconductors.

V.3 Role of self-energy peaks to superconductivity

To understand the significance of the peaks at $\omega_{\rm PEAK}$ in ${\rm Im}\Sigma^{\rm ano}$, we show the contribution of the peaks to ${\rm Re}\Sigma^{\rm ano}(k,\omega=0)$ estimated from the normalized partial Kramers-Kronig relation [30] (Cauchy relation) (see also Eqs. (30) and (31)) defined by

[TABLE]

${\rm Re}\Sigma^{\rm ano}(k,\omega=0)$ is a measure of the superconducting amplitude, because the gap $\Delta(k,\omega=0)$ is proportional to $\Sigma^{\rm ano}(k,\omega=0)$ (Eq. (9)). Since $I_{\Sigma}(\Omega=-\infty)=1$, the contribution of the peak in ${\rm Im}\Sigma^{\rm ano}$ to the superconductivity can be estimated from the increment in $I_{\Sigma}(\Omega)$. Figure 9(a) shows that the inner energy peak at $\omega_{\rm PEAK}\sim-70$ meV ($-45$ meV) for Bi2212 (Bi2201) both contribute to more than 90% of ${\rm Re}\Sigma^{\rm ano}(k_{\rm AN},\omega=0)$ (note that ${\rm Im}\Sigma^{\rm ano}$ is an odd function of $\omega$). Namely, these peaks are indispensable in the emergence of the superconductivity.

To further demonstrate the crucial role of the inner peak in ${\rm Im}\Sigma^{\rm ano}(k_{\rm AN},\omega)$, we have hypothetically eliminated the peaks in ${\rm Im}\Sigma^{\rm nor}$ and ${\rm Im}\Sigma^{\rm ano}$ as shown in Fig. 9(c) for Bi2212. The resultant $A(k,\omega)$ in Fig. 9(d) shows the gap disappearance and switching to a normal metal with a sharp quasiparticle peak, confirming the crucial role of the peaks to superconductivity.

Through our Boltzmann-machine analyses, $\Sigma^{\rm nor}$, $\Sigma^{\rm ano}$ and $W$ are revealed to have prominent (positive or negative) peaks, while they cancel each other in the sum $\Sigma^{\rm tot}$. It is important that this conclusion is obtained directly from experimental data without assuming any theoretical models aside from mathematical (causal) requirement for the Green function. A recent self-energy analyses of ARPES data [13, 31] did not identify the present prominent structure. However, our analyses on the momentum dependence suggest that the results are not necessarily inconsistent each other for the case of Ref. 13, because the momentum range in Ref. 13 is limited in the nodal region far from the antinodal point, where we also see that the prominent peak is missing because of the $d$-wave symmetry of the superconducting gap. It is crucially important to see the antinodal region to see the prominent peak as we discuss in Appendix A and  L. Further, we show in Appendix A that physically inappropriate assumptions for strongly correlated electron systems posed in Refs. 13 and 31 lead to failures of identifying the existence of the peaks.

If a large superconducting gap is open around $\omega=0$ as in the experimental $A({k,\omega})$, it requires the corresponding famous gap structure in ${\rm Re}\Sigma^{\rm ano}$ around $\omega=0$, where inside the two peaks at the gap energy $\omega=\pm\Delta$, ${\rm Re}\Sigma^{\rm ano}$ shows plateau as shown in Fig. 8(a). Then consistency requires a prominent peak structure around $\omega_{\rm PEAK}$ in ${\rm Im}\Sigma^{\rm ano}$ through the Kramers-Kronig relation. The peak of ${\rm Im}\Sigma^{\rm ano}$ in turn naively anticipates prominent structures outside the gap in $A(k,\omega)$ through Eq. (6). However, such structures are missing. This is possible only when $\Sigma^{\rm nor}$ plays a role to cancel the prominent structure in ${\rm Im}W$. This is corroborated by the vanishing superconductivity in Fig. 9(d). Furthermore, the superconducting order accompanied by coherent quasiparticle excitations observed in experiments can be generated from ${\rm Im}\Sigma^{\rm ano}(\omega)$ only when electrons at $\omega$ become coherent, signaled by the reduction of $|{\rm Im}\Sigma^{\rm tot}|$ (or more precisely $|{\rm Im}\Sigma_{\rm N}|$) (see Figs. 8(c) and (d), their captions, Sec. III, and Appendix K for the definition of $|{\rm Im}\Sigma_{\rm N}|$) seen in the region $|\omega|<\omega^{\ast}\sim 0.07$ eV. This restricts $\omega_{\rm PEAK}$ to this range. The present machine learning indeed reproduced this natural expectation in a physically transparent and reasonable way.

We show some analyses on the temperature dependence of the self-energies including a case above $T_{\rm c}$ and the momentum dependence away from the antinodal point in Appendix M and L, respectively, as supporting data of the present analyses. The results confirm the validity of the present conclusion.

Although we found that the prominent peak in ${\rm Im}\Sigma^{\rm ano}$ is crucial for the high critical temperature of the curates, full understanding and mechanism of prominent peaks in $\Sigma^{\rm nor}(\omega)$ and $\Sigma^{\rm ano}(\omega)$, which are absent in $\Sigma^{\rm tot}(\omega)$ are open to further analyses. Our result is significant because the cancellation poses a severe constraint on theories of curate superconductors and calls for further consistent studies based on this finding. In Discussion, we refer to one possible explanation.

VI Discussion

VI.1 Insights into intrinsic (Planckian) dissipation

In this and next subsections, we discuss further possible connection gained from the main findings of the peaks and their cancellation to other more open issues of the cuprate superconductors in regard to the damping (incoherence) of quasiparticle and the factors that determines $T_{\rm c}$ to emphasize the significance of the findings with outlook. The present machine learning is also useful to separately extract other theoretically fundamental quantities such as the momentum resolved superconducting order parameter (the density of Cooper pairs or the superfluid density) $F(k)$, mass renormalization factor $z_{\rm qp}(k)$ and the single-particle relaxation time $\tau$, which had been inferred only indirectly or only in combinations of more than one quantity in experiments in the literature, although these quantities play crucial roles independently of each other below in understanding physics. (see Sec. II.2 for precise definition of the above quantities).

How frequently the single-particle excitations are scattered is encoded in the imaginary part of the normal self-energy ${\rm Im}\Sigma^{\rm nor}$. Landau’s Fermi-liquid-like behavior characterized by ${\rm Im}\Sigma^{\rm nor}(\omega)\propto\omega^{2}$, is satisfied only in a small region ($|\omega|<0.03$ eV) for Bi2212, and looks even linear ($\propto|\omega|$) in the same region for Bi2201, implying non-Fermi liquid (marginal Fermi liquid) behavior [32] (see Figs. 8(c) and (d)), which can be fit by ${\rm Im}\Sigma^{\rm nor}(k,\omega)\sim c_{0}(k)+{\rm sign}(\omega)c_{1}(k)\hbar\omega$ in the range 15 meV $<\omega<40$ meV with a dimensionless marginal-Fermi-liquid coefficient $c_{1}(k)$.

The $\omega$-linear component $c_{1}(k)\omega$ is disruptive to the quasiparticle picture, and manifests emergent inelastic dissipation absent in Landau’s Fermi liquids. As supporting information, tiny quasi-particle renormalization factor $z_{\rm qp}$ corroborating the non-Fermi liquid together with its effects on pair breaking is also shown in Appendix B (see definition of $z_{\rm qp}$ in Eq. (12) ).

The single particle relaxation time $\tau$ is defined by $\tau(k,\omega)^{-1}=z_{\rm qp}(k){\rm Im}\Sigma^{\rm nor}(k,\omega)/\hbar$. When the carrier relaxation time is estimated from $\tau$, the $\omega$-linear term, $z_{\rm qp}(k)c_{1}(k)\omega$, is associated with the universally observed $T$-linear resistivity in the cuprates [33, 34] through the $\omega$-$T$ correspondence $\tau(\hbar\omega)\leftrightarrow\tau(k_{\rm B}T)$ transformed to the self-energy of two-particle Green function for the conductivity. The temperature-insensitive $z_{\rm qp}(k)c_{1}(k)$ shown in Appendix N also supports the correspondence. (See Fig. 25 for each $z_{\rm qp}(k)$ and $c_{1}(k)$.)

A remarkable property of the inelastic relaxation rate $\Gamma(k)=z_{\rm qp}(k)c_{1}(k)$ is its high value ($\sim 1$) with only weak dependence on the doping, momentum (see Fig. 10(a)) and temperature. This universal behavior of $\Gamma\sim 1$-$1.5$ seems to support a local and universal mechanism of the relaxation, for instance, the Planckian dissipation mechanism of the hydrodynamic state, which claims $\tau^{-1}(T)=\Gamma k_{\rm B}T/\hbar$ or $\tau^{-1}(\omega)=\Gamma\omega$ with a universal constant $\Gamma$ of the order unity [35, 36].

Although simple version of Planckian mechanism expects only an extended broad self-energy structure due to “unparticle physics,” the self-energy has a broad but prominent peak structure around $\omega=\omega_{\rm PEAK}$ which is responsible for the superconductivity through the Kramers-Kronig transformation as we discussed. At the same time, the actual line shape is rather broad with the width around 0.05 eV (see Figs. 8(c) and (d)), which is comparable to $\omega_{\rm PEAK}$ itself. More importantly, the peak is smoothly connected in the tail with the $\omega$-linear behavior near the zero energy, implying that the “Planckian dissipation and hydrodynamic behavior” associated with the strange metal [35] is caused by the source of the superconductivity. The broad prominent peaks could be due to the damped pole but it could also be ascribed to “unparticle object” generated by entangled bare electron and dark object.

The marginal fermi liquid behavior (${\rm Im}\Sigma^{\rm nor}(\omega)\propto\omega$) needs to be understood with care. Since the present photoemission data could include extrinsic background effect, our analysis may not clarify the high-energy part of intrinsic $\omega$-linear behavior. In fact, in relation to the $T$-linear resistivity, the related $\omega$-linear behavior should show up around the gapless nodal region, while the peak of the normal self-energy vanishes at the nodal point (see Fig. 22). Therefore, the $\omega$-linear coefficient observed as the tail of the peak is not necessarily the same as the $T$-linear coefficient in the resistivity. In fact the high-energy $\omega$-linear component in Fig. 8 (c) and (d) (black squares) has substantially smaller slope than the present $\omega$-linear component directly associated with the prominent peak. This smaller slope at the high-energy region ($\omega<-0.1$ eV region) is consistent with the high-temperature $T$-linear resistivity [34] through the correspondence relation $\hbar\omega\leftrightarrow\pi k_{\rm B}T$ and the $\omega$-linear self-energy in the high-energy part identified in an earlier study [13]. We need further studies on the relation between these two $\omega$-linear components.

The $T$-linear scaling may remind the readers of quantum critical behaviors. In the present study, we examine only two sets of data for different compounds with different doping concentrations. Therefore, we could not exclude the possibility that these samples are by chance both close to quantum critical points. In fact, the underdoped sample shows more linear behavior than the optimum doped sample, which might imply that the underdoped sample is closer to the quantum critical point. Alternatively, the $T$-linear scaling behavior could emerge in a distinct phase covering a finite range of doping concentrations as discussed above as the Planckian fluid. However, it is impossible to draw a conclusion from these two samples only and it is left for future studies.

VI.2 Factors that determine the superconducting critical temperature

Fundamental quantities revealed by the machine learning provide further insight into the superconductivity through the scaling among experimental observables: The linear relation $F\propto T_{\rm c}$ between $T_{\rm c}$ and the superfluid density $F$ measured from the muon-spin relaxation rate $R$

(theoretically proportional to $F(k)z_{\rm qp}(k)$ averaged over Fermi surface momentum) has been examined through the Uemura plot [40] in high-$T_{\rm c}$ superconductors as in an example of the purple triangles in the inset of Fig. 10(b) for Bi2201 [37]. The linearity should be satisfied for attractive interaction stronger than the effective Fermi energy scale $E_{\rm F}$, which is proportional to the carrier density in two spatial dimensions. Here, $E_{\rm F}$ is roughly the effective bandwidth of the dispersion $z_{\rm qp}\epsilon_{k}$. This proposal interprets the linearity as a signature of the Bose-Einstein condensation (BEC) regime.

Homes et al. [39] proposed empirical but more universal fitting as plotted in an example of Bi2201 by blue upside-down triangles in the inset of Fig. 10(b) [37, 38], where the dc conductivity $\sigma_{\rm dc}$ at $T_{\rm c}$ enters as $R\sim CT_{\rm c}\sigma_{\rm dc}(T_{\rm c})$ with a material independent constant $C$. Since $\sigma_{\rm dc}$ is proportional to the momentum relaxation time, the Homes relation proposes qualitatively different physics involving dissipation and scattering effects beyond the naive BEC regime. However, since $\sigma_{\rm dc}$ is believed to be proportional to both the carrier density and the relaxation time, it is not easy to single out the relaxation effect. Related scaling of the superfluid density $Fz_{\rm qp}$ proportional to the quasiparticle peak weight was also proposed [41, 42].

Here, we heuristically propose a better scaling for $T_{\rm c}$ by utilizing the present finding to show the power of the machine learning and the significance of the peak. The amplitude of the self-energy peak discovered here responsible for the superconductivity has to represent the scale of the effective attractive interaction for the Cooper pair in analogy to the Eliashberg formalism and should enter the $T_{\rm c}$ scaling.

Our proposal for $T_{\rm c}$ is given by

[TABLE]

where $\Gamma$ is the damping introduced before and $F$ is the superfluid density. The factor $g$ is the scale of the effective attractive interaction as will be discussed later. It is reasonable that $T_{\rm c}$ is scaled by the mean field acting on the formation of the Cooper pair given by the product of the attractive interaction $g$ and the order parameter given by $F$ (superfluid density).

Let us first discuss how the characteristic effective attractive interaction $g$ is extracted from our self-energy analysis. We first introduce the bare attraction $\Omega_{0}(k_{\rm AN})$, which is represented by the ratio of the peak intensity $\overline{W_{\rm PEAK}}(k_{\rm AN})$ (peak intensity of $W$) in Fig. 8 to the absolute value of the peak energy $\omega_{\rm PEAK}(k_{\rm AN})$. Note that $\overline{W_{\rm PEAK}}(k_{\rm AN})/\omega_{\rm PEAK}(k_{\rm AN})$ is proportional to the gap through Kramers-Kronig relation. (Through Eq. (8), the residues of the poles of $W$, $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ follow the same scaling at the self-energy peak). Here we assume that the antinodal value is the representative of this estimate, because the gap is the maximum. Now $\overline{W_{\rm PEAK}}(k_{\rm AN})=\int d\omega{\rm Im}W_{\rm PEAK}(k,\omega)$ is the integrated intensity of the peak ${\rm Im}\Sigma_{\rm PEAK}=-{\rm Im}W_{\rm PEAK}$ plotted as the yellow and pink areas in Figs. 8(c) and (d).

We note that the quasiparticle weight constrains the coherent pairing so that the peak intensity $\overline{W_{\rm PEAK}}(k_{\rm AN})$ should involve the averaged quasiparticle weight at the peak energy. For the quasiparticle weight, instead of $z_{\rm qp}$ defined in the $\omega\rightarrow 0$ limit, we employ the renormalization factor $Q$ at $k_{\rm AN}$ (see Eq. (10)) averaged in the self-energy peak region, namely, $\overline{Q}=\int d\omega{\rm Im}W_{\rm PEAK}(k,\omega)Q(k,\omega)/\overline{W_{\rm PEAK}}(k)$ with the integration over the interval $\omega<0$. Then

[TABLE]

gives the scale of the attractive interaction.

The damping $\Gamma$ in Eq. (43) looks counterintuitive because it tells that the strongly damped electrons would have higher $T_{\rm c}$. However, in this strong coupling superconductor, the strong damping is originated from the quantum entanglement as is discussed in the previous subsection, which may promote the quantum mechanical singlet pairing. In Fig. 10(b), the fit by Eq. (43) is shown (the main panel of Fig. 10(b)).

Although the Homes plot does not offer how $T_{\rm c}$ is determined because $\sigma(T_{\rm c})\propto 1/T_{\rm c}$ cancels in the relation to $F$, the present result indeed shows $T_{\rm c}$ linearly scaled by $\Gamma(k_{\rm N})g(k_{\rm AN})F(k_{\rm AN})$. The linearity is crucially different from the Uemura plot as well because of the dependence on the relaxation rate $\Gamma$. Intuitively, $\tau=\hbar/(\Gamma k_{\rm B}T)$ or $\hbar/(\Gamma k_{\rm B}\omega)$ is related to the characteristic length scale $\lambda$ for the extension of the quantum mechanically entangled area through $\lambda\sim v_{\rm F}\tau$ [35], where $v_{\rm F}$ is the characteristic electron velocity (“Fermi velocity”). The larger attraction generates stronger self-energy peak. It necessarily generates the steeper $\omega$-linear tail of ${\rm Im}\Sigma^{\rm nor}$ near zero energy, further enhancing more local and stronger pairing adiabatically continued to the BEC limit beyond the realistic cuprate regime, and raises $T_{\rm c}$ through Eq. (43). The strange (dissipative) metal and high $T_{\rm c}$ with the strong attraction represent the two sides of the same coin. It is also interesting to note that Eq. (43) looks compatible with the scaling $E_{\rm c}\propto\gamma_{\rm c}T_{\rm c}^{2}$,where $E_{\rm c}$ is the condensation energy and $\gamma_{\rm c}$ is the Sommerfeld constant of the specific heat [43, 44], because $gF$ plays the role of the gap, which generates the energy gain.

Note that $\Gamma(k)$ should be analyzed around $k_{\rm N}$, while $\Omega_{0}$ and $F$ contribute to $T_{\rm c}$ at $k_{\rm AN}$ for better fitting. The in-plane transport and the quantum entanglement are dominated by the contribution around the nodal region, while the pairing looks driven in the antinodal region, both of which contribute to raise $T_{\rm c}$.

Equation (43) is the best scaling among various attempts we have made. To convince readers, we just show two examples of plot in Figs. 10(c) and (d). The first example is $z_{\rm qp}(k_{\rm AN})F(k_{\rm AN})$ vs. $T_{\rm c}$ plotted in Fig. 10(c), which, though not perfectly equivalent, apparently mimics the Uemura plot. The second is $z_{\rm qp}(k_{\rm N})F(k_{\rm AN})$ vs. $1/z_{\rm qp}(k_{\rm N})c_{1}(k_{\rm N})$ which mimics the Homes plot, because $T_{\rm c}\sigma(T_{\rm c})\propto 1/z_{\rm qp}(k_{\rm N})c_{1}(k_{\rm N})$ is expected (see Fig. 10(d)). The standard deviation is by far best for the present fit in Fig. 10(b) with Eq. (43).

The primary origin of the large difference of $T_{\rm c}$ between Bi2201 and Bi2212 is identified as the difference in $\overline{W_{\rm PEAK}}$ (namely the coupling strength of the electron with the dark object which makes the prominent self-energy peaks), supplemented by the difference in $\omega_{\rm PEAK}$. (see Table 1. $\overline{W_{\rm PEAK}}(k_{\rm AN})\sim 7.6\times 10^{-3}$ eV2 and $\omega_{\rm PEAK}(k_{\rm AN})\sim 0.07$ eV for the optimum Bi2201 and $\overline{W_{\rm PEAK}}(k_{\rm AN})\sim 1.4\times 10^{-2}$ eV2 and $\omega_{\rm PEAK}(k_{\rm AN})\sim 0.045$ eV for the optimal Bi2212 at $k=k_{\rm AN}$). (More quantitative details, see also the list of $\Omega_{0}(k_{\rm AN})$ in Table 1 and the angle and doping dependences of $\overline{W_{\rm PEAK}}(k)/\omega_{\rm PEAK}(k)$, and $\overline{Q}(k)$ in Fig. 25.)

A recent intensive study [45] on clean thin films of (La,Sr)2CuO4 has revealed that the Uemura plot shows a linear scaling for $T_{\rm c}\gtrsim 10$ K and a quadratic scaling, where $T_{\rm c}$ is proportional to the square root of the relaxation rate, for $T_{\rm c}\lesssim 10$ K. Consequently, the higher $T_{\rm c}$ linear scaling extrapolates $T_{\rm c}$ to a finite value at the zero relaxation rate limit. The present scaling modifies the higher $T_{\rm c}$ linear scaling by taking into account $g\Gamma$. If $g$ is approximately a constant for a given crystal structure and, $\Gamma$ increases by decreasing the doping, as already reported in the overdoped region [46], the linear scaling in Fig. 10(b) is consistent with the sub-linear scaling of $T_{\rm c}$ to $F$ observed in Ref. 45. We note that, if $\Gamma$ diverges as $T_{\rm c}$ decreases in the underdoped region, intrinsic doping dependence of $\Gamma$ may explain the square-root scaling at the underdoped limit qualitatively. However, the quantitative estimate of the power of the scaling is beyond our scope. On the other hand, ARPES measurements on clean samples at the low $T_{\rm c}$ limit in the overdoped region are not available.

The present result is significant because the whole analyses are obtained solely from the single ARPES line shape of $A(k,\omega)$ and contains much less ambiguity than before. The present machine learning purely from experimental data sheds new light on understanding the superconducting mechanism, where the energy dissipation plays a role through the extension of quantum entanglement. For detailed doping concentration and momentum dependences of $T_{\rm c}$, $F(k)$, $c_{1}(k)$, $z_{\rm qp}(k)$, and the superconducting gap $\Delta_{0}(k)$ for Bi2201 at 11K are found in Appendix O.

VI.3 On cancellation of peaks

On the cancellation of the two self-energy contributions which make the superconducting temperature high and the role of dissipation in determining $T_{c}$, it is desired to examine the present results in other cuprate compounds by measuring $A(k,\omega)$ at high accuracy and resolution.

The present finding also call for research to find the microscopic origin of the canceling self-energy peaks. A two-component model was proposed [47, 48] for the cancellation of the self-energy poles, but even if it is the case, based on the present experimental evidence, further pursuit for the microscopic description of the hidden-fermion excitation and its interaction is highly desirable. In the two-component fermion theory, electrons are fractionalized into the bare electrons and dark fermions (hidden fermions) consistently with the cluster dynamical mean-field theory (cDMFT) [47, 48]. However, the prominent tail of the peak extending to $\omega\sim 0$ as $\omega$-linear feature suggests that the dark fermion must have strong interaction effect. See also Ref. 49 for the prediction of the fractionalization expected in other spectroscopic data such as the resonant inelastic X-ray scattering.

If the spontaneous symmetry breaking such as the stripe order coexists with the superconductivity, the cancellation may be accounted for as an alternative interpretation [47, 48]. A phenomenological resonating valence bond theory [50] also accounts for the cancellation because of the same fractionalized description as is discussed in Ref. 51. The present result poses severe constraints on possible theories. Whether there exist other origins of the peaks and their cancellation rather than the above possibilities would be equally intriguing based on the present finding.

VI.4 Kink

A kink structure was observed in the dispersion of momentum distribution curve (MDC) peak of Bi2212 and other compounds mainly near the nodal point [52, 53, 54]. In addition, a kink-like structure was identified in the study on the Hubbard model [55]. It is an intriguing future issue to study whether this kink in the MDC peak dispersion has any connection to the peak and accompanied sudden change of the slope and sign in the real part of normal self-energy found here in Figs. 8(a) and (b), because the energy scale $\sim-0.07$ eV of the sign change for Bi2212 is similar to the kink energy scale.

VII Summary and Outlook

We have formulated a reliable way of extracting the normal and anomalous self-energies separately from the ARPES data of superconductors by taking advantage of recently developed machine learning method. Careful benchmark tests including simple metals, conventional BCS superconductors and model systems that have established solutions indicate that the method offers an accurate and reliable regression of the self-energies only from the ARPES data even for challenging strongly correlated electron systems.

Then the method has been applied to cuprate superconductors, Bi2201 and Bi2212. We have successfully extracted the normal and anomalous components of the self-energy from the ARPES spectra.

In contrast to previous studies, the result shows that the imaginary part of the normal and anomalous self-energies have prominent peak structures as a function of the frequency. However, the contributions of the normal and anomalous components cancel in the spectral function, which accounts for the failure to identify the prominent structure for long time. Nevertheless the peak in the anomalous self-energy has been shown to generate more than 90% of the superconducting gap and thus turned out to be the primary source of the superconductivity. Therefore, the discovered these peak structures and their cancellation pose a severe constraint with insight on the mechanism of the high transition temperature of superconductivity.

The origin of the failure in identifying the peak structures in the previous studies is elucidated in Appendix A to be primarily the assumptions in the previous studies: The previous studies assumed linear energy dispersion, and/or self-energies that are momentum independent along the direction perpendicular to the Fermi surface, for example in Ref. 13. These assumptions are not justified in the strongly correlated electron systems.

Thus newly obtained quantities hidden in the direct experimental measurements in the past allow us to show that the superconducting transition temperatures are well scaled by the product of the superfluid density $F$, the effective attractive interaction $g$ and the Planckian dissipation $\Gamma$.

Present successful examples of insight obtained purely from the machine learning analysis of experimental data indicates an opening of a promising field which allows understanding physics hidden in experiments, without relying on involved and specific theoretical assumptions and constraints that are not shown to be justifiable in strongly correlated electron systems. At the same time, we have shown that very accurate experimental data are required to extract hidden quantities. For instance, in Refs. 56, 57, and 58, the signal-to-noise ratio in the ARPES measurements for Bi2212 seems to be already sufficiently small while available momenta in the Brillouin zone are limited. Although Ref. 59 for Bi2212 and Ref. 60 for Bi2201 reported the ARPES spectra in a wide range of momenta, the signal-to-noise ratio does not seem to be small enough. It is important to improve the experimental resolution and suppress errors to increase the reliability of the machine learning inference.

In addition combining with other independent measurements such as the quasiparticle interference obtained from the scanning tunneling microscope in this case is important to reach better statistics. By combining with other experimental data and indisputable theoretically basic constraints such as symmetry, much more powerful tool will be provided for understanding physics of complex phenomena.

The present study will stimulates studies on the origin of the peak structures. Indeed, there have been studies [61, 62] on the origin of the peak structures, which are inspired and motivated by our results [63].

Acknowledgements.

We thank Takeshi Kondo and Adam Kaminski for providing us ARPES data published in Refs. 12 and 28. We also thank Takeshi Kondo for discussions on the experimental results. We are grateful to Shiro Sakai for discussions and comments on the manuscript and Chandra Varma for clarification of the procedure of the analysis in Ref. 13. This research was supportd by MEXT as “Priority Issue on Post-K computer” (Creation of New Functional Devices and High-Performance Materials to Support Next-Generation Industries (CDMSI)) and “Basic Science for Emergence and Functionality in Quantum Matter

• Innovative Strongly-Correlated Electron Science by Integration of Fugaku and Frontier Experiments -” (JPMXP1020200104) as a program for promoting researches on the supercomputer Fugaku, supported by RIKEN-Center for Computational Science (R-CCS) through HPCI System Research Project (Project ID: hp170263, hp180170, hp190145, hp200132 and hp210163). Y. Y. was supported by PRESTO, JST (JPMJPR15NF). Y. Y. and M.I. were supported by JSPS KAKENHI (Grant No. 16H06345). A.F. was supported by KAKENHI (Grant No. 19K03741). The present regression is performed by using our house code. The numerical code will be available upon request.

Appendix A Comparison with Previous Studies on Self-Energy

Here, we discuss the comparison with the analysis by Bok et al. [13], which has not clearly identified a prominent peak structure in the normal and anomalous self-energies. We first point out that the primary origin of the discrepancy may be the momentum region they studied. They have analyzed mainly only around the nodal region and at most up to $\theta=20^{\circ}$ measured from the nodal point, which is far away from the antinodal point. This makes the identification of the prominent peak difficult. In Appendix L, we show the momentum dependence of the EDC curve for the optimally doped Bi2201 (see Fig. 22). In this notation, $20^{\circ}$ from the nodal point in [13] nearly corresponds to the point between Fig. 22(h). It is natural that the peak structure is not clearly visible there. However, if we take a close look, aside from the clear difference of the featureless slope arising from the instrumental difference and presumable different background effect, the peak-like (or shoulder-like) structures at -0.06 eV in the imaginary part of the normal and anomalous self-energies in Fig. 3 of Ref. 13 shares a common feature with our result in Fig. 22(h). Unfortunately, this tiny signature and the cancellation in the spectral function is hardly identified conclusively because of the momentum far from the antinodal point.

More importantly, a crucial origin of the underestimate of the peak is the usage of the Dynes function as explained below. It leads to the underestimation of peak-like structures in the imaginary part of the self-energies. The underestimation inevitably leads to the difference in the self-energies obtained in the present paper and Ref. 13.

Before going into the explanation of the underestimation, we review the method used in Ref. 13 to make the discussion self-contained. The method assumes that the real part of the Dynes function $N(\omega)$ defined as,

[TABLE]

is equal to the ratio of the integrated spectral function at the normal and superconducting states as

[TABLE]

where the integrated spectral functions are obtained with respect to the bare dispersion $\epsilon(k_{\perp})$ along momentum perpendicular to the Fermi surface as $\mathcal{A}_{S}(\omega)=\int dk_{\perp}A_{s}(k_{\perp},\phi,\omega)$ and $\mathcal{A}_{S}(\omega)=\int dk_{\perp}A_{n}(k_{\perp},\phi,\omega)$. Here, $A_{s}$ ($A_{n}$) is the spectral function of the superconducting (normal) state. The relation, however, holds only when the bare dispersion $\epsilon(k_{\perp})$ is linear at the large bandwidth limit. Practically, the method approximately works around the nodal direction where the bare band is linear within a certain energy range. However, even around $\phi=20^{\circ}$, the approximation does not work, where $\epsilon(k_{\perp})$ shows strong deviation from the linear dispersion, as shown below.

Here, we concretely demonstrate that the relationship ${\rm Re}N(\omega)=\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)$ indeed underestimates the imaginary part of the self-energies. As the most striking example, we would like to demonstrate that, even if there are peak structures in the genuine self-energies, the self-energies estimated by using ${\rm Re}N(\omega)=\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)$ collapse to shoulder structures instead of the peak structures.

When there is no background $b(\omega)$ in the experimentally observed spectral function, the scheme to extract the self-energies used by Bok et al. is summarized as follows:

(1)

${\rm Re}N(\omega)$ is (initially) estimated by $\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)$.

(2)

${\rm Im}N(\omega)$ is obtained from ${\rm Re}N(\omega)$ by using the KK transformation. Then, $\Delta(\omega)$ is obtained from ${\rm Re}N(\omega)+i{\rm Im}N(\omega)=\{\omega/\sqrt{\omega^{2}-\Delta^{2}}\}$.

(3)

By inputing $\Delta(\omega)$ into the Green function in the superconducting phase, $Z(\omega)=1-\Sigma^{\rm nor}(\omega)/\omega$ and $\Sigma^{\rm ano}(\omega)=\phi(\omega)=\Delta(\omega)\Sigma^{\rm nor}(\omega)$ are obtained.

Below, we estimate the self-energies by following Bok et al. [13].

As a model self-energies that show peak structures, we take the self-energies at $\phi=21.24^{\circ}$ obtained in 22(h), which is shown in the left panel of 11. We also assume a typical bare dispersion for bismuth cuprates, $\epsilon(k)=\mu-2t_{1}(\cos k_{x}+\cos k_{y})+4t_{2}\cos k_{x}\cos k_{y}-2t_{3}(\cos(2k_{x})+\cos(2k_{y}))$, where $\mu=405$ meV, $t_{1}=360$ meV, $t_{2}=108$ meV, and $t_{3}=36$ meV. Then, we integrated the spectra $\mathcal{A}_{S}(\omega)$ with $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ along the momentum cuts (in the first quadrant of the Brillouin zone) specified by the angle $\phi$ measured from the antinode and obtain $\mathcal{A}_{N}(\omega)$. To estimate $\mathcal{A}_{N}(\omega)$, we generate the normal state self-energy $\Sigma^{\rm nor}$ by eliminating peak structures that cancel the peak structure in $\Sigma^{\rm ano}$, which is shown in the middle panel of Fig. 11 by assuming that the effect of the superconductivity appears at the peak structure only in the normal self-energy. Then, in the right panel of Fig. 11, we obtain the estimated $\Sigma^{\rm nor}$ (black curve) by assuming ${\rm Re}N(\omega)=\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)$, in comparison with the original $\Sigma^{\rm nor}$ taken from Fig. 22(h). Although we started from $\Sigma^{\rm nor}$ and $\Sigma^{\rm ano}$ that contain prominent peaks, in the resultant $\Sigma^{\rm nor}$ derived by assuming the above Dynes function shows the disappearance of the peak. This shows that the usage of the assumption of the Dynes function leads to the self-contradiction with the underestimate of the peak structure.

As shown in the right panel of Fig. 11, the estimated ${\rm Im}\Sigma^{\rm nor}$ shows the shoulder structure around -50 meV, instead of the peak around -70 meV in the original ${\rm Im}\Sigma^{\rm nor}$. The shoulder structure of ${\rm Im}\Sigma^{\rm nor}$ found in Fig.3D of Ref. 13 around -50 meV is interpreted as a remnant of the original peak structure in the genuine self-energies, which is consistent with our present results. Therefore, the difference between our results and those of Ref. 13 is attributed to the artificial reduction in the amplitude of ${\rm Im}\Sigma^{\rm nor}$ due to the assumption ${\rm Re}N(\omega)=\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)$ in Ref. 13.

The self-energies obtained in Ref. 13 also show a remnant of the cancellations of the peak structures of ${\rm Im}\Sigma^{\rm nor}$ and ${\rm Im}W$ found in the present paper. As shown in 12, the shoulder structure of ${\rm Im}\Sigma^{\rm nor}$ obtained in Ref. 13 is canceled by ${\rm Im}W$. Then, ${\rm Im}\Sigma^{\rm tot}={\rm Im}\Sigma^{\rm nor}+{\rm Im}W$ does not show any shoulder structure.

The above general underestimate (or elimination) of the self-energy peak is a universal failure of the assumption employed by the Dynes function and the assumption of the wide band limit. The physical origin is the following: In the wide band width limit, $\mathcal{A}_{N}(\omega)$ only contains the information of the bare band, which does not depend on whether the system is normal or superconducting. Then, the left-hand side of $\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)={\rm Re}\{\omega/\sqrt{\omega^{2}-\Delta^{2}}\}$ does not contain not only the normal-state self-energy but also the bare band information, which are absent in the right-hand side of the relation. However, even at $\phi=20^{\circ}$, the assumption of the wide band width is invalid. Therefore, the normal-state self-energy information remains in the denominator of the ratio $\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)$,

[TABLE]

where $Z_{\rm n}(\omega)$ ($Z_{\rm s}(\omega)$) is $Z(\omega)=1-\Sigma^{\rm nor}(\omega)/\omega$ in the normal (superconducting) state. To compensate the contribution of $\Sigma^{\rm nor}$ in the left-hand side of the relation $\mathcal{A}_{S}(\omega)/\mathcal{A}_{N}(\omega)={\rm Re}\{\omega/\sqrt{\omega^{2}-\Delta^{2}}\}$, the normal-state self-energy $\Sigma^{\rm nor}$ below $T_{\rm c}$, whose contribution is contained in $\mathcal{A}_{S}(\omega)$, may resemble the self-energy $\Sigma^{\rm nor}$ above $T_{\rm c}$, whose contribution is contained in $\mathcal{A}_{N}(\omega)$. When the normal state self-energy does not show significant structures and, thus, does not introduce pseudogap and anomalies, such as peak-dip-hump structures, in the normal-state spectral functions, the estimated superconducting self-energies indeed do not show significant structures. We note that the shoulder structure consistently present in Ref. 13 and 12b in ${\rm Im}\Sigma^{\rm nor}$ without a signature of the pseudogap around 20 degrees momentum could be contributed from the background. However, the above analysis does not change even when the background is subtracted, because it appears both in $\mathcal{A}_{S}(\omega)$ and $\mathcal{A}_{N}(\omega)$ and the effect of the subtraction cancels.

In addition, Bok et al. [13] have assumed that the normal and anomalous self-energies are momentum independent along the direction perpendicular to the Fermi surface. The assumption imposes crucially restrictive condition when one infers the self-energy and superconducting gap function from the momentum distribution curve (MDC) as in Ref. 13 (see Eq. (S8) in the Supplementary information of their paper). This assumption is adequate in the BCS superconductors in conventional weakly correlated systems, while it is questionable in the present strongly correlated cases such as the cuprates, where even the normal self-energy can be singular and strongly dependent on the momentum. In addition to the restricted momentum dependence assumed in [13], they have assumed the quasiparticle representation along the momentum perpendicular to the Fermi surface (or the Lorentzian form of the MDC) and the spectral function that cannot be represented by the Lorentzian form are interpreted as the background, while these have not been assumed in the present paper, because, though it is satisfied in the Fermi-liquid normal state, it is not clear whether these constraints are satisfied in the strong coupling cuprate superconductors. Indeed, there exist a number of numerical evidences for the violation of the assumption: For instance, the singular momentum dependence of the normal self-energy with emergence of the coexisting zeros and poles of the Green function [64] and non-quasiparticle spectral function [65, 66] in doped Mott insulators. The machine learning is more fit in solving the present problem, because the flexible fitting of the self-energy function is required at least away from $\omega=0$ particularly in the antinodal region, where the breakdown of the quasiparticle picture is apparent. Concerning the difference in the high-energy part ($\omega<-0.1$ eV) of the original two ARPES data (namely in [13] and [12]), we have already analyzed the effect of the possible extrinsic high-energy part in F and has shown that it does not affect the structure of the peaks as clarified in 19 and it cannot be the origin of the difference in the peak structure.

Next we discuss the origin of discrepancy in the result by Li, et al. [31]. They assumed momentum independent self-energy in the MDC analysis as one sees in Eq. (3) of their supplementary note 2, whose basis is unclear. More crucially, they assumed the imaginary part of self-energy in the forms Eq. (5) or (6) in Supplementary Information of Ref. 31 for the normal part and

[TABLE]

with constant fitting parameters $W_{3},W_{4},E_{3}$ and $E_{4}$, for the anomalous part, which is unjustified. Particularly, the assumed form Eq. (45) does not allow the formation of peak or dip and it does not allow the cancellation with the structure in the normal contribution in the spectral function as we discovered. It is crucially important to allow the flexibility of the self-energy form and the machine learning is one of the best way to incorporate it while the attempt by Li, et al. failed in implementing the flexibility. We have attempted to fit the self-energy with the constraint of Eq. (45) for the anomalous part and found the resultant optimized $\chi^{2}$ is {\color[rgb]{0,0,0}\chi_{\rm MLF}^{2}}=6.1\times 10^{-6}, which is much higher than the present result {\color[rgb]{0,0,0}\chi_{\rm ML}^{2}}=2.1\times 10^{-6}. Because the experimental resolution is $\chi_{\rm exp}^{2}=1.4\times 10^{-6}$ as mentioned above, the intrinsic $\chi^{2}$ defined by $\delta\chi_{\rm MLF}^{2}=\chi_{\rm MLF}^{2}-\chi_{\rm exp}^{2}$ is $4.7\times 10^{-6}$. Then the intrinsic machine learning error in the unit of the experimental standard deviation is {\color[rgb]{0,0,0}\delta\chi_{\rm MLF}/\chi_{\rm exp}=1.8}, which is close to the twice of the standard deviation. Namely, the probability that this constrained choice is true is less than 7%. The resultant self-energy does not show any appreciable peak as it should be in contrast to the present result as one sees in Fig. 13. Note that, despite the constrained anomalous part, due to the unbiased choice of the normal self-energy here, the self-energies in Fig. 13 should be much better fit than the more constrained ones (including the normal part) in Li, et al [31]. The spectral function obtained from the additional constraint in the normal part, Eq. (5) or (6) in Supplementary Information of Ref. 31 must give even higher $\chi^{2}$ than $\chi_{\rm MLF}^{2}$.

In addition, the form of the imaginary part of the anomalous self-energy is unphysical. The form assumes the attractive interaction ranging to the infinite frequency scale (namely, instantaneous attractive interaction), which can not exist in the real experiments.

In summary, constraints or assumptions unjustified a priori such as momentum independent self-energy, strictly constrained but unjustifiable self-energy form, and Eliashberg equation in the previous works do not allow prominent peak structures and at the same time give higher errors in the regression of the experimental data than our result indicating the superiority of our analysis. It supports that the cancellation of the normal and anomalous contribution in the total self-energy found in our analysis and not in the former studies should be considered seriously.

Appendix B Pair Breaking Effect

The renormalization factor (quasiparticle residue) estimated from the expression

[TABLE]

theoretically equivalent to $z_{\rm qp}$ defined in Eq. (12) is the weight of the quasiparticle, which can be substantially reduced from the noninteracting value $z(k)=1$ due to the interaction effects. The renormalization factor estimated from the fitting of Eq. (46) is $z_{\rm qp}\sim 0.1$ for Bi2212 and $z_{\rm qp}\sim 0.03$ for Bi2201 (see Fig. 14) supporting the non-Fermi liquid behavior especially in the underdoped case.

As shown in Figs.  14(c) and (d), the non-Fermi-liquid-like ${\rm Im}\Sigma^{\rm nor}(\omega)$ affects the gap function $\Delta(\omega)$ (in Eq. (9) through $Q(\omega)$ (see Eq. 10). In general, negative ${\rm Im}\Delta(\omega)$ for $\omega<0$ enhances ${\rm Re}\Delta(\omega=0)$ through the Kramers-Kronig relation and is indeed negative in most of $\omega$ in Fig. 14(c) and (d). However, ${\rm Im}\Delta(\omega)$ is positive at $|\omega|<0.04$ eV ($|\omega|<0.06$ eV) for Bi2212 (Bi2201). Because ${\rm Im}\Sigma^{\rm ano}$ is found to be always negative for $\omega<0$, it is ascribed to the pair breaking effect of $Q$, arising from poles of $\Sigma^{\rm nor}$ inside the superconducting gap as already pointed out [47]. The pair breaking is much more prominent for underdoped sample, Bi2201.

Although a similar conclusion for the underdoped Bi2201 suggests a universal nature, the prominent non-Fermi liquid behavior and the pair breaking could be accounted for by an alternative at $k=k_{\rm AN}$, namely the pole of $\Sigma^{\rm nor}$ shifts to the energy $\omega\sim 0$ and destroys $\Sigma^{\rm ano}$ accompanied by an insulating gap. Although such a solution gives worse $\chi^{2}$ in our analysis, a momentum selective insulating behavior at the antinodal point deserves to be explored further together with the full momentum and temperature dependences.

Appendix C Gap Functions

C.1 Resolution of gap functions

The gap function $\Delta(\omega)$ defined in Eq. (9) can show significant $\delta$ dependence near the small $\delta$ limit around $\omega\sim 0$. The $\delta$ dependence originates from the finite imaginary part of the normal self-energy ${\rm Im}\Sigma^{\rm nor}(k,\omega=0)$ inevitable in the experimental data. When we modify $Q$ as

[TABLE]

we obtain stable behaviors of $\Delta(\omega)$ for $|\omega|>10$ meV by keeping $\delta=10$ meV and restricting to $\delta^{\prime}<\delta$. In Fig. 14, we use $\delta^{\prime}=2.5$ meV.

C.2 Peak shift in gap function

Fig. 14(c) and (d) show the gap function $\Delta$ defined in Eqs. (9) and (10). It reveals that the peak positions in $\Delta(k_{\rm AN},\omega)$ are different from those in ${\rm Im}\Sigma^{\rm ano}$, which is consistent with the hidden fermion theory [47]. In fact, the peak positions of ${\rm Im}\Delta(k_{\rm AN},\omega)$ ($\sim\pm 80$ and $\pm$220 meV for Bi2212 and $\sim\pm 80$ and $\pm$210 meV for Bi2201) are nearly the same as the peak positions of ${\rm Re}\Sigma^{\rm ano}$, while the peak positions of ${\rm Re}\Delta(k_{\rm AN},\omega)$ ($\sim\pm 50$ and $\pm 180$ meV for Bi2212 and $\sim\pm 50$ and $\pm 160$ meV for Bi2201) are nearly the same as the peak positions of ${\rm Im}\Sigma^{\rm ano}$. This is because the imaginary part of $Q$ is dominant in the relevant frequency region ($\sim 100-200$ meV) as shown in Fig. 14. The shift of the peak positions indicates the strong renormalization effect in the normal quasiparticle contained in $Q$. In any case, in the contribution to the real order parameter of the superconductivity $\Delta_{\rm AN}={\rm Re}\Delta(k_{\rm AN},\omega=0)$ is expected to be contributed mostly from the two peaks in ${\rm Im}\Delta(k_{\rm AN},\omega)$ through the Kramers-Kronig relation. The $d$-wave gap amplitude $\Delta_{\rm AN}$ is 30 meV for the optimally doped Bi2212 while it is around 10 meV. However, at small energy ($\sim 60$ meV), the gap amplitude is both around 40 meV, which is comparable to the peak energy of ${\rm Im}\Sigma^{\rm ano}$ and ${\rm Im}\Sigma^{\rm nor}$ suggesting the similar pseudogap energy for optimum and underdoped samples.

Appendix D

Accuracy, Stability and Robustness Tested by Solvable Benchmarks

In this section, we employ exactly solvable models as benchmarks. We examine whether our machine learning correctly reproduces the exact self-energies (with prominent peak structures), if the exact solution indeed shows the cancellation of the normal and anomalous self-energy contributions in the total self-energy and the spectral function $A(\omega)$ shows only a weak peak-dip-hump structure. To our knowledge, exact solution, which shows such a cancellation is not found except for the case of the two-component fermion model. Then, as a benchmark, we inferred the self-energy of a superconducting two-component fermion model defined by the following Lagrangian,

[TABLE]

which is essentially the same form as that introduced in Ref. 47 and discussed in Ref. 48. In the following discussion, we assume that the noninteracting density of states determined from $\epsilon_{c}(k)$ is a constant $N_{F}$ and focus on a specific momentum $k$ at the Fermi momentum just for simplicity. Because of the momentum independence, this consideration at a specified momentum does not cause loss of generality. Here, we add $\Sigma^{(0)}(\omega)$ at the above momentum defined by

[TABLE]

in addition to $\epsilon_{c}(k)$ to mimic the additional normal Fermi-liquid-like component seen in the experimental result arising from interaction effect for the part not represented by the coupling to the $d$ fermion, where $a$, $b$ and $\alpha$ are constants. The self-energy in the exact solution of this two-fermion model is given as,

[TABLE]

For simplicity, we have dropped the momentum dependence in the solutions (LABEL:Snor) and (51). In our calculation, we set $a=0.008$ eV2, $b=0.2$ eV2, $\alpha=0.08$ eV3, $V_{1}=0.075$ eV, $D_{1}=0.0375$ eV, and $\epsilon_{d}=\epsilon_{c}=0$. The present choice of the parameters is enough to generate the spectral function observed at the Fermi momentum we focus, and the dependence on the doping and dimension of the system etc. are implicitly contained in $N_{F}$. By using the exact solution for the spectral function $A_{2f}(\omega)$, we add small but finite noise, where $\sigma^{2}$ of the noise is set to be $6\times 10^{-4}$ to simulate the role of noise in experiments and perform the machine learning using this noisy $A_{2f}(\omega)$. Although it is irrelevant to the inferred self-energies, the Fermi-Dirac distribution with $T=40$ K is introduced in the spectral function used in the machine learning just by following the scheme with finite-temperature experimental data. In Figs. 15(a)-(c), the spectral function of the two-component fermion model $A_{2f}(\omega)$ and the self-energies are shown for the exact solutions (solid and dashed curves) and the machine learning results (symbols). In the self-energy inference, we choose $n_{0}=0.45$ without fine tuning and use the spectrum within $-0.55\ {\rm eV}\ <\omega<0.05\ {\rm eV}$. The peak position in ${\rm Im}\Sigma^{\rm nor}$ and ${\rm Im}\Sigma^{\rm ano}$ and the Fermi liquid-like normal contribution in the exact solution (shown with the index $2f$ such as $A_{2f}(\omega)$ illustrated by solid and dashed curves) are well reproduced by our machine learning results (symbols). The peak cancellation in $\Sigma_{2f}^{\rm tot}$ in the exact solution is also well reproduced.

In Figs. 15(d)-(f), we show an artificial case, where the pole of $\Sigma^{\rm ano}$ is shifted 0.03 eV from the solution (51), where the pole cancellation in the total self-energy does not occur any more and the spectral function shows weird two peaks. The concrete representation of the self-energies of the modified two-component model is given by the following form,

[TABLE]

where the pole of $\Sigma^{\rm ano}$ is shifted from that of $\Sigma^{\rm ano}$ because of $\Delta D_{1}=0.03$ eV and a factor $r=2$ is introduced to avoid singular spectrum. The spectrum and self-energies of the modified two-component model are denoted with the index m as $A_{\rm m}(\omega)$. Even in this case, the machine learning results well reproduce all the line shapes. This indicates that our machine learning flexibly and accurately reproduces the exact solution irrespective of the presence or absence of the peak cancellation.

Appendix E Numerical Optimization Procedure

E.1 Details in minimization of training error

The parameters in the Boltzmann machine, \mbox{\boldmath\alpha}^{\rm nor}=(b,\{W_{\ell m}\}) and \mbox{\boldmath\alpha}^{\rm ano}=(w_{\lambda},\{b_{\ell}^{\lambda}\},\{V^{\lambda}_{\ell m}\}), are optimized by using the standard gradient method. The parameters at the $k+1$th step, \mbox{\boldmath\alpha}^{\rm nor}_{k+1} and \mbox{\boldmath\alpha}^{\rm ano}_{k+1}, are updated from the $k$-th values as

[TABLE]

where

[TABLE]

and $\|\cdots\|_{1}$ represents $L_{1}$ norm. The factors \left(\|S^{-1}\mbox{\boldmathg}^{\rm nor}_{k}\|_{1}\right)^{-1/2} and \left(\|\mbox{\boldmathg}^{\rm ano}_{k}\|_{1}\right)^{-1/2} are introduced to accelerate the optimization. Here, we use the natural gradient method to optimize the variational parameters in ${\rm Im}\Sigma^{\rm nor}(\omega_{j})$ because of its efficiency, [20, 67, 68] while the simple steepest descent method is employed to optimize the part of ${\rm Im}\Sigma^{\rm ano}(\omega_{j})$ because the natural gradient method assumes that the optimized distribution is positive or negative definite, while ${\rm Im}\Sigma^{\rm ano}(\omega_{j})$ does not satisfy this condition. During the optimization of the Boltzmann machine, we may introduce a regularization term by $L_{1}$ norm of the mixture of the Boltzmann machines as $\lambda_{w}\sum_{\lambda}|w_{\lambda}|$. While $\lambda_{w}=10^{-3}$ will accelerate the optimization, the results of the optimization is confirmed to be insensitive if $\lambda_{w}\leq 10^{-3}$. In the actual fitting, we employed $\lambda_{w}=10^{-3}$.

E.2 Parameters in optimization

In the present paper, first, we optimize the Boltzmann machine with $L=8$ visible nodes and $2L=16$ hidden nodes for the part ${\rm Im}\Sigma^{\rm nor}$ and, then, we enhance the resolution with $L=9$ visible nodes and 18 hidden nodes to obtain better resolution with reasonable numerical cost. In the optimization with $L=9$, we skip the outer loop (the update of the center of mass by the Bayesian process) to reduce the computational cost and perform longer minimization steps up to $2\times 10^{4}$. We employ the broadening factor $\delta=10$ meV throughout this paper. We show in Appendix E.3 that the result does not sensitively depend on the choice of $\delta$.

E.3 Effects of resolution $\delta$

In the present study, the small imaginary part $i\delta$ utilized in the Green functions is chosen to be equal to the experimental resolution. When substantially larger resolution $\delta$ is taken, the detailed spectra are trivially not reproducible. On the other hand, when smaller resolution $\delta$ is taken, the spectra may be easily fitted. Here, we examine how the smaller $\delta$ affects the inferred self-energy. As a typical example, we take $\delta$=5 meV, which is a half of $\delta$ used in the main article, and confirmed that the smaller $\delta$ does not change the qualitative structure of the self-energy. As shown in Fig. 16, the peak structures in ${\rm Im}\Sigma^{\rm nor}$ and ${\rm Im}W$, and the cancellation between them are reproduced.

E.4

Gaussian distribution represented by Boltzmann machine

When we choose the parameters as

[TABLE]

in Eqs (26) and (27), the Boltzmann machine easily represents the Gaussian distribution with the center $x_{\lambda}$, variance $s_{\lambda}^{2}$, and weight $w_{0\lambda}$, which is a localized sparse distribution. Superposition of the Gaussian distribution can easily be expressed by Eq. (28) by taking $M$ larger than 1 (typically we take $M$ several).

Appendix F Robustness of Machine Learning

The present use of machine learning is categorized to a general class of regression analysis as addressed in the first paragraph of Sec. II.1. In the standard simple case of the regression task, training data set is simply given by the observed $A$ at discrete number of $x$ and we infer the functional form of $A(x)$. In the present case, it is more involved and the training data is the experimentally measured discrete and limited number of $A$ and $\omega$, and the regression task is to determine $\Sigma$ as a continuous function of $\omega$. In terms of the optimization with the machine learning, our task is to minimize the difference between the measured data $A$ and that obtained from the inferred $\Sigma(\omega)$, which is a continuous function of $\omega$. Therefore, our work is categorized to the machine learning application to a regression task, one of the most widely applied machine learning fields. Our regression scheme is illustrated in Figure 5.

In the regression analysis, it is helpful to examine the reliability of the machine learning by using solvable cases as the benchmark, as in other type of the regression task found in the problem of solving quantum many-body problems and classical statistical physics problems[10] It is also important to test the stability of the procedure by adding noises. In this section we show the robustness against the noise and in D, we show several benchmark tests for solvable models.

F.1 Stability against noise

We examine stability of the present machine-learning scheme. By using $A^{\rm fit}$ and $\sigma_{n}^{2}$ introduced in Sec. III.5, we can generate synthetic experimental data with the same or larger amplitude of noise than the original data. Here, we use the synthetic data to examine the input data dependence of the present scheme.

Here, $\chi^{2}$ of the optimized $A(k,\omega)$ by the machine learning from the synthetic ARPES spectrum $A^{\rm syn}$ generated by the maximally-likelihood inference of the ARPES spectrum is given by $\chi_{\rm ML}^{2}\equiv\overline{\chi^{2}}=\sum_{s}^{N_{\rm d}}\sum_{r}^{N_{r}}(A^{\rm ML}(\omega_{s}^{(r)})-A^{\rm syn}_{r}(\omega_{s}^{(r)}))^{2}/N_{\rm d}N_{r}=2.1\times 10^{-6}$ (defined in Eq. (34)), which is the same level as the experimental $\chi^{2}$, namely $\chi_{\rm exp}^{2}\equiv\sigma_{n}^{2}=\sum_{i}^{N_{\rm d}}(A^{\rm exp}(\omega_{i})-A^{\rm fit}(\omega_{i}))^{2}/N_{\rm d}=1.4\times 10^{-6}$ obtained in Eq. (33). The same level of $\chi^{2}$ value indicates that the machine learning optimization to fit the experimental $A(k,\omega)$ is successfully achieved within the limit of the level of the experimental noise. The standard deviation of the experimental uncertainty is around $\chi_{\rm exp}=1.2\times 10^{-3}$. To estimate the likelihood (degree of certainty) of the present solution as the experimental interpretation, we used a standard index (for noise, variance and bias decomposition, see Ref. 21) expressed as $\delta\chi_{\rm ML}=\sqrt{\chi_{\rm ML}^{2}-\chi_{\rm exp}^{2}}=0.8\times 10^{-3}$. This is nothing but the pure generalization error/test error derived after subtracting the experimental noise. Here $\delta\chi_{\rm ML}/\chi_{\rm exp}=0.7$ is the intrinsic machine learning error in the unit of the experimental standard deviation. This is well within the experimental error bar. We show in A that other example of optimization without peak structure shows much larger standard error. If we assume that the inferred $A(k,\omega)$ follows the probability distribution $P=\exp[-\chi_{\rm ML}^{2}/2\chi_{\rm exp}^{2}]$ given from the maximum likelihood inference (see Ref. 21), one can estimate the corresponding variance of the inference for the self-energy by sampling the variation of the peak structure. The variance is plotted in Fig. 17(b) for the peak part of ${\rm Im}W$. This indicates that the variance for the peak position and the weight is small and the existence of the peak is robust.

To further examine the reliability of the emergence of the peak, in Fig. 17(a), we first show the estimated $A^{\rm fit}$ and amplified noise $\tilde{\sigma}_{n}=4\sigma_{n}$ for the optimum doped Bi2212. With $\tilde{\sigma}_{n}$ we generate many synthetic experimental samples. The reason why we take $\tilde{\sigma}_{n}$ instead of $\sigma_{n}$ is to secure the stability of the peak structure in the presence of the experimental noise with the safety factor 4. Then we perform the machine learning and extract the self-energy from the synthetic $A(k,\omega)$, which provides us with the error bars of the estimated self-energies in our machine learning. As shown in Fig. 2(c), the variance of ${\rm Im}\Sigma^{\rm nor}$, and ${\rm Im}\Sigma^{\rm ano}$ thus obtained from the synthetic data is reasonably small with the peak structure in the imaginary part of the self-energy, which indicates that our solution of the inverse problem is numerically stable. Note that the error bars are somewhat overestimated here (namely, larger error bars than those of Fig.2(c)) because of the factor 4 above, but still the peak structure is reasonably retained. However, further increase of the noise to several times of $\tilde{\sigma}_{n}$ smears out the peak structure, implying that very accurate experimental data in the present ARPES quality are required to reveal the peak structure.

The stability in the present inverse problem shown here clearly indicates the difference from notorious ill-conditioned problems such as the analytic continuation from the imaginary time (Matsubara frequency) variable to the real frequency typically studied by the maximum entropy method. In contrast to the non-sparse and involved nature of the analytic continuation from the Matsubara frequency, the present stability numerically shown here is consequences of the sparse structure of the transformation between the spectral function and the self-energy. The mapping between $A$ and $\Sigma$ is, though strongly nonlinear, diagonal in the variable $\omega$ and transformation is sparse confined in a limited frequency range. Imposed physical requirement further ensures the stability. The machine learning including the Boltzmann machine in known to be powerful to strongly nonlinear transformation[9, 14, 69, 70, 7] such as Eqs. (1) and (2).

Of course if the noise is too high, the peak structure is smeared out. Our synthetic data analysis tells that 64 times higher noise level than the estimated experimental noise washes out the peak structure (not shown) and the factor 10 to the present experimental level would be the limit for the meaningful quantitative analysis.

F.2 Dependence on initial guess

The stability of the optimal solution is examined in the previous subsection. While a single local minimum of the cost function is analyzed above, in the present regression scheme, a multi-valley structure of the cost function in the parameter space may appear due to the non-linear nature of the regression. To explore the nature of the possible multi-minima, we perform the outer loop optimization, which is illustrated in Figs. 1 and 3. In the outer loop optimization, we update the initial condition of the imaginary part of the anomalous self-energy to prepare for the next iteration of the inner loop optimization. During the practical optimization, the initial guess for the (R)BM parameters at the first stage of the optimization shown in Fig. 3 will affect the optimized self-energies.

To illustrate the initial guess dependence, here, we examine typical solutions for Bi2212 at 12K obtained from the randomly chosen initial guesses for ${\rm Im}\Sigma^{\rm ano}$. As explained in Sec. III.1, we initialize ${\rm Im}\Sigma^{\rm ano}$ as a linear combination of the Gaussian distributions. We randomly chose the center of mass, height, and width of these Gaussian distributions. From the physical constraint, we choose the center of mass within $|\omega|<$ 0.3 eV. In Fig. 18, the typical examples of the solutions are shown.

There are two kinds of the solutions: Superconducting and pseudogap solutions are found. When an initial guess for ${\rm Im}\Sigma^{\rm ano}$ generates a large enough superconducting gap, superconducting solutions are obtained. The superconducting solutions show a minimum of the amplitude of ${\rm Im}\Sigma^{\rm nor}$ around $\omega=0$, while the amplitude of ${\rm Re}\Sigma^{\rm ano}$ around $\omega=0$ is substantial enough to generate a quasiparticle gap in the spectral function. In contrast, if an initial guess for ${\rm Im}\Sigma^{\rm ano}$ cannot generate a large enough superconducting gap, the amplitude of the normal component ${\rm Im}\Sigma^{\rm nor}$ shows a (negative) peak around $\omega=0$ to generate a gap in the spectral function. We call such a solution the pseudogap solution. A typical example of the pseudogap solution is shown in Fig. 18(a), which shows a three order of magnitude larger cost function as illustrated in Fig. 18(d). Two superconducting solutions with cost functions larger than the minimum value $\chi^{2}_{\rm ML}=1.6\times 10^{-6}$ are shown in Figs. 18(b) and (c). Similarly to the optimal solution shown in Fig. 8(a), the superconducting solution shown in Fig. 18(b) exhibits peak structures of ${\rm Im}\Sigma^{\rm ano}$. However, due to the shift in the peak position, the solution gives an one order of magnitude larger cost function. A featureless ${\rm Im}\Sigma^{\rm ano}$ also generates a superconducting solution with a large cost function as shown in Figs. 18(c) and (d).

Here we note that the larger cost functions originate from systematic deviation of the regression model $A(\omega)$ from the experimental data $A^{\rm exp}(\omega)$. The difference between them, $A(\omega)-A^{\rm exp}(\omega)$, is shown for the three solutions in Fig. 18(e). The pseudogap solution (Fig. 18(a)) and the superconducting solution with the featureless ${\rm Im}\Sigma^{\rm ano}$ (Fig. 18(c)) overestimate the spectral function within the quasiparticle gap: $A(\omega)$ is larger than $A^{\rm exp}(\omega)$ for $\omega\sim 0$ eV. In contrast, the solution with a peak of ${\rm Im}\Sigma^{\rm ano}$ at a higher energy scale shows stronger superconducting gap, which results in $A(\omega)<A^{\rm exp}(\omega)$ around $\omega=0$ eV.

As examplified by the solutions from randomly chosen initial guesses in Fig. 18, the optimal self-energies shown in Fig. 8 indeed give the spectral function closer to the experimental data. Within our many attempts, we found the unique solution that has a cost function value comparable to the estimated experimental error, as shown in Fig. 8.

F.3 Stability against energy cut-off and background

The present machine-learning scheme is based on the imaginary parts of the self-energy within a finite frequency range $-\Lambda<\omega<\Lambda$, where $\Lambda\simeq 0.4$ eV, because the experimental data observed within $-0.4\lesssim\omega\lesssim 0.2$. Therefore, in the genuine self-energy, there is a possible unknown contribution from the outside of the cutoff energy $\Lambda$. However, as explained below, such a contribution is a monotonic and bounded function of $\omega$, and, thus, possible errors due to the lack of information can be estimated.

Due to the Kramers-Kronig relation, the real part of the self-energy can be affected by the cutoff energy $\Lambda$. Because the imaginary part of the normal self-energy is expected to extend over the cutoff energy, the real part of the normal self-energy has a monotonic and bounded contribution from the outside of the cutoff energy. On the other hand, because the anomalous self-energy is finite only within the cutoff energy scale, the real part of the anomalous self-energy can be affected by the cutoff only through the normal self-energy.

In the main text, we ignored the contribution of the high-energy part of normal self-energy. To critically examine the possible contribution from the outside of the cutoff energy, here, we assume a possible distribution of the imaginary part of the normal self-energy outside the cutoff: The imaginary part of the normal self-energy outside the cutoff is assumed to be confined within $\Omega^{\prime}-W^{\prime}/2\lesssim\omega\lesssim\Omega^{\prime}+W^{\prime}/2$ centered at $\Omega^{\prime}$, where $|\Omega^{\prime}|>\Lambda\gtrsim W^{\prime}$ and the amplitude of the imaginary part is approximately constant within this energy range. Then, contribution to the real part is given by

[TABLE]

which is monotonic for $-\Lambda<\omega<\Lambda$. When $|\Omega^{\prime}|>W^{\prime}$ is assumed, the contribution $\Delta{\rm Re}\Sigma^{\rm nor}(\omega=0)$ and its derivative $\left.\partial\Delta{\rm Re}\Sigma^{\rm nor}(\omega)/\partial\omega\right|_{\omega=0}$ are approximately estimated as $\left[-{\rm Im}\Sigma^{\rm nor}(\Omega^{\prime})/(\pi\Omega^{\prime})\right]\left[W^{\prime}+\mathcal{O}({W^{\prime}}^{2}/{\Omega^{\prime}})\right]$ and $\left[{\rm Im}\Sigma^{\rm nor}(\Omega^{\prime})/(\pi{\Omega^{\prime}}^{2})\right]\left[W^{\prime}+\mathcal{O}({W^{\prime}}^{2}/{\Omega^{\prime}})\right]$, respectively. If we consider formation of a lower Hubbard band, for instance, we may assume that $|\Omega^{\prime}|\sim|{\rm Im}\Sigma^{\rm nor}(\Omega^{\prime})|\sim\mathcal{O}(1)$ eV and $W^{\prime}\lesssim\Lambda$. Then, the contribution from the outside of the cutoff is bounded as $|-{\rm Im}\Sigma^{\rm nor}(\Omega^{\prime})W^{\prime}/(\pi\Omega^{\prime})|\lesssim\Lambda/\pi$ and ${\rm Im}\Sigma^{\rm nor}(\Omega^{\prime})W^{\prime}/(\pi{\Omega^{\prime}}^{2})<1/(2\pi)$.

As we show in Fig. 19(a), the qualitative feature of the peak structure at $\omega<-0.4$ eV does not change even when we add an artificial high-energy part for the normal self-energy, where we compare the self-energy inferred from the experimental ARPES data with the self-energy inferred with an additional high-energy part $\Sigma^{\rm H}$. Here, we added the following fixed high-energy part $\Sigma^{\rm H}$ in the optimization process,

[TABLE]

where $s=0.15$eV3, $a=0.3$eV2, $b=0.025$eV2, and $c=0.01$eV2, whose imaginary part becomes substantial for $\omega>-0.4$ eV as shown in the inset of Fig. 19(a). Then we fit the experimental ARPES data within the experimentally measured energy range using this self-energy form with the added high-energy tail. Indeed, in the solution, the artificial high-energy part does not affect the prominent peak structures at all in ${\rm Im}\Sigma^{\rm nor}$ and ${\rm Im}\Sigma^{\rm ano}$. When the energy range covered by the measured $A(k,\omega)$ becomes wider after excluding the extrinsic background, the uncertainty becomes of course further reduced.

To examine effects of high-energy structure of ${\rm Im}\Sigma^{\rm nor}$ further, especially effects of a possible peak structure responsible for the waterfall structure around -0.5eV [72, 73, 74, 75, 76, 77], we conducted the following analysis by adding $\Sigma^{\rm H}$ that has a peak around -0.5eV to the Boltzmann-machine representation of ${\rm Im}\Sigma^{\rm nor}$. Here, we assume a Lorentzian peak around $\omega=$ -0.5eV as $\Sigma^{\rm H}$ instead of monotonic $\omega$ dependence at high energies and find optimized $\Sigma+\Sigma^{\rm H}$ to fit the experimentally observed ARPES data of the optimally doped Bi2212. As shown in Fig. 19(b), the presence of the peak structure around -0.5eV does not essentially change the structures of the low-energy self-energy in the region $\omega>$ -0.2eV. Here, we note that larger amplitude and/or wider width of the Lorentzian cannot fit the experimental spectral function, and we also note that the anomalous cusp in ${\rm Im}\Sigma^{\rm nor}$ near the peak structure around -0.5eV is inevitable to minimize the cost function if one keeps the Lorentzian peak structure. Aside from details, thus, the potential peak structure does not alter the peak structure of our interest. As examined in theoretical and numerical studies (for example, Ref. 78), we also note that the realistic waterfall structure can be derived from broader structures in ${\rm Im}\Sigma^{\rm nor}$ than $\Sigma^{\rm H}$ in Fig. 19(b). In this case the effect of that broad structure becomes smaller than the present critical test.

As an alternative way, one can also examine whether the peak structure is insensitive to the possible slowly varying extrinsic energy dependence or not, by studying the effect of the possible background onto the spectral function. The origin of the background in the experimental spectral function may be the electronic incoherent part or the extrinsic experimental setup extended in the low- and high-energy regions. The background is expected to have a broad (slowly $\omega$ dependent) structure. Such a background has been studied before [71] and we mimic such background structure as a possible effect to see the sensitivity to the peak (dip) structure around $\omega\sim-0.07$ eV. In Fig. 19(c), we show the optimized solution when a model background $b(\omega)$ shown in the inset is subtracted from the spectral function by hand. The high energy offset in the experimental $A(k,\omega)$ does not appreciably depend on temperature and momentum, implying such an extrinsic origin. The way of optimization to minimize $\chi^{2}$ in Eq. (32) is the same as before except that we fit $A^{\rm exp}(k,\omega)-b(\omega)$ with a $\omega$ independent rescaling of the amplitude to satisfy the optimal $n_{0}=0.4$, instead of $A^{\rm exp}(k,\omega)$. The result shows again that such a broad structure does not affect the prominent peak structure and the peak cancellation between the normal and anomalous self-energies is retained. While the cancellation continues to be retained even after the subtraction of the background, the amplitude of the peak structures becomes smaller after the subtraction of the background, as shown in 19(d). Although the reduction of the peak intensity is obtained as a consequence of the regression, it has a simple physical interpretation originating from the renormalization effect: The inferred renormalization factor is affected by the amplitude of the imaginary part of the normal self-energy for $\omega<$ -0.1 eV. After the background subtraction, the amplitude of the imaginary part becomes smaller and, thus, the renormalization factor becomes larger. While the superconducting gap amplitude estimated in the quasiparticle peak does not depend on the background, the normalization factor depends. To reproduce the gap amplitude, the smaller amplitude of ${\rm Im}\Sigma^{\rm ano}$ is required when the renormalization becomes larger. Therefore, the peak structures of ${\rm Im}\Sigma^{\rm nor}$, which have to cancel the peaks in ${\rm Im}\Sigma^{\rm ano}$ to reproduce the spectral function become smaller after the subtraction of the background. By the subtraction of the hypothesized background which is essentially constant for $\omega<-0.1$ eV, we see that the high-energy imaginary part of the normal self-energy has a marginal fermi liquid feature ${\rm Im}\Sigma^{\rm nor}(k,\omega)\propto\omega$ consistent with the coefficient of the $T$-linear resistivity in experiments.

Here we remark that the realistic background estimated from the ARPES data has a similar form to $b(\omega)$, but with much smaller amplitude. As a standard method to estimate the background, we utilize the EDC curve at the momentum far outside the Fermi surface provided by Kondo [79]. The estimated background at the antinodal point can be well fitted by a modified sigmoid function, $w(1-a\omega)/\{1+\exp[(\omega-b)/c]\}$, where $a=2.24$ eV*-1*, $b=-0.02$ eV, $c=0.014$ eV, and the weight of the background $w$ is determined to constitute 20% of the original spectrum for -0.4 eV $<\omega<$ 0.1 eV, which is nearly a quarter of $b(\omega)$ in amplitude employed in Fig. 19(c). Therefore, the analysis in Fig. 19(c) is regarded as the effect of unrealistically exaggerated background as a extreme case. In the realistic background, the modification of the self-energy peak by the background must be much smaller.

ARPES data are basically not available in the unoccupied part and therefore the inferred behavior in the positive energy side of the spectral function has relatively larger uncertainty. In fact, electron-hole asymmetry was suggested in the scanning tunneling microscope data [82]. If any kind of data in the positive energy side can be analyzed together, the asymmetry can be analyzed more quantitatively. However, this issue is beyond the scope of the present work. Nevertheless, the asymmetry expected in the optimally doped Bi2212 (our main target of the analysis) is modest or weak if it exists as one sees from the STM data and the level of quantitative uncertainty in our analysis does not essentially alter the peak structure in the negative energy side. Technically this is analogous to the case of the uncertainty at the tail part of negative energy side beyond the accessibility by ARPES or the effect of the background discussed above. In fact, consideration of the background assumed only in the negative energy side employed here necessarily introduces electron-hole asymmetry in the measurement and the insensitivity to the background is interpreted for the asymmetry as well.

F.4 Bilayer nature

In the bilayer cuprate, the bare band structure consists of the bonding and antibonding band. However, we have analyzed the ARPES data by using the single-band description of the Green function, which describes the bonding band, because it crosses the Fermi energy. We here examine the effects from the antibonding band, which is located above the bonding band in the $\omega$ axis, and show that the contribution from the antibonding band does not change the qualitative results. Below, we explicitly take into account the two-band nature of Bi2212 and demonstrate that the results from the multiband treatment are consistent with those from the single-band treatment.

Due to the symmetry of the Bi2212 crystal structure, the bare band structure is diagonalized by using the bonding and antibonding Wannier orbitals that are given by the bonding and antibonding combination of the $d_{x^{2}-y^{2}}$ orbitals in two adjacent CuO2 planes. Even when the self-energies are taken into account, the Green functions for the bilayer cuprates are diagonal if the interlayer components of the normal and anomalous self-energies are negligible.

When we neglect the interlayer components of the normal and anomalous self-energies, the self-energies in the two adjacent CuO2 planes are identical. Therefore, we only introduce the self-energies, $\Sigma^{\rm nor}(k,\omega)$ and $\Sigma^{\rm ano}(k,\omega)$, for both of the adjacent CuO2 planes. The Green functions for the bonding and antibonding band are obtained as follows, respectively:

[TABLE]

Here, we note that the self-energies are identical even after diagonalization to obtain the bonding and antibonding band, irrespective of the amplitude of the interlayer hopping constituting the bare band.

Then, the spectral function, $A(\omega)$, is decomposed into the contribution from the bonding orbitals, $A^{\rm B}(\omega)$, and antibonding orbitals, $A^{\rm AB}(\omega)$, as follows:

[TABLE]

Here, the orbital dependent matrix elements $\rho_{\rm B}$($=1-\rho_{\rm AB}$) and $\rho_{\rm AB}$ are taken into account.

Then, we extract the self-energies by using the two band decomposition of $A(\omega)=\rho_{\rm B}A^{\rm B}(\omega)+\rho_{\rm AB}A^{\rm AB}(\omega)$. As same as in the single band picture, we train the self-energies by minimizing the training error between the theoretical spectral function $A(\omega)$ and experimental spectral function $A^{\rm exp}(\omega)$. By following Kordyuk et al. [81], we assume that the integrated weight of the contribution from the antibonding band is 25% as $\rho_{\rm AB}=0.25$. In Fig. 20, the optimized spectral function $A(\omega)$ and self-energies are shown.

Even when we take the multiband nature into account, we obtain the self-energies that are essentially same as those obtained in the single band analysis. The cancellation between ${\rm Im}\Sigma^{\rm nor}$ and ${\rm Im}W$ is also observed in Fig. 20b.

Appendix G Details of Wavelet Analysis and Improved

${\bf Im}\Sigma$ to Reduce Test Error

Here, we discuss why we employ this wavelet formalism. The idea of using the binary representation Eq. (16) is to represent a complex function of $\omega$ by successive coarse graining. If a function has $\omega$ dependence with various frequency scales, this hierarchical structure can be efficiently picked up by wavelet with different scales and each wavelet is represented by each digit of the binary number $\sigma_{j}$ $(j=0,1,\dots,L-1)$. For example, the last digit $\sigma_{L-1}$ represents the slowest nonzero modulation of frequency dependence (or in other words, short real-time value), namely with the period of the half of our frequency range. $\sigma_{0}$ picks up the most rapid modulation in frequency (or in other words, long real-time value) alternating in the period of the frequency grid mesh $\Delta\omega=\omega/2^{L}$ etc.. One can have an analogy to Fourier series analysis, where the first digit of the binary number $\sigma_{0}$ corresponds to the largest time component and the last digit $\sigma_{L-1}$ corresponds to the shortest nonzero time component in the form of $e^{i\omega t}$. It was shown that the wavelet can represent the $\omega$ dependence in the orders of magnitude different scales simultaneously and has a flexible representability in the regression problem with small number of parameters because of the logarithmic description (in the present case, $C$ and $D$ with only $L$ arguments, where each of the $L$ components represent logarithmically different scales of frequency dependences) for any discrete data-point set [16, 17, 18].

The grid mesh $\Delta\omega$ is chosen to be smaller than or comparable to the experimental energy resolution ($\sim 10$meV [28, 12]) to fully reproduce the experimental $A(k,\omega)$ within the resolution of the grid size $\Delta\omega$.

If the experimental noise level is comparable to the energy resolution, the step-wise representation is sufficient. However, when the experimental spectral data contains only small noise, the step-wise representation for the imaginary part of the self-energies may introduce a systematic increase in the test errors. To reduce the possible increased error, we introduce a piecewise-linear representation instead of the step-wise representation ${\rm Im}\Sigma^{\rm nor/ano}$. Namely, we interpolate the self-energies between $\omega=\omega_{I}$ and $\omega=\omega_{I+1}$ linearly as

[TABLE]

for $\omega_{I}\leq\omega<\omega_{I+1}$, where $\omega_{I}$ is the midpoint of the I(\mbox{\boldmath\sigma})th interval defined by $\omega_{I}\equiv\Lambda(I+1/2)/2^{L}-\Lambda/2$. Since the estimated noise is very small for the experimental ARPES spectra of the optimally doped Bi2212 ($T_{\rm c}=90$K) at 11K (analyzed in Figures 6(a), 7(a), 7(c), 8(a), 8(c), 8(d) and 9 of the main text; Figs. 10(a), 10(c), 12, 13, 14, 15, 16, 19, and 23 in Appendix), the piecewise linear representation is helpful to achieve the comparable size of test errors with the noise in the experimental data (see Appendix F for the quantitative discussion).

Appendix H Details in Real Part of Self-Energy

The real part of the retarded self-energy is obtained through the Kramers-Kronig relation as Eqs. (28) and (29). For example, the real part of normal self-energy in the stepwise representation derived from Eq. (22) through the discretized Kramers-Kronig relation is given as

[TABLE]

where we introduce a broadening factor $\delta$ to represent a principal value,

[TABLE]

by

[TABLE]

When the piecewise linear representation in Eq. (64) is employed, the real part of the self-energy is corrected by $f_{I}$ as

[TABLE]

where the correction term $f_{I}$ from the stepwise representations in each interval $[\omega_{I},\omega_{I+1})$ is calculated as

[TABLE]

where the width of the interval is given by $\Delta\omega=\Lambda/2^{L}$ and the increment of the imaginary part is $\Delta{\rm Im}\Sigma(\omega_{I})={\rm Im}\Sigma^{\rm nor}(\omega_{I+1})-{\rm Im}\Sigma^{\rm nor}(\omega_{I})$. The amplitude of the contribution $f_{I}$ has the extremum significantly smaller than the increment $\Delta{\rm Im}\Sigma(\omega_{I})$, because

[TABLE]

for $\delta>\Delta\omega$ and thus

[TABLE]

is satisfied when $\delta=10$meV and $\Delta\omega=\Lambda/2^{L}\sim 3.2$meV. Thus, the piecewise linear representation introduces a negligible correction to the real part.

Appendix I Optimization of $n_{0}$

We employed $n_{0}=0.3$ in the present analyses. This choice is justified by the least square fit. We have examined the optimum choice by the least square fit of $\overline{\chi^{2}}$ by taking several choices of $n_{0}$. Here, we note that $\overline{\chi^{2}}$ is trivially scaled by the square of the amplitude of $A^{\rm exp}(\omega)$, and thus is scaled by the square of $n_{0}$. Therefore, we need to optimize $\overline{\chi^{2}}/n_{0}^{2}$. Fig. 21a shows that $\overline{\chi^{2}}$ normalized by $n_{0}^{2}$ for optimally doped Bi2212 at the antinodal momentum has indeed minimum at $n_{0}=0.3$, which indicates that the machine learning suggests that this choice is the optimized value of $n_{0}$. The obtained self-energies do not sensitively depend on the choice of $n_{0}$ as one sees in Figs. 21(b)-(d) and we see no qualitative change in the feature of the pronounced peaks in ${\rm Im}\Sigma^{\rm ano}$ and ${\rm Im}\Sigma^{\rm nor}$ at the same energy together with their cancellation in ${\rm Im}\Sigma^{\rm tot}$.

Appendix J Unoccupied states, bilayer nature, and matrix element

J.1 Unoccupied states

One might be concerned with the effect of the unoccupied states, especially the effect of the particle-hole asymmetry in the scanning tunneling spectroscopy reported in the literature [83, 84]. Our present scheme does not assume the particle-hole symmetry and is capable of the asymmetry if it exists, although the asymmetry at the energies far from the Fermi level may have some uncertainty because of the lack of the data. Aside from the possible origin of the asymmetry [84, 85], the cancellation of the peak structures in the self-energies turns out to be robust.

First of all, in the superconducting phase, the low-energy spectra around the Fermi level are plausibly particle-hole symmetric. Therefore, due to the lack of the information about the unoccupied states, the optimized self-energies are almost symmetric. Then, the asymmetric behaviors could originate from the self-energies away from the Fermi level. However, these high energy asymmetric behaviors hardly affect the cancellation. We have already demonstrated the robustness of the cancellation under the high-energy perturbation in Appendix F. Thus, the essence of the present results is not affected by the asymmetry, which may originate from the high-energy or unoccupied parts of the spectrum. See also Appendix F2.

J.2 Matrix elements

One might also be concerned with the effects of the matrix elements. Even when the energy dependent matrix elements exist, the self-energies obtained by using our method is robust as far as the matrix elements are smooth functions of $\omega$. This is because, when the energy dependent matrix elements are treated as a smooth function of $\omega$, the effects of the matrix element can be simply taken into account as a smooth background. As already examined in Appendix F.3, a smooth background does not change the qualitative and essential conclusion of the present paper. Only if the high resolution ARPES measurements on the entire Brillouin zone, especially at the antinodal region, become available from different experimental setups, quantitative examination of the effects of the matrix elements and background will be conducted and these are highly desirable. However, the data are not available yet and clearly beyond the scope of the present paper.

Appendix K Decomposition of Self-Energy

While the normal and superconducting components of the total self-energy, ${\rm Im}\Sigma^{\rm nor}(k,\omega)$ and ${\rm Im}W(k,\omega)$, show the prominent peak structures, which are absent in the Bardeen-Cooper-Schrieffer (BCS) mean-field theory [86], there are an extended background and a BCS-like superconducting contribution, in addition to the peaks. To highlight the peak part, we decompose ${\rm Im}\Sigma^{\rm nor}(k,\omega)$ and ${\rm Im}W(k,\omega)$ into the peaks and other components. As proposed in Ref. 29, $\Sigma^{\rm tot}(k,\omega)$ may consist of a single pole that generates a superconducting gap and a smooth normal state component. Then, we decompose ${\rm Im}\Sigma^{\rm tot}(k,\omega)$ as

[TABLE]

Here, while the BCS-like superconducting contribution is represented by a Lorentzian,

[TABLE]

where $\Delta_{0}$ and $\Gamma$ are phenomenological parameters that correspond to a BCS-like superconducting gap and life time of quasiparticles, respectively. The background ${\rm Im}\Sigma_{\rm N}(k,\omega)$ can be represented by a linear combination of many Gaussian distributions, where its large amplitude signals either the incoherence of electrons at that energy or some extrinsic origin arising from the experimental setup or background. Then, the peak contribution canceled in ${\rm Im}\Sigma^{\rm tot}(k,\omega)$ is, if it exists, obtained from ${\rm Im}\Sigma^{\rm nor}(k,\omega)$ as

[TABLE]

and from ${\rm Im}W(k,\omega)$ as

[TABLE]

where ${\rm Im}\Sigma_{\rm PEAK}(k,\omega)=-{\rm Im}W_{\rm PEAK}(k,\omega)$ holds.

Appendix L Momentum Dependence

We have shown the machine learning results in the main text at the antinodal point, because the remarkable structure of the pronounced anomalous self-energy peak, coexisting with the normal self-energy peak is most clearly identified with its dominant contribution to the superconductivity. However, the momentum dependence of the peak structure provides us with useful insight. We here show the momentum dependence of self-energy structure. Fig. 22 shows the imaginary part of the normal and anomalous self-energies together with $W$ for the ARPES measurement angle 1.0*∘, 11.1∘, 21.2∘, 31.3∘, and 41.9∘* obtained by the machine learning result of Bi2201 at optimum doping at 11.3 K [12]. Note that 0*∘* is the antinodal and 45*∘* is the nodal points. Although the peaks in the normal and anomalous self-energies become less significant with approaching the nodal point as is expected, the cancellation of ${\rm Im}\Sigma~{}{\rm nor}$ and $W$ always holds and the prominent peak is missing in ${\rm Im}\Sigma^{\rm tot}$ similarly to the case at the antinodal point. This result further corroborates the universal mechanism of the peak cancellation and the dominant contribution to the superconductivity.

Appendix M Machine Learning Results above $T_{\rm c}$

We have shown the machine learning results in the main text for the superconducting phase well below $T_{\rm c}$ to show the remarkable structure of the pronounced anomalous self-energy peak with its dominant contribution to the superconductivity. However, the question how the cancellation of the normal and anomalous self-energy evolves with raising temperatures provides us with further insight on its role in the superconductivity. Fig. 23 shows machine learning results of the normal and anomalous self-energies together with $W$ for the ARPES measurement obtained by the machine learning of Bi2212 at optimum doping at 40, 80 and 120K [12]. The anomalous self-energy peak vanishes above $T_{\rm c}$ as it should be, which further confirms the validity of the present machine learning scheme. ${\rm Im}\Sigma^{\rm nor}(\omega)$ at 120 K (red curve overlapped with the black curve, ${\rm Im}\Sigma^{\rm tot}(\omega)$) shows small signature of pseudogap (small dip around $\omega=0$). The pseudogap, though not a standard behavior as observed in the underdoped region, is clearly seen in the spectrum (23(c)) if one would perform the electron-hole symmetrization and was analyzed in detail in the original experiments (Figs. 17c and 1d in Ref. 12). In Ref. 12, it was even argued that the (incoherent) electron pairing is formed below 150K. The energy of the dip of ${\rm Im}\Sigma^{\rm nor}(\omega)$ shifts with raising temperatures and crosses the quasiparticle peak energy when $T$ crosses $T_{\rm c}$ consistently with the results in Ref. 47. Then the pseudogap is interpreted as generated by the peak of the normal self-energy above $T_{\rm c}$, which is continued from the peak below $T_{\rm c}$, while it is hidden in the spectral function in the superconducting phase below $T_{\rm c}$ because of the cancellation with the anomalous self-energy. It turns out that the $d$-wave superconducting gap has an entirely different origin, which emerges as the sharp drop around $\omega=0$ in ${\rm Im}W$, namely as arising from a pole of $W$, distinct from the pole of the self-energies $\Sigma^{\rm ano}$ and $\Sigma^{\rm nor}$. Although the superconducting gap has a different origin, the main contribution to the superconducting order is attributed to the peak (ideally pole) of ${\rm Im}\Sigma^{\rm ano}$ around $\omega=-0.07(0.04)$ eV for Bi 2212 (Bi2201) (see Figs. 2 and 3). This indicates a tight relation of the pseudogap (manifested by the normal self-energy peak) to the superconductivity (contributed from the peak in ${\rm Im}\Sigma^{\rm ano}$) through their cancellation.

Appendix N

Local and Temperature Insensitive Scattering Rate $z_{\rm qp}(k)c_{1}(k)$

Angle (momentum) dependence of $z_{\rm qp}(k_{\rm F})c_{1}(k_{\rm F})$ plotted for Bi2201 in Fig.4a and resultant $\tau(k)$ shows that it is only weakly dependent around the unity on the angle and doping concentration. (See also Fig. 25 for the plots for each $z_{\rm qp}(k_{\rm F})$ and $c_{1}(k_{\rm F})$.) Even for the optimal Bi2212 at the antinodal point, despite the large difference in $T_{\rm c}$, the value of $z_{\rm qp}(k)c_{1}(k)$ is similar ($\sim 1.4$). (Note that the value $z_{\rm qp}(k_{\rm F})c_{1}(k_{\rm F})$ is somewhat large ($\sim 1.5$) at the nodal point for the underdoped Bi2201 sample, consistently with the increasing slope of the $T$-linear resistivity in the underdoped region [38]. This could be related to the effect of competing insulating behavior.)

With raising temperatures, experimentally observed $T$-linear resistivity has to lead to the temperature insensitive $z_{\rm qp}(k)c_{1}(k)$ at least near the node, because the transport is governed by the nodal region. Temperature dependence of $z_{\rm qp}(k)c_{1}(k)$ in the normal state at $k_{\rm AN}$, for instance for the underdoped Bi2201 is also weak with a large constant offset shown in Fig. 24, which implies that the $T$-linear dependence ($\propto z_{\rm qp}(k)c_{1}(k)T$) is preserved irrespective of the momentum. It supports the local nature of dissipation saturated against temperature below and above $T_{\rm c}$ and intrinsically quantum mechanical.

Appendix O

Momentum and Doping Dependences of $F(k),z_{\rm qp}(k),c_{1}(k)$ and $\Delta_{\rm qp}(k)$

In Fig. 25, doping concentration dependences of the superconducting carrier density $F(k)$, the mass renormalization factor $z_{\rm qp}(k)$, $c_{1}(k)$ defined as the $\omega$-linear component of ${\rm Im}\Sigma^{\rm nor}(k,\omega)$, and $\Delta_{0}(k)$ at the Fermi momentum $k_{\rm F}$ are plotted. Here, $c_{1}(k)$ is defined by

[TABLE]

obtained from the linear fitting of ${\rm Im}\Sigma^{\rm nor}(k,\omega)$ in the range of 15meV$<\omega<40$meV.

Discrepancy between the doping dependence of $T_{\rm c}$ and quasiparticle gap amplitude, established in the literature [42, 43], is further examined by the present self-energy learning. In Table 1, the doping dependences of the density of Cooper pairs $F(k)$, the gap amplitude estimated from the peak position in EDC, the superconducting gap $\Delta_{0}$ obtained by the Boltzmann machine learning, the quasiparticle renormalization factor $z_{\rm qp}$, $c_{1}(k)$, $\bar{Q}$, $P$, $\omega_{\rm PEAK}$ and $\Omega_{0}$ obtained from ARPES data of underdoped (UD), optimally doped (OP), and overdoped (OD) Bi2201 samples [12] are summarized. While the order parameter $F$ and superconducting gap $\Delta_{0}$ show dome-like doping dependence as $T_{\rm c}$ does, the gap amplitude estimated from the peak position in EDC monotonically decreases upon increasing doping [42, 43]. On the other hand, $z_{\rm qp}(k_{\rm AN})$ monotonically increases. They are all consistent with the observed trend in the cuprates.