Identification of non-Fermi liquid fermionic self-energy from quantum Monte Carlo data

Quantum Monte Carlo (QMC) simulations of correlated electron systems provide unbiased information about system behavior at a quantum critical point (QCP) and can verify or disprove the existing theories of non-Fermi liquid (NFL) behavior at a QCP. However, simulations are carried out at a finite temperature, where quantum critical features are masked by finite-temperature effects. Here, we present a theoretical framework within which it is possible to separate thermal and quantum effects and extract the information about NFL physics at T = 0. We demonstrate our method for a specific example of 2D fermions near an Ising ferromagnetic QCP. We show that one can extract from QMC data the zero-temperature form of fermionic self-energy Σ(ω) even though the leading contribution to the self-energy comes from thermal effects. We find that the frequency dependence of Σ(ω) agrees well with the analytic form obtained within the Eliashberg theory of dynamical quantum criticality, and obeys ω2/3 scaling at low frequencies. Our results open up an avenue for QMC studies of quantum critical metals.


INTRODUCTION
Understanding non-Fermi liquid (NFL) behavior near a metallic quantum critical point (QCP) remains one of the most ambitious goals of the studies of interacting electrons. Examples of systems evincing metallic quantum criticality include fermions in spatial dimensions D ≤ 3 at the verge of either spin density-wave, charge density-wave, or nematic-order, 2D fermions at a half-filled Landau level, quarks at the verge of an instability to color superconductivity, and several Sachdev-Ye-Kitaev (SYK)-type models with either electron-electron or electron-phonon interaction  . At a QCP, fluctuations of the corresponding bosonic order parameter become soft. The fermion-fermion interaction, mediated by these soft fluctuations, yields a fermionic selfenergy Σ(ω) ∝ |ω| a with a < 1. The real and imaginary parts of this self-energy are comparable in magnitude and both are larger than ω at low frequencies. This implies that the damping of quasiparticles remains comparable to their energy even infinitesimally close to the Fermi surface, in variance with the central paradigm of Landau's theory of a Fermi liquid (FL). Studies of NFL became the mainstream of research on correlated electrons after a series of discoveries of high-temperature superconductors, which display unconventional metallic properties in the normal state 15,24,44,45 . In most of these materials, superconductivity borders other ordered phases with either spin or charge order. There are also multiple overlaps between the behavior of fermions at a QCP and high-energy physics and string theory 44,46 .
In recent years, several analytical approaches have been developed to study NFL behavior at a QCP. These approaches are based on effective fermion-boson models, in which soft fluctuations of a specific order parameter serve as the source of NFL behavior. The long-standing goal of these studies is to find the functional form of Σ(ω) at a QCP and extract the exponent a < 1 from its small ω behavior. One-loop calculations show that Σ(ω) does become singular at a QCP, for example, in 2D at a transition to nematic or Ising ferromagnetic (FM) order with momentum Q = 0, it scales at the lowest frequencies as ω 2/3 (a = 2/3). Whether this behavior extends beyond one loop is a more tricky issue. Power counting arguments indicate that higher-order terms in the loop expansion for the self-energy reproduce the ω 2/3 scaling form 7 . However, detailed calculations reveal that additional ðlog ωÞ n factors appear, and that n increases with the loop order 25,27,29,31,32,35 . Such logarithms imply that at low enough frequencies, ω ( ω mod , Σ(ω) gets modified from its one-loop form. As a further complication, the same interaction that gives rise to NFL behavior also gives rise to superconductivity at a non-zero T c , so normal state self-energy holds only at ω > T c . It is difficult to extract from analytical studies whether ω mod is larger or smaller than T c , that is, whether the modification of Σ(ω) from higherorder processes is relevant for a metal, which displays superconductivity near a QCP, or only for a putative normal state at T = 0. This uncertainty has triggered an interest in independent numerical studies of the behavior of fermions at a metallic QCP.
Numerical methods for itinerant fermions near a QCP have witnessed great progress in recent years, and at present one can analyze quantum criticality via reliable large-scale numerical simulations 47,48 . In particular, it has been found that designer models of fermion-boson models offer a pathway to access fermionic QCPs while avoiding the notorious sign problem in large-scale quantum Monte Carlo (QMC) simulations. Such models have been implemented in several simulations, studying nematic 49,50 , ferromagnetic 51 , antiferromagnetic 52-56 , gauge field 57-62 , and Yukawa-SYK-type 41 QCPs. The focusing on a particular soft boson offers an unbiased numerical test for either a Q = 0 or a finite Q analytical theory of metallic quantum criticality. The mutual inspiration and dialog between numerical and theoretical communities, arising from these studies, has also stimulated progress along the numerical front (SLMC 63 and EMUS 64 are successful examples of this).
Sign-problem-free QMC has its own limitations as well. To avoid superconductivity and finite size effects, simulations are done at a finite T, which is not the smallest energy scale in the system, such that on a Matsubara axis the fermionic self-energy Σ(ω n ) is a function of a discrete Matsubara frequency ω n = (2n + 1)πT. (The self-energy also has a momentum dependence Σ = Σ(ω n , k), but here and henceforth we suppress this notation for clarity, except where needed.) At non-zero T, it can generally be expressed as Σ (ω n ) = Σ T (ω n ) + Σ Q (ω n ), where the "thermal" part Σ T (ω n ) is the contribution from static thermal fluctuations and the "quantum" part Σ Q (ω n ) is the contribution from dynamical bosonic fluctuations. At T = 0, ω n is a continuous variable, Σ T = 0, and Σ = Σ Q (ω n ) is an NFL self-energy at a QCP. However, at a finite T, the selfenergy differs from its T = 0 form, and the presence of Σ T (ω n ) can mask the behavior associated with Σ Q (ω n ). Besides, at a finite T, Σ Q (ω n ) also generally differs from its T = 0 form. We note that in the Yukawa-SYK model these finite-temperature effects have recently been analyzed using an emergent conformal (reparametrization) symmetry of the low-energy theory, which automatically incorporates thermal and quantum effects 41 . However, the treatment of the present critial FS model without conformal symmetry requires separate analyses of Σ Q and Σ T .
The main purpose of this paper is to provide the method to disentangle Σ T and Σ Q from QMC data for the self-energy. Our approach is based on three observations: • First, to study QC behavior one should avoid the effect of fluctuations from fermions with energies of order of the bandwidth, such as would lead to, or example, Mott physics. For this, the effective femion-boson coupling (labeled g in the text) should be much smaller than the bandwidth W. In systems with a large Fermi surface, W is comparable to the Fermi energy, so the necessary condition is g ( E F .

•
Second, at small g, there is a wide range of frequencies ω n ≪ E F , for which ω n is much larger than Σ(ω n ). In this range, the thermal self-energy has a simple form, valid for finite temperatures and frequencies ω n ≫ Σ, Σ T (ω n ) ≈ α(T)/ω n up to logarithmic corrections, that is, ω n Σ T (ω n ) = α(T) is approximately independent of ω n .
• Third, in the same range, Σ Q (ω n ) still has NFL form and is well approximated by the one-loop, T = 0, expression, modulo that ω n is discrete.
By considering the above points, one arrives at the following conclusion: if a QMC study is performed at g ( E F and provides data for Σ(ω n ) for a substantial number of Matsubara points in the range ω n ≫ Σ(ω n ), it is possible to extract Σ Q from the data by the following simple procedure. First, extract (the approximately constant) α(T) from the data by fitting ω n Σ(ω n ) by a continuous function of frequency and extrapolating to zero frequency, where it is equal to α(T) because ω n Σ Q (ω n ) extrapolates to zero. Once α T is known, subtract Σ T (ω n ) = α(T)/ω n from the full Σ(ω n ) and obtain Σ Q (ω n ), which, as we said, should have the same form as T = 0 selfenergy. For a more accurate separation of Σ Q from Σ T include the slow frequency dependence of α(T) in the fitting procedure, which is still quite straightforward to do, as we will show later.
We apply our strategy to a metal near an Ising FM-QCP. We show the schematic phase diagram in Fig. 1a. It contains regions of a paramagnetic metal (PM) and an ordered Ising FM, separated by a QCP. Right above the QCP, there is a region of small T, where the system displays truly NFL behavior, that is, Σ(ω n ) is nonanalytic and larger than ω n . At higher T, Σ(ω n ) becomes smaller than ω n , yet the self-energy still has non-FL form, and, by our reasoning, its quantum part, Σ Q (ω n ) should be almost the same as at T = 0. In Fig. 1b we show the full self-energy, obtained in QMC simulation, and in Fig. 1c we show Σ Q (ω n ), extracted using the approximate procedure outlined above. The black line in Fig. 1c is the analytical one-loop result for the self-energy at a QCP at T = 0. We see that the data for all ω n nicely fall onto this curve. At small ω n , the analytic one-loop self-energy behaves as ω 2=3 n , and the fact that QMC data fall onto the T = 0 curve implies that the QMC data are consistent with ω 2=3 n scaling at the lowest ω n at a QCP. The deviation from ω 2=3 n scaling in the analytical formula (Eq. (16) in the text) is due to two reasons. First, for the model used for QMC simulations, the bosonic propagator D(q, Ω) contains a regular Ω 2 term along with the Landau damping term, Ω/q. When this term becomes relevant, Σ Q (ω n ) tends to saturate. Second, even when the Landau damping term dominates, the ω 2/3 form is the lowfrequency limit of a more complicated function Σ Q ðω n Þ / ω 2=3 n U ω n =ω b ð Þ, and ω 2/3 behavior holds only when ω n ≪ ω b , that is, UðzÞ % Uð0Þ. The crossover frequency ω b $ ðgE F Þ 1=2 (see Eq. (12) below). In our simulations, this ω b is much larger than the upper boundary of NFL behavior, ω F $ g 2 =E F , but is still much smaller than E F . Accordingly, most of our ω n fall into ω n > ω b , where Σ Q (ω n ) differs from ω 2=3 n . We emphasize that Σ Q (ω n ) has an NFL form regardless of the ratio ω n /ω b . Figure 2 presents a summary of the relevant energy scales in our QMC study.
It is instructive to compare our results with recent analysis of QMC data for similar models. Reference 65 demonstrated that a rather flat dispersion of Σ(ω n ), obtained in QMC simulations, is reasonably well reproduced by Σ(ω n ) = Σ T (ω n ) + Σ Q (ω n ), where both are computed analytically within a metallic QC theory. For that study, a larger coupling g $ E F was used to increase the Fig. 1 Identification of the non-Fermi liquid. a Schematic phase diagram of an FM (2 + 1)D QCP, adapted from ref. 51 . b Fermionic self-energy at an FM (2 + 1)D QCP, calculated by QMC simulation, adapted from ref. 51 . Here, we focus on the point on the Fermi surface with Fermi wavevector along the x-direction. The QMC self-energy appears to have a leading term of the form 1/ω n . c Quantum part of fermionic selfenergy at an FM (2 + 1)D QCP. The black dashed line shows the theoretical prediction of the zero-temperature Fermi self-energy, while the red dashed line marks the low-frequency asymptotic form. We emphasize that the theory is parameter free, and all system parameters, for example, Fermi velocities, are determined separately from the model or QMC measurement. X.Y. Xu et al. magnitude of the self-energy. The discrepancy between the analytic and QMC self-energies in ref. 65 was~20%. This was small enough to see that analytic and QMC self-energies have similar dispersion, but still too high to reliably extract Σ Q (ω n ) from the QMC data. For the current study, g is smaller, and typical Σ(ω n )/ω n is roughly five times smaller than in that work. In this situation, we argue that the QC form of Σ Q (ω n ) can be extracted from the data.
The structure of the paper is the following. In section "The lattice model, phase diagram and QMC self-energy", we describe the lattice model for which the QMC simulations have been performed, and present the numerical results for the self-energy. In section "Analytic self-energy at Ising FM-QCP", we present the analytical results for the self-energy within the self-consistent oneloop analysis. In section "Analysis of QMC data", we extract Σ Q (ω n ) from QMC data and show that for all n > 0 it falls onto the analytic, T = 0 form of Σ Q (ω n ). In section "Discussion", we summarize the results and discuss the implication of this work to other QC cases studied in QMC simulations. We argue that the computational scheme that we proposed can be used as a generic method to extract NFL self-energy at a QCP and can be further extended to study more subtle effects, for example, the flow of the dynamical exponent z.

RESULTS
The lattice model, phase diagram, and QMC self-energy As shown in Fig. 3, we consider a model describing Ising FM fluctuations coupled to a Fermi surface 51 . The model is implemented on a square lattice with HamiltonianĤ ¼Ĥ f þĤ s þ H sf and each part readŝ (1) whereĤ f describes two layers (or orbitals, λ = 1, 2) of spinful (σ = ↑, ↓) fermions with nearest-neighbor hopping on a square lattice, and the chemical potential μ tunes the size of bare Fermi surface. The bare fermion dispersion dictated byĤ f is ϵðkÞ ¼ À2tðcosðk x Þ þ cosðk y ÞÞ À μ and the bandwidth is W = 8t.Ĥ s represents a transverse field Ising model on the same lattice, where by tuning T and h/J an Ising FM to PM transition can be obtained. The onsite coupling termĤ sf between the fermions and Ising spins mediates a fermion-fermion interaction, establishing a metallic system with ferromagnetic fluctuations. We present the schematic phase diagram in Fig. 1a. In the analysis that follows, we focus on the model parameters {t = 1, μ = − 0.5t, J = 1, ξ/t = 1}, for which we find an FM-QCP at h c /J ≈ 3.270 (6). The parameters associated with the fermiology for these parameters are listed in Table 1.
As shown in ref. 51 , our model gives rise to an FM-QCP. However, the bare numerical fermionic self-energy data from QCP, as shown in Fig. 1b, shows a behavior distinctively different from the expected NFL, Σðω n Þ / ω 2=3 n . At low frequency, the self-energy shows an unusual upturn instead of going to zero. Such a upturn in the imaginary part of fermionic self-energy, in the usual numeric setting, implies a gap opening on the Fermi surface. However, our data of the fermionic Green's function does not show a well-formed gap on the FS. Similar behavior of the numerical NFL self-energies has also been observed in other cases, including nematic and AFM-QCPs 49,53 . As discussed in the "Introduction", the rest of this paper is devoted to an analysis of the self-energy data in Fig. 1b, and to understand how to disentangle the thermal and quantum parts of the self-energy, as shown in Fig. 1c.
To understand the situation described in Eq. (1) of itinerant electrons coupled to critical bosonic fluctuations, we can encode the dynamics of bosons and fermions in their propagators, and where k = (ω n , k), q = (Ω m , q) are three-vectors with ω n = (2n + 1) πT and Ω m = 2mπT, the fermionic and bosonic Matsubara frequencies, respectively, ϵ(k) is the dispersion from section "The lattice model, phase diagram and QMC self-energy", M 2 0 represents the bare distance to the QCP before the interaction is turned on (in the QMC it is controlled by the transverse magnetic field, M 0 = M 0 (h)), and Σ, Π are, respectively, the fermionic and bosonic selfenergies. Both self-energies are represented by a diagrammatic series in g ¼ ð ξ 2 Þ 2 D 0 . The series is depicted pictorially in Fig. 4, where solid and wiggly lines are the full propagators G(k), D(q) and the triangles are fully dressed vertices. In general, it is not justified to neglect the vertex corrections. However, it is customary to split the corrections into two types: those coming from fermions away from the Fermi surface ("high-energy" fermions on the scale of the bandwidth W), and those coming from near the Fermi surface. The high-energy contributions just give some static corrections to an effective low-energy theory, which can be absorbed into an effective renormalized coupling g. The condition for the smallness Table 1. Parameters of the fermiology. of these corrections is weak coupling, This condition is valid away from the QCP, that is, when M 2 0 $ k 2 F . In the low-energy theory, at low enough temperatures and frequencies, it is not justified to neglect vertex corrections. However, those vertex corrections that contribute to Σ(k) can be neglected if we are in a regime where |Σ(ω n )| ≪ ω n . As shown in Fig. 1b, the lowest fermionic frequency in our QMC simulation is ω 0 = πT = 0.157 with T = t/20 and the corresponding fermionic self-energy |Σ(k F , ω 0 )| = 0.058, so this condition is satisfied. A longer discussion on this is presented in another work by some of us 65 .
In our QMC study, we are always in the regime |Σ(ω n )| ≪ ω n , and Eq. (4) is obeyed, so without further discussion we will assume that vertex corrections are negligible. Then, Π, Σ are described by the coupled self-consistent equations, ΠðqÞ Here N f is the number of fermion flavors (N f = 2 in the model of section "The lattice model, phase diagram and QMC self-energy") and the factor 2 in Π comes from spin summation.
In principle, Eqs. (5) and (6) have momentum integrals over the entire Brillouin zone, which means that they still include contributions to the self-energies that come from high energies. One of these is a static contribution to Π. This contribution just renormalizes the mass towards the QCP, that is, M 2 0 in Eq. (3) is replaced by Thus, M 2 can be tuned to a QCP by varying g, or alternatively by varying M 2 0 (this is what is done in the QMC simulations). An additional static contribution renormalizes D 0 and we absorb it into g. There are also static contributions to Σ, but they do not change the critical dynamics so we absorb them into the fermionic dispersion. Then, there are dynamical contributions that we will now compute.
Beyond neglecting vertex corrections, we further assume that the fermionic dispersion can be linearized near the FS, which means that the theory describes a low-energy effective theory near the FS. Then, integrating over linearized fermionic dispersion we obtain, where Θ(x) is the step function, the density of states V F ðθÞ ¼ k F ðθÞ=υ F ðθÞ and k F , υ F are the Fermi vector and velocity at an angle θ on the FS, as given in Table 1. In Eq. (8), as |Σ(ω n )| ≪ ω n , we neglected contributions from self-energy and assumed that the k F and υ F vector are approximately parallel. For the fermionic selfenergy we get, ΣðkF; ωnÞ % gT X l Z p dp 2π where σ(x) is the sign function andnðθ k Þ ¼ ðÀυF y ;υF x Þ υF θ¼θ k is an unit vector pointing parallel to the FS at the angle θ k . In a C 4 symmetric system, we can replacen by υ F =υ F , since the unit vector only determines the value of υ F (θ) in Eq. (8). We first evaluate the bosonic self-energy which to leading order is, where the C 4 symmetry of the lattice is used to replace υ F (θ ± π/2) = υ F (θ) and similarly for V F . Next, we turn to the fermionic self-energy. Plugging Eq. (10) into Eq. (9) yields, In Eq. (11) we rescaled momentum to frequency ω = υ F p, and hid the explicit angular dependence υ F = υ F (θ k ) for conciseness. The frequency scale introduced by Π is From Eq. (11) we can read the relevant frequency scales for Σ (ω n ). The typical scale of the ω l sum is ω l~ωn due to the sign function, that is, typical internal frequencies are constrained to be on order of the external frequency.
We now show that at finite-temperature, but as long as |Σ(ω n )| ≪ ω n the fermionic self-energy in Eq. (5) splits into two parts: thermal and quantum (for detailed derivations and discussions see, e.g., refs 16,21,34,43,65 ). The quantum part recovers the zero-temperature fermionic self-energy, while the thermal part takes on a very simple form and scales as 1/ω n . Thus, after simply deducting this 1/ω n term, the finite-temperature self-energy directly provides the zero-temperature behavior of fermions, although the measurement is done at finite temperature, at which thermal fluctuations has a significant contribution. This is one of the key conclusions of this work. We separate the summation in Eq. (11) into two parts where Σ T is the ω l = ω n piece of the sum in Eq. (11), namely where As SðxÞ vanishes rapidly at large x, it predicts that Σ T only contributes significantly at finite temperature and close enough to the QCP (πT ≳ υ F υ F M). In that regime, as noted in the introduction, α(T, ω n ) = ω n Σ T (ω n ) depends at most logarithmically on frequency at the smallest ω n , α(T, ω n ) ≈ α(T).
The quantum part includes all other terms in the Matsubara sum. This sum can be approximately replaced by an integral, which immediately recovers the T = 0 form of the fermionic Fig. 4 The diagrammatic representation of bosonic self-energy Π(q) and fermionic self-energy Σ(k). X.Y. Xu et al. self-energy, that is, with UðzÞ ¼ The scaling function UðzÞ has the following asymptotics, : (18) where z À1 c ¼ ðυ F =cÞ 3=2 and u 0 is a constant which depends on υ F =c. For υ F =c ≪ 1, u 0 ≈ 1/8; while for parameters of section "The lattice model, phase diagram and QMC self-energy" (υ F =c ≈ 0.42), u 0 ≈ 0.1. Note in the case of our QMC study z c ≈ 3.7, so that the intermediate regime cannot really be seen. Equation (17) is exactly the formula we used to generate the black line in the Fig. 1c, and is the quantum NFL self-energy Σ Q of an FM-QCP. It saturates in the large frequency region as shown in the figure, as predicted by Eq. (18) in the z c ≪ z limit. Combining Eqs. (14) and (16), we indeed see that the self-energy has a thermal 1/ω n term plus the zero-temperature quantum self-energy.
Let us briefly elaborate on the physics behind the scaling function UðzÞ. In Eq. (17), the part in the square brackets correspond to the fermionic propagator and the other part in the integral corresponds to the bosonic propagator. Consider the limit ω ≪ ω b corresponding to z ≪ 1. Expanding for z ≪ 1 we find that to leading order the terms in the square brackets are a constant, and the dx integral is limited to 0 < x < 1. Physically this is the statement that the momentum integration (∫dy) is only on bosonic momentum parallel to the FS. In addition, the ðυ F =cÞ 2 x 2 y term in the boson propagator is also negligible, which corresponds to the fact that the bare Ω 2 part of the boson dynamics is irrelevant at low frequency. Evaluating Eq. (16) for ω n ≪ ω b we find, where Equation (19) is the formula used to generate the red dashed line in Fig. 1c, as an asymptotic line of the quantum part of the selfenergy predicted by Eq. (17). The analysis of Σ leading to Eqs. (19) and (20), as well as analogous analysis for superconducting selfenergy, is conventionally termed "Eliashberg theory" (ET), due to its similarity to ET of superconductivity from electron-phonon interactions 73 . Now consider the opposite limit, ω ≫ ω b , corresponding to z ≫ 1. For simplicity let us assume υ F /c ≪ 1. In that case the term in the square brackets, corresponding to the fermionic propagator, is not constant, and the bulk of the contribution to u 0 is given by the range 1 < x < ∞. Physically this means that scattering is not confined to be parallel to the FS and is two dimensional, although it is still confined to be near the FS. It is instructive to compute the subleading term for small z. After some algebra, one finds that this contribution is also given by 2D scattering, and gives ð2π= ffiffi ffi 3 p ÞUðzÞ % 1 À 0:73z 1=3 . This means that for z~1, Σ 0 is reduced by a factor of almost 4 from the expected value if one considers only the leading contribution. This is the reason that the deviation from the asymptotic red line in Fig. 1c is so large. At even larger z, the Ω 2 term in the bosonic propagator begins to contribute, which just modifies the high-frequency behavior of the self-energy. However, the deviations from ω 2/3 scaling occur already at z~1. We term the theory which accounts for both highfrequency modifications and the finite-temperature corrections of Eq. (14) a modified Eliashberg theory (MET).
Analysis of QMC data Now we turn to the QMC data analysis. We study the fermionic self-energy from the FM-QCP model described in section "The lattice model, phase diagram and QMC self-energy" and compare the QMC data with the MET in section "Analytic self-energy at Ising FM-QCP".
Let us begin by going through the relevant physical parameters in the QMC data. We normalize all quantities by the hopping energy t = 1, see section "The lattice model, phase diagram and QMC self-energy" for details. By tuning the transverse field h, ref. 51 was able to extract QMC data for different T above the QCP, and also deep in the disordered phase, where the self-energy should have an FL form Σ(ω n ) ∝ ω n at low frequencies. We concentrate on the data at the QCP. The parameters from section "The lattice model, phase diagram and QMC self-energy" imply a bare g ¼ 1=4, which is much smaller than E F ≈ 1.6, implying the QMC is in the weak coupling regime (we remind that all energies are quoted in units of the hopping). The bosonic propagator in QMC was found to agree well (as shown in ref. 51 , it turns out the bosonic propagator is dominated by the bosonic self-energy part, with a small finite anomalous dimension in q 2 and Ω 2 terms, it will not change the main results of this paper) with Eqs. (3), (7), and (10), that is, with D 0 = 1, M 2 (T) = 0.13T 1.48 , c = 3.16, that is, the measured D 0 agrees with the bare one, all obtained from the bosonic propagator data in ref. 51 . In addition, it was found that Σ(ω n ) ≪ ω n for all temperatures and Matsubara frequencies that were obtained. Thus, we may expect that corrections to the bare g are small, and the renormalized g, which is an input to MET is at the order of the bare one. Under this condition, the relevant scales for Σ are The temperatures we analyze are T = 0.05, …, 0.1, which implies the first Matsubara frequency is πT = 0.16, …, 0.31, see the schematics of energy scale in Fig. 2. Thus, ω F is completely irrelevant as is verified by the fact that the self-energy is always small. As we discussed in section "Introduction", the QMC selfenergy appears to have a leading term of the form as shown in Fig. 1b. This is consistent with the prediction of MET, see Eqs. (14) and (15). We analyze the data in two ways. First, we extract the quantum self-energy and compare it to the T = 0 prediction. To do this, we need to remove the thermal part. This is most conveniently done simply by studying the product ω n Σ(ω n ). As discussed in the "Introduction" and the previous section, according to Eqs. (13) and (16) we have, providing we treat α(T) as constant, neglecting its slow frequency dependence, see Eq. (15). In Eq. (24), α(T) includes both the contribution from Σ T (ω n , T), and corrections from finite size effect (such as a possible small gap due to the mismatch of finite size h c X.Y. Xu et al. and the thermodynamic h c ). The second part, that is, ω n Σ Q (ω n ), comes from the MET prediction for Σ Q (ω n ), Eq. (16), which recovers ET prediction Eq. (19) in the low-frequency limit (ω n ≪ ω b ). We fit Eq. (24) to the data for all T simultaneously. Importantly, in Eq. (24), the fitting parameters are only the constants α(T) and g. This is because Σ Q is a function only of g and system parameters, see Eq. (17). Figure 1c from the beginning of our paper depicts the result of our fit. We obtain a fitting of g ¼ 0:245 ± 0:023 for 95% confidence intervals, in excellent agreement with the theory. Regarding α(T), we find that α(T) ≈ 8 × 10 −3 is almost a constant, in disagreement with the expected ∝ T behavior of ω n Σ T (ω n , T). Clearly, part of this discrepancy is due to our neglecting the frequency dependence of α. We therefore repeat the analysis using the following fitting procedure, which takes the full frequency behavior of Σ T into account. Guided by the previous fit, we set g ¼ 0:25 to be the bare one to reduce the number of fitting parameters. We show the result of this fit in Fig. 5 and the extracted α 0 ðTÞ in Fig. 6. The agreement is very good, and we checked that the data collapse can be made even better by allowing g to vary somewhat (equivalent to about 13% change in the bare vertex ξ). The extracted α 0 ðTÞ indicates the formation of a small gap forming at around T = 0.1, which is expected to yield a self-energy contribution of the form α'(T)/ω n = Δ 2 (T)/ω n . The gap size Δ corresponding to α 0 ðTÞ is much less than the numerical inverse reciprocal lattice spacing, so the appearance of this gap is actually an expected effect, which however is beyond the resolution of the standard methods for verifying the appearance of long-range order. Thus, our analysis of the self-energy yields a method for more accurately finding the QCP in our system.
Here we add a word of caution. Previous work has shown 33,38,74 that the first Matsubara frequency does not obey the quantum critical scaling Σ Q (πT) ∝ (πT) 2/3 , and therefore should not be included in the fitting procedure. We verified that dropping the first Matsubara point does not change our results. Also, note that within the error range in Fig. 5, it is possible that Σ Q (πT) < 0. In fact, it can be verified that Σ Q (πT) is always negative 21,65 .
To avoid this issue, we also numerically computed Σ T (ω n ) and Σ Q (ω n ) by performing the Matsubara sum in Eq. (9), using g ¼ 0:25. This procedure takes into account the full frequency dependence of Σ T (ω n ) as well as finite mass effects and the first Matsubara frequency issues. Figure 7 depicts a comparison of the QMC self-energy with the numerical summation. There is an excellent agreement between the two, except for a T-dependent constant offset between the MET and QMC results. The result is consistent with the first analysis we performed above. For completeness, we also performed a comparison between the MET and QMC data for the data in the disordered phase (the FL regime). Figure 8 shows this comparison, again with very good agreement.
We therefore conclude that we have extracted the quantum self-energy from the QMC data, and that it shows excellent agreement with the expected QC behavior.

DISCUSSION
NFLs play a crucial role in a wide range of quantum many-body phenomena, such as quantum criticality, high-temperature superconductivity in correlated materials, unconventional transport in strange metals, and have been a key focus in the study of modern condensed matter physics [1][2][3][5][6][7]12,[14][15][16][17][18][19]22,24,25,[27][28][29]33,35,36,44,45,47,48,53,[66][67][68][69][70][71][72]75 . Despite of the intensive research efforts, key questions remained open and Fig. 5 Extraction of the quantum self-energy. The solid dots correspond to the QMC ω n Σ(ω n ) for different T, while for each T dataset, the thermal part and a constant α'(T) has been deducted (see Fig. 6). The dashed line corresponds to ω n Σ Q (ω n ), computed at T = 0, for the bare g from the parameters of section "The lattice model, phase diagram and QMC self-energy". The gray shaded area is the 95% confidence interval. Fig. 6 The extracted gap contribution to ω n Σ(ω n , T). See Eq. (25) for details. Fig. 7 Comparison of the full self-energy between MET and QMC at the QCP. The solid dots correspond to the QMC data, and the hollow dots correspond to a numerical summation of the Matsubara sums in Eq. (11). Fig. 8 Comparison of the full self-energy between MET and QMC in the FL region (h/J = 3.6 > h c /J). The solid dots correspond to the QMC data, and the hollow dots correspond to a numerical summation of the Matsubara sums in Eq. (11).
the problem of NFLs is still one of the most challenging topics in many-body physics, even with the most sophisticated field theoretical treatments 16,25,27,29,35 , powerful numerical many-body algorithms and high-performance supercomputers 47,48,63,64 .
Our work provides a pathway to address a key challenge in the study of NFLs, that is, the fact that the smoking-gun signature of NFLs (the predicted unconventional low-temperature fermion self-energy), has never been directly observed or verified in large-scale unbiased numerical methods. Through combined numerical and theoretical efforts, we proved that this key signature of NFLs can be accessed through QMC simulations, by simply deducting a ∝ 1/ω n thermalfluctuation background. This technique enabled us to directly compare numerical results with theoretical predictions, providing a bridge between theoretical, numerical, and experimental studies.
Although this paper mainly focuses on the itinerant ferromagnetism QCP as an example to demonstrate the physics, the technique is universal and can be easily generalized to other itinerant QCPs, such as nematic and AFM-QCPs 49,52,53,56 . Furthermore, this technique can also be used to explore the predicted non-trivial effects from higher-order corrections 16,25,27,29,35,76 , and thus open up a pathway towards a full understanding about this challenging subject of NFLs.

Numerical calculations
The numerical results for fermionic and bosonic self-energies have been obtained using state-of-art determinantal QMC simulations as reported in ref. 51 .

Analytical calculations
Analytical calculations have been carried out diagrammatically within ET and MET, by solving the set of self-consistent equations for fermionic and bosonic self-energies.

DATA AVAILABILITY
The data that support the findings of this study are available from the first author upon reasonable request.