Hierarchy of energy scales in an O(3) symmetric antiferromagnetic quantum critical metal: a Monte Carlo study

We present numerically exact results from sign-problem free quantum Monte Carlo simulations for a spinfermion model near an O(3) symmetric antiferromagnetic (AFM) quantum critical point. We find a hierarchy of energy scales that emerges near the quantum critical point. At high energy scales, there is a broad regime characterized by Landau-damped order parameter dynamics with dynamical critical exponent z = 2, while the fermionic excitations remain coherent. The quantum critical magnetic fluctuations are well described by HertzMillis theory, except for a T −2 divergence of the static AFM susceptibility. This regime persists down to a lower energy scale, where the fermions become overdamped and concomitantly, a transition into a d−wave superconducting state occurs. These findings resemble earlier results for a spin-fermion model with easy-plane AFM fluctuations of an O(2) SDW order parameter, despite noticeable differences in the perturbative structure of the two theories. In the O(3) case, perturbative corrections to the spin-fermion vertex are expected to dominate at an additional energy scale, below which the z = 2 behavior breaks down, leading to a novel z = 1 fixed point with emergent local nesting at the hot spots [Schlief et al., PRX 7, 021010 (2017)]. Motivated by this prediction, we also consider a variant of the model where the hot spots are nearly locally nested. Within the available temperature range in our study (T ≥ EF /200), we find substantial deviations from the z = 2 Hertz-Millis behavior, but no evidence for the predicted z = 1 criticality.


I. INTRODUCTION
Quantum criticality in itinerant many-fermion systems is of great importance in the study of strongly correlated materials. Due to strong inherent correlations, many classes of materials, such as the cuprates [1], iron-pnictides [2], organic superconductors [3], and heavy-fermion compounds [4] host rich phase diagrams with magnetic order, non-Fermi liquid regimes, and strange metal transport. It is generally believed that quantum critical points (QCPs) at T = 0 are playing a role in these systems [5][6][7][8]. Perhaps the most prominent feature in those materials is an extended phase of high temperature superconductivity, which is generally observed in the vicinity of antiferromagnetic order. It has been proposed [9][10][11][12][13][14][15][16] and numerically demonstrated [17][18][19][20][21] that nearly-quantum critical fluctuations can mediate unconventional superconductivity and anomalously enhance the transition temperature T c . However, an in-depth understanding of this universal mechanism still remains highly desirable.
Despite of numerous analytical attempts over the past decades [5,7,15,[22][23][24][25] few definite statements about the critical properties could be made due to the inherently non-perturbative nature of the interactions between gapless bosonic (magnetic) and fermionic (electronic) degrees of freedom. An audacious line of research [26][27][28] carried out by Sung-Sik Lee and collaborators, however, has recently made notable progress in identifying a route to a controlled expansion in an emergent control parameter to compute the properties of a strongly coupled fixed point at the O(3) antiferromagnetic QCP in an exact manner. One concrete prediction of this non-perturbative analysis is a low temperature regime exhibiting a dynamical scaling exponent z = 1 [26][27][28]. At this fixed point, the tendency towards superconductivity is strongly suppressed. This general prediction calls for a complementary numerical study inspecting the hierarchy of energy scales in the vicinity of an antiferromagnetic metallic QCP.
It is the purpose of this manuscript, to provide precisely this type of comprehensive numerical insight for an elementary microscopic model -the spin-fermion model [24,28,29] capturing the interplay of itinerant electrons and quantum critical fluctuations of an antiferromagnetic O(3) spin density wave (SDW) order parameter. We employ large-scale quantum Monte Carlo simulations, using a controlled and unbiased Monte Carlo flavor, determinant quantum Monte Carlo (DQMC) [30][31][32][33][34], which does not suffer from the infamous sign problem [35]. On a technical level, it is the presence of an The black lines correspond to the Fermi surfaces of the two fermion flavors ψx (solid), ψy (dashed). One band has been shifted by Q = (π, π) such that hot spot pairs (red) occur at crossing points. The energy of the ψx band is shown as background (color shading).
antiunitary symmetry [36] (similar to time reversal) that renders the effective "Boltzmann weights" in the partition function strictly positive semidefinite. Building upon this insight, we calculate numerically exact results for the O(3) antiferromagnetic QCP in a two-dimensional metal.
Mapping out the phase diagram (shown in Fig. 5 below), we find an extended region of d-wave superconductivity while other electronic orders, like charge density wave (CDW), do not proliferate. Above T c , we find a broad temperature regime characterized by quantum critical scaling with dynamical exponent z = 2 and overdamped dynamics of the SDW order parameter, while the fermionic excitations remain underdamped. As in the case of the recently studied easy-plane XY order [17,37], the antiferromagnetic susceptibility fits a Hertz-Millis form [24,25], except for its temperature dependence, which is consistent with T −2 . This is remarkable as it is well understood that Hertz-Millis theory is a formally uncontrolled approach [15,23]. At temperatures directly above the QCP, the imaginary part of the fermionic self-energy is approximately constant at as a function of frequencies and temperatures, indicating a breakdown of Fermi liquid theory.
Finally, we search for the regime described by Schlief et al. [26], which is characterized by a z = 1 dynamical scaling and strongly spatially anisotropic SDW order parameter correlations. To that end, we conduct an ultra-low temperature (down to T ∼ E F 200) DQMC study for a Fermi surface tuned close to local nesting at the hot spots -a scenario that is expected to bring us closer to the z = 1 fixed point [26,27]. Indeed, we observe strong deviations in the Matsubara frequency and momentum dependence of the SDW propagator from the Hertz-Millis form. However, we neither observe a dynamical critical exponent z = 1 nor a momentum-anisotropy in the propagation of magnetic modes for our model parameters. This leaves us with the conclusion that, within the resolution of our DQMC simulations for systems of linear size L ≤ 14 and temperatures T ≥ E F 200, we cannot numerically observe the new fixed point discovered by Schlief et al. [26] in our O(3) spin-fermion model.
Our discussion in the remainder of the manuscript is organized as follows. In Section II we introduce a spin-fermion model for a metal hosting an antiferromagnetic QCP. In Section III we provide an overview of the theoretical background from the perspective of a perturbative expansion in the spinfermion coupling. We discuss a hierarchy of relevant energy scales separating different physical scaling regimes and highlight, in particular, the role of vertex corrections in this context. We then turn to an extensive DQMC analysis in Sec. IV in which we map out the phase diagram of the metal for FIG. 2. Electron hopping and local interactions on the square lattice. Illustration of nearest (green), next-nearest (blue), and nextnext-nearest neighbor (red) hopping processes (a) and local Yukawa interaction (dashed lines) of the two fermion flavors ψx, ψy and the magnetic order parameter φ (b) present in model (1). different coupling strengths before examining magnetic and metallic correlations near the critical point in detail. In Section V we investigate the effect of local Fermi surface nesting on scaling properties by means of ultra-low temperature DQMC simulations. Finally, we close by discussing our findings and commenting on potential future research directions in Section VI.

II. MODEL
In the vicinity of a metallic SDW quantum critical point, the relevant low-energy degrees of freedom [7,23,50] are the antiferromagnetic collective modes and electrons at the hot spots, a set of points on the Fermi surface where spin-fermion scattering is resonant, see Fig. 1. To describe the interplay between these gapless fermionic and bosonic modes, we consider a low-energy effective action S = S ψ + S λ + S φ on the two-dimensional square lattice, in which two flavors of spin-1 2 fermions are coupled to an isotropic three-component SDW order parameter ⃗ φ which tends to order [35] at a commensurate antiferromagnetic ordering wave vector Q = (π, π), Here r, r ′ are coordinates on the square lattice, s ∈ {↑, ↓} denotes spin, α ∈ {x, y} is a fermion flavor index, ⃗ σ is a vector of the Pauli matrices, and ∇ is the lattice gradient. Imaginary time is denoted as τ and the corresponding integrals implicitly go from zero to inverse temperature β = 1 T .
The first part, S ψ , describes the kinetic energy of the fermions ψ αs . SDW ordering is modelled by a canonical φ 4 -theory, S φ , which couples to the fermion spin density via a Yukawa coupling in S λ . We set the bare velocity of the bosonic collective mode to c = 3, the quartic coupling to u = 1, and the fermion-boson Yukawa coupling to λ = 1 (unless stated otherwise). The parameter r, related to the mass of the order parameter field, allows us to tune through a SDW transition.
Specifically, unless otherwise noted, we consider nearest neighbor hopping amplitudes t xh = t yv = 1, t yh = t xv = 0.5 where subscripts h, v indicate horizontal and vertical lattice directions, respectively. The chemical potential has opposite sign for the two fermion flavors, µ x = −µ y = −0.5. The resulting Fermi surface is displayed in Fig 1. Note that for this choice of parameters model (1) has a two-fold rotational symmetry, consisting of a combination of a π 2 real-space rotation, shifting the momentum by Q, a particlehole transformation, and interchanging fermion flavors, i.e. ψ αrs → se iQ⋅r ψ † −αR π 2 (r)−s . This symmetry implies that hot spots occur on the momentum diagonals, k x = k y .

III. THEORETICAL BACKGROUND
We start our discussion of the spin-fermion model (1) by providing a brief overview of its relevant energy scales, summarizing some essential aspects of its theoretical background along the way. To identify a number of crossover regimes and to understand the basic impact of the spin-fermion interactions, it is instructive to perform an analysis for small Yukawa couplings λ -despite the fact that, strictly speaking, the Yukawa coupling is a relevant parameter in the sense that conventional perturbation theory becomes unreliable at sufficiently small energy scales.
Focusing on the antiferromagnetic order parameter first, we can extract an energy scale Ω b , at which the nature of the dynamics of collective spin excitations changes due to interactions with surrounding fermions. To leading order in λ, the propagator of the bosonic field φ takes the renormalized form in whichr is the renormalized tuning parameter and is a damping constant. Here, N = 4 is the number of hot spot pairs, θ ∈ [0, π] is the angle between the Fermi velocities (of magnitude v F ) at the hot spots (Fig. 4). By estimating when the dynamically generated contributions to the action dominate over bare ones one then finds [5,15,18,22] for the energy scale Ω b At frequencies small compared to Ω b the dynamics of collective magnetic modes is overdamped. Because of the emergence of the Landau-damping term γ ω n in Eq. (2) the dynamical critical exponent z is found [5,7,15,18,22] to increase from 1 (at the Wilson-Fisher fixed point) to 2. Next, we consider the lowest-order effect of the renormalized bosons on the fermion dynamics at the hot spots. To leading order in λ one finds [15,22] that, at zero temperature, the fermions acquire a self-energy where n b = 3 is the number of order parameter components. This implies that the fermions become strongly damped by the coupling to SDW fluctuations, with a damping rate that scales as √ ω, indicating a distinct deviation from ordinary Fermi liquid character (ω 2 log(1 ω) in two spatial dimensions). Taking the same approach as above, we estimate the energy scale Ω f at which this breakdown occurs as For frequencies ω ≪ Ω f the feedback from 'dressed' bosons on fermion excitations at the hot spots is strong and leads to quasi-particle decoherence. Apart from entering a non-Fermi liquid regime, the system may eventually become unstable against superconductivity at sufficiently low temperatures. Close to the critical point collective spin fluctuations can generate attractive pairing interactions and facilitate the formation of Cooper pairs [17,28]. At the hot spots, and for small λ, one can estimate the energy scale of the onset of superconductivity, T c , by solving the Eliashberg equation for the superconducting vertex in the Landau-damping regime, where self-energy corrections are negligible. For our model (1) this predicts [18] the superconducting susceptibility to diverge at a scale Note that the energy scales T c and Ω f are both of the order O( λ 2 N ) and thus there is a priori no parametrically large  separation between a non-Fermi liquid and a superconducting regime.
We summarize this hierarchy of energy scales obtained from a small λ analysis in Fig. 3. We note that the scales Ω b and Ω f are both of order λ 2 in the Yukawa coupling. This is in stark contrast to the case of an Ising-nematic QCP (where Ω b ∼ λ and Ω f ∼ λ 4 ) which shows a parametrically large window of Landau-damped physics (with z = 3) [19,38]. However, the energy scales are separated by their dependence on the number of hot spot pairs, N = 4, which may serve as a control parameter for an extended Fermi liquid regime with relaxational boson dynamics.
Finally, we inspect first-order corrections δΓ to the spinfermion vertex itself. Following a similar procedure, one can extract an energy scale Ω λ at which the dynamically generated higher-order interactions becomes comparable to the bare vertex, δΓ ∼ λ. One obtains [15] where Ω uv is an ultra-violet cutoff and 0 < θ < π. Note Ω λ does not depend on the coupling strength λ -the naive λ 2 is compensated by the Landau-damping constant γ. Instead, for a fixed number of hot spot pairs, the scale Ω λ is only sensitive to the relative angle θ between the Fermi surfaces at the hot spots. We illustrate the dependence of Ω λ on θ in Fig. 4. In particular, we see a moderate increase of Ω λ in the limit of (local) nesting, θ → 0. Let us remark that Eq. (8) is special to the case of an O(3) symmetric SDW order parameter. Generically, the leading order correction to the boson-fermion vertex scales as δΓ ∼ 2 − n b , and thus precisely vanishes for easy-plane antiferromagnetism.
Despite the general guidance provided by such a perturbative treatment of the Yukawa coupling, this approach is formally uncontrolled and must be complemented by a more sophisticated analysis. Many intensive analytical efforts have been made to get a handle on strong interactions and gain control over calculations. These include various 1 N -expansions [5,15,22] and extensions of the problem to fractional dimensions [27,51]. One of the more striking developments has been a recent, self-consistent and non-perturbative study which has been put forward by Schlief et al. [26]. It utilizes an emergent control parameter -the degree of local nesting at the hot spots -to compute critical exponents in an exact manner. Notably, the identified strong coupling fixed point is characterized by the following form of the SDW susceptibility: where γ ′ is a parameter. Note that the SDW correlations are highly anisotropic in space, and have a dynamical exponent z = 1. In this regime, the fermions remain underdamped. Clearly, a confirmation of those predictions with an unbiased method would be highly desirable and the motivation for the numerical studies to be reported in the following. A central open question is whether the revealed z = 1 fixed point is a generic property of a SDW QCP, in the sense of an extensive basin of attraction, or rather limited to a narrow parameter range.

IV. DETERMINANT QUANTUM MONTE CARLO
We investigate the physics of model (1) by means of extensive finite-temperature DQMC [30][31][32][33][34] simulations. The key feature of these quantum Monte Carlo methods is that they are numerically exact -given sufficient computation time both the statistical Monte Carlo error and the systematic Trotter error can be made arbitrarily small -and do not rely on the smallness of any expansion parameter. Provided the absence of the famous sign-problem, they allow us to explore the relevant region of the exponentially large configuration space in polynomial time. For our two-flavor spin-fermion model, it has been proven that probability weights are strictly positive because of the presence of an antiunitary symmetry similar to time reversal [35].

A. Phase diagram
Before turning to quantum critical properties, we show the phase diagram of model (1) for coupling, λ = 2, in Fig. 5. In two dimensions, the O(3) symmetry present in S φ cannot be spontaneously broken at finite temperatures according to Mermin and Wagner's theorem [52]. Magnetic correlations, indicated by the SDW susceptibility in which we measure momenta relative to the antiferromagnetic ordering wave vector Q = (π, π), are therefore bound to be finite-ranged at any nonzero T and there is no classical phase transition originating from the QCP at T = 0. Anticipating that, at large length scales, the susceptibility is of  1), for λ = 2 (left) and λ = 1 (right). The correlation length ξAFM is shown for a L = 12 system (color coding, paramater grid indicated by grey points). In accordance with the Mermin-Wagner theorem there is no magnetic phase transition at finite temperatures but only a crossover originating from the QCP at r = rc. For large tuning parameter values r ≫ rc the system is in an ordinary Fermi-liquid phase. In the opposite limit, r ≪ rc the SDW correlation length diverges as T → 0. For λ = 2, an extended dome-shaped d-wave superconducting phase is masking the QCP (blue). The maximal Tc is of the order of E F 20. For λ = 1, the system stays non-superconducting down to the lowest considered temperature, T = 0.025. Depending on the coupling strength, the quantum critical point, marking the onset of SDW order at T = 0, lies at rc ≈ 3.85 (λ = 2) or rc ≈ −1.89 (λ = 1), respectively.
Ornstein-Zernike form [53], we can define a correlation length of antiferromagnetic fluctuations [53,54] ξ AFM = L 2π with q = 2π L, the smallest non-vanishing momentum of a finite system. In Fig. 5 we display ξ AFM across the (r, T )-plane using a linear interpolation between a discrete set of r values and temperatures (indicated by dots in the figure). As can be clearly seen, the data manifestly reveals a finite-temperature magnetic crossover culminating into a quantum critical point marking the onset of SDW order at T = 0. While the Mermin and Wagner theorem also holds for the continuous U (1) symmetry associated with superconductivity there is the possibility of a Berezinskii-Kosterlitz-Thouless (BKT) transition for the two dimensional superconducting order parameter. We determine the superconducting transition temperature T c (r) as the temperature where the superfluid density of the system surpasses the universal BKT value ∆ρ s = 2T π [17,32]. The resulting phase transition is shown in Fig. 5. For λ = 2, we observe an extended dome-shaped superconducting phase that shows its peak critical temperature T max c ≈ 0.09 close to the QCP. We find that the value for T max c is in good agreement with an independent self-consistent Eliashberg calculation which gives T El c ≈ 0.08. The nature of the superconducting state reveals itself in the uniform pair- Note that under a π 2 rotation ∆ η → η∆ † η . Hence, ∆ − has d-wave character while ∆ + is of s-wave nature. In Fig. 6 we show the momentum resolved pairing susceptibilities P η across the first Brillouin zone. The strong peak in the sign-changing symmetry channel is evidence for a d-wave superconducting state.
In principle, quantum critical fluctuations can also promote other electronic orders like CDW formation [16,55,56]. For particle-hole symmetric electron dispersions, d-wave superconductivity and d-wave CDW are degenerate due to an emergent SU (2) symmetry [57]. Inspecting the relevant susceptibilities (see App. C), we find that CDW correlations are mildly enhanced at low temperatures, most strongly in the vicinity of the QCP, but are only weakly dependent on system-size. Interestingly, this amplification seems to be irrespective of the onset of superconductivity at T = T c , indicating the absence of a competition between CDW and SC.

B. Quantum critical correlations
To extract the critical properties of model (1) we suppress the superconducting transition by studying the system for smaller Yukawa coupling λ = 1 where the critical temperature T c < 1 40 is below the lowest accessible temperatures. Simulations are then performed in the vicinity of the QCP and above any possible ultralow-temperature superconducting phase.

Magnetic correlations
We investigate the system's magnetic correlations by analyzing the bosonic SDW susceptibility at momenta q, taken relative to the ordering wave vector Q, and Matsubara frequencies ω n = 2πnT . First, the dependence of the inverse susceptibility χ −1 (q, iω n ) on Matsubara frequency in the vicinity of the QCP is illustrated in Fig. 7. The data very visibly has linear character for small Matsubara frequencies ω n , both for q = 0 and small finite momenta q > 0, with an apparent cusp at ω n = 0. To establish the presence of a ω n -term in χ −1 we perform a linear regression for the small frequency data. The resulting fits to the form a 1 ω n + a 0 , displayed in Fig. 7 as black lines, are in good agreement with the data and confirm the linear Matsubara frequency dependence. This finding suggests overdamped dynamics of the boson field φ due to interactions with the fermions. At finite Matsubara frequencies, finite-size effects are negligibly small (within error bars), as evident in the data collapse of χ −1 (q = 0, iω n ) for different system sizes, 8 ≤ L ≤ 14.
Turning to the momentum dependence of χ −1 (q, iω n ) next, we find that the momentum dependence shown in Fig. 8 is consistent with a quadratic form q 2 for small momenta q. This holds both for ω n = 0 and small non-vanishing Matsubara frequencies ω n . Similar to above we establish the presence of a q 2 term in χ −1 (q, iω n ) by fitting the data to the form a 1 q 2 + a 0 . The results are indicated in Fig. 8 as black lines, showing good agreement with the DQMC data. In combination with the observed linear Matsubara frequency dependence, this provides strong evidence for a dynamical critical exponent z = 2. Next, we illustrate the dependence of the inverse susceptibility χ −1 (q = 0, iω n = 0) on the tuning parameter r in Fig. 9. For tuning parameter values r ≥ r 0 ≈ −1.89 we find that the data for different system sizes follows a linear dependence. Due to finite-size and finite-temperature effects the onset of a finite value of χ −1 does not appear to be as abrupt but instead is moderately continuous [37].
Finally, we inspect the temperature dependence χ −1 (T ) close to the critical point in Fig. 10. Within the limits of our numerical accuracy, we find that the DQMC results are consistent with a T 2 term at low temperatures. We highlight this point by fitting the data to a second degree polynomial (black line) which is able to adequately capture the temperature trend over a broad range 0.025 < T ≲ 0.5. At higher temperatures, the situation changes and we find a linear T -dependence, as shown in the inset of Fig. 10.
Taken together, these findings indicate that the inverse SDW susceptibility near the metallic QCP (at low temperatures) has the form, Remarkably, this is precisely the functional dependence predicted by Hertz-Millis theory, in which fermion degrees of freedom are integrated out -despite their gapless natureto obtain an effective critical theory for φ [24,25]. However, our DQMC simulations also suggest certain deviations from Hertz-Millis theory, specifically in terms of the observed quadratic temperature dependence (in lieu of a linear scaling). Interestingly, this general comparison to Hertz-Millis theory is in full analogy to what has been observed for an O(2) order parameter in Refs. [17,37].
Let us close this section by noting that it is always a possibility that the character of the onset of magnetism is first order rather than continuous. When we turn off the quartic boson coupling in model (1), i.e. u = 0, we indeed observe a region of coexistence of magnetic and non-magnetic states. This is indicated by a jump-like tuning parameter dependence (App. A) of the magnetic susceptibility χ −1 (q = 0, iω n = 0) and an apparent double-peak structure in its distribution function. Keeping u finite and positive should, however, generally favor a continuous transition.

Single fermion correlations
Upon approaching the QCP from the disordered phase, the fermions at the hot spots lose their coherence. This is manifested in a substantial growth of the imaginary part of the (Matsubara) self-energy, defined as Σ k (ω n ) = iω n − k −  = (k f , 0). The data points are averaged over different twisted boundary conditions, while the width of the line is the standard deviation of ImΣ between the different boundary conditions, thus representing the uncertainty due to finite-size effects, see text. Fig. 11. At the QCP, r ≈ r c , the selfenergy at the hotspots is very weakly dependent on both temperature and frequency. We identify the lowest temperature studied here (T = 0.05) as the non-Fermi-liquid crossover temperature, at which the self-energy begins to dominate the bare frequency dependence of the Green's function, namely ImΣ k (ω n = πT ) = πT . For r < r c the imaginary part of the self-energy tends to diverge, indicating the gapping out of the hot spots. The self energy away from the hotspots remains small, and tends to vanish at low frequencies.
On a technical note, we mention that the leading cause of uncertainty in the numerical estimation of the self-energy are finite size effects. To minimize these effects we used data from several sets of twisted boundary conditions, as explained in more detail in Appendix B. The time discretization error in G k (ω n ) is reduced by using the 'Filon-Trapezoidal' rule [58].

V. LOCAL NESTING
The recent non-perturbative calculations by Schlief et al. [26] revealed a z = 1 fixed point with characteristics very much different from our findings above. As an attempt to promote a flow to this self-consistent critical theory, we modify our original spin-fermion model to host a Fermi surface with (almost) perfect local nesting near the hot spots. Some of the low-energy properties of this system close to criticality are then extracted in ultra-low temperature DQMC simulations.
Specifically, we introduce higher order hopping processes into Eq. (1) and choose t xh = 1 = t yv , t ′′ xv = 0.45 = −t ′′ yh , where t ′′ xh and t ′′ yv are second neighbor hopping amplitudes for ψ x in horizontal and ψ y in vertical direction, respectively (Fig. 2). All other hopping amplitudes, as well as the quartic boson coupling u are set to zero, whereas µ x = −µ y = −0.47. The resulting Fermi surface is illustrated in Fig. 12. Note that the relative angle between the Fermi velocities for this set of parameters is θ ≈ 8.5 ○ whereas the Fermi surface in Fig. 1 has Fermi surface close to local nesting near the hot spots across the Brillouin zone. The black lines correspond to the two fermion flavors ψx (solid), ψy (dashed). One band has been shifted by Q = (π, π) such that hot spot pairs (red) occur at crossing points. θ ≈ 36.9 ○ (the magnitude of v F does not change significantly).
We start off by discussing a potential superconducting transition. In an almost locally nested system antiferromagnetic fluctuations will only slightly enhance T c [5,18,26]. As indicated in Eq. (7), the energy scale associated with superconductivity is expected to decrease as a function of the angle, T c ∼ sin θ. This relation has been numerically confirmed for a generic Fermi surface coupled to an O(2) order parameter in Refs. [18,37]. Anticipating that superconductivity arising from fermions at the hotspots, Eq. 7, is independent of the precise band structure [18] we expect a fourfold reduction of T c for the Fermi surface shown in Fig. 12 compared to the one considered in Sec. II. Indeed, we do not observe superconductivity in our DQMC simulations for temperatures as low as T c = 0.01.
In Fig. 13, we show the Matsubara frequency dependence of the inverse SDW susceptibility χ −1 (q = 0, iω n ). In the interval 0 < ω n < 1 a distinct curvature can be seen, which is more pronounced at smaller Matsubara frequencies. To resolve this nonlinearity, simulations have been performed at ultra-low temperatures, β = 100. For ω n > 0, data points obtained from DQMC simulations of differently sized systems fall on top of each other, indicating the absence of significant finite size effects. To understand the origin of the curvature, we consider the non-interacting fermionic susceptibility Π 0 . As illustrated in Fig. 13, Π 0 (ω n ) qualitatively shows the same trend as χ −1 (ω n ) for small Matsubara frequencies. This suggests that the observed nonlinearity in the SDW susceptibility is likely due to low-energy features of the band structure.
Finally, we show the dependence of the inverse SDW susceptibility on squared momentum in Fig. 14. For finite q > 0 the DQMC results for χ −1 (q, iω n = 0) are consistent with a q 2 term, while a noticeable drop is visible at q = 0. We perform a linear regression to establish the quadratic momentum dependence, illustrated in Fig. 14, and find good agreement with the numerical data over the range 0 < q 2 ≤ 2. Importantly, all finite-q data points collapse onto a single line and do not branch out as a function of squared momentum. The isotropic behavior of the boson propagator lies in stark contrast to an anisotropic q x + q y + q x − q y term predicted at the z = 1 fixed point [26]. However, the apparant discontinuity at q = 0 indicates there are features at low momenta that cannot be resolved for the available system sizes. In particular, the z = 1 anistropic behavior may emerge at smaller values of q. We note that the momentum dependence is similar at higher temperature (T = 0.025, see inset to Fig. 14).

VI. DISCUSSION
In this work, we have deployed large-scale DQMC simulations to explore the phase diagram of itinerant fermions coupled to an isotropic antiferromagnetic order parameter in an unbiased and rigorous manner. By inspecting fermionic and bosonic correlations in the vicinity of the QCP, which marks the onset of SDW order, we were able to unveil the critical low-energy behavior up to numerical accuracy.
Our main finding is that, over a broad range of temperatures and the tuning parameter near the QCP, the critical SDW fluctuations are remarkably well described by Hertz-Millis theory, in which one naively integrates out the fermions [24,25]. The dynamical exponent in this regime, within our accuracy, is z = 2. This is surprising in view of the fact that Hertz-Millis theory neglects infinitely many marginal couplings that are generated when integrating over the gapless fermion degrees of freedom at the Fermi surface [23]. Qualitatively similar observations have been made for metals with easy-plane XY [37] and Z 2 antiferromagnetic order [39], despite noticeable differences in the perturbative structure [15], which hints towards a generic property of the SDW QCP.
A further important aspect of our DQMC results is that, away from perfect local nesting at the hotspots, the QCP always appears to be shadowed by high-temperature superconductivity. We find that the superconducting order parameter is unambiguously of d-wave character, i.e. changes sign under π 2 rotations. The observed maximal critical temperature T max c ≈ E F 20 for λ = 2 is in good agreement with Eliashberg theory. Anticipating that the energy scale associated with superconductivity due to fermions at the hotspots, Eq. 7, is independent of the band structure parameters [18] our T c can be compared to Ref. [17,37], which considers a very similar model with easy-plane AFM fluctuations. We find T 2, consistent with a quadratic dependence of T c on the number of order parameter components.
In much of the "quantum critical" regime above the superconducting T c , the fermions at the hot spots remain underdamped, in the sense that their self-energy is smaller than their energy, except at the lowest temperatures approaching T c . Above the QCP, the imaginary part of the self-energy at the hot spots is only weakly dependent on Matsubara frequency, in stark contrast to the behavior expected within Fermi liquid theory. On the ordered side of the QCP the self-energy seems to be diverging, which is consistent with a gapping out of the hot spots.
Motivated by the recent prediction of a novel z = 1 fixed point [26][27][28] by Sung-Sik Lee and collaborators, we consid-ered a variant of the model where the hot spots are (nearly) locally nested. In this case, the superconducting T c is strongly suppressed below the lowest temperature available in our study, T ≈ E F 200. Above this temperature, we find substantial deviations from the z = 2 Hertz-Millis behavior, but no evidence for the predicted z = 1 criticality. Notably, the momentum dependence of the SDW order parameter is isotropic up to numerical accuracy.
A further investigation of the possibility of a z = 1 fixed point manifestation at even lower temperatures and with different band structures is desirable. Moreover, directly extracting the renormalized vertex function from DQMC simulations and analyzing its Matsubara frequency dependence seems like a viable approach to shed important light onto why Hertz-Millis theory is capable of capturing large parts of the critical physics correctly, despite being formally uncontrolled. On a different front, the recent advances in the field of quantum machine learning [59][60][61][62][63] provide a new exciting route to accessing transport properties [64] and might provide a computationally efficient way to identify and map out an extended non-Fermi liquid regime.
To reduce finite size effects in our simulations we thread a single magnetic flux quantum φ 0 through our system by introducing Peierls phases into the fermion kinetic energy  [17,37,66], Here, A αs ij = 2π φ0 ∫ j i A αs ⋅ dl and we choose a vector potential in Landau gauge, A αs (r) = −B αs yx. To retain the antiunitary symmetry, which renders our model sign-problem free, we choose flavor and spin dependent amplitudes, Note that this choice explicitly breaks the SU(2) spin rotation symmetry of the model. However, the symmetry breaking field scales as L −2 , rendering it unimportant as long as the antiferromagnetic correlation length is ξ ≲ L. Numerically, we find that the parallel and perpendicular components of the magnetic susceptibility agree (up to error bars) for r > r QCP and only start to diverge at low temperatures on the ordered side of the QCP.
Explicitly, we implement the following phases for nearest, next-nearest, and next-next nearest neighbor hoppings (Fig. 2 in the main text).
Nearest neighbor (B3)  To study the fermionic Green's function in momentum space, we cannot apply such a flux, but rather use twisted boundary conditions. Here, too, to avoid the sign problem, the twist angles ϕ have an orbital and spin dependence ϕ x↑ = ϕ y↓ = −ϕ x↓ = −ϕ y↑ .

Next-nearest neighbor
(B6) We use ϕ = 2πn 4 for n = 1...4 to increase the set of available momenta 16 fold. To study the self-energy on the Fermi surface, as in Fig. 11 we combine the data from the different boundary conditions, as mentioned above, and identify the local maxima of the G k (τ = β 2) at T = 0.2 as the Fermi surface. The hot spots are the closest point to the intersection of the Fermi surfaces of the two bands ψ x and ψ y . We take the self energy to be the average of the self-energies at the hotspot and the 4 adjacent momenta (which correspond to different boundary conditions). The error is taken to be the standard deviation across the aforementioned momenta. We follow a similar procedure for Σk x .
Appendix C: Charge-density correlations In order to investigate the possibility of the presence of charge-density wave (CDW) fluctuations we inspect the susceptibilities C ± (r) = dτ ⟨∆ † ± (r i , τ )∆ ± (0, 0)⟩, In Fig. 17 we show the temperature dependence of the maximum of C ± (q), omitting q = 0, for three different tuning parameter cuts through the phase diagrams shown in Fig. 5. Both C + and C − are peaked close to (but not quite at) the corners of the Brillouin zone at momentum q = (π, q max ≈ 0.83π). First, we notice that C + < C − consistently for all chosen parameters. Second, we observe that while C + is suppressed with decreasing temperature for r ≪ r c , indicating a competition with SDW fluctuations, C − is enhanced for all tuning parameter values. Focusing on C − (T ), we observe that the amplification is most pronounced close to the QCP, r ≳ r c , and seems insensitive to the onset of superconductivity at T = T c . This behavior is different from that of an easy-plane antiferromagnetic QCP, where the C − was found to drop sharply below T c [17]. Despite this mild enhancement of CDW correlations at low temperatures in the O(3), we find that their system size dependence across the phase diagram is very weak (not shown), indicating that there is no CDW long-range order.