Non-Hertz-Millis scaling of the antiferromagnetic quantum critical metal via scalable Hybrid Monte Carlo

A key component of the phase diagram of many iron-based superconductors and electron-doped cuprates is believed to be a quantum critical point (QCP), delineating the onset of antiferromagnetic spin-density wave order in a quasi-two-dimensional metal. The universality class of this QCP is believed to play a fundamental role in the description of the proximate non-Fermi liquid behavior and superconducting phase. A minimal model for this transition is the O(3) spin-fermion model. Despite many efforts, a definitive characterization of its universal properties is still lacking. Here, we numerically study the O(3) spin-fermion model and extract the scaling exponents and functional form of the static and zero-momentum dynamical spin susceptibility. We do this using a Hybrid Monte Carlo (HMC) algorithm with a novel auto-tuning procedure, which allows us to study unprecedentedly large systems of 80 × 80 sites. We find a strong violation of the Hertz-Millis form, contrary to all previous numerical results. Furthermore, the form that we do observe provides good evidence that the universal scaling is actually governed by the analytically tractable fixed point discovered near perfect “hot-spot’" nesting, even for a larger nesting window. Our predictions can be directly tested with neutron scattering. Additionally, the HMC method we introduce is generic and can be used to study other fermionic models of quantum criticality, where there is a strong need to simulate large systems.

Once r c (T ) is obtained for several T values, we extrapolate them to T = 0 using a simple polynomial fit, r c (T ) = r c − a T b , which gives us our estimate of the quantum critical point r c ≡ r c (T = 0). The values of r c that we obtain for v i , i = 1, . . . , 5 are all given by r c ≈ 1.2554. Putting everything together gives the phase diagram in Fig. 4. The error bars denote the one sigma confidence interval.
To probe for superconductivity we use the method of Ref. [1]. We observe no superconductivity at r c down to the lowest temperature that we study, T = 1/80. We give more details on this computation in Supplementary Note 2.
In Supplementary Fig. 2 we show plots of the occupation function for more nesting parameter values, namely v = v 5 .

B. Scaling of spin susceptibility at criticality
Having established the phase diagram and located the QCP, we examine the boson correlation function at criticality. We set r = 1.26 ≳ r c , which is very slightly further into the disordered phase than the extrapolated value r c . This is done to stay on the safe side of the transition, i.e. within the critical fan.
The log-log plots of χ −1 (ω) − χ −1 (0) for β = 80, L = 20 and of χ −1 (q x ) − χ −1 (0) for β = 20, L = 80 are shown in Fig. 5a and Fig. 5b, respectively. As noted there, the regions of intermediate algebraic scaling are chosen by eye. For all nesting parameter values we define the regions as 10π/β ≤ ω < 2.0 and 8π/L ≤ q x < 2.3. The slope of the linear fits in these regions determine the critical exponents z and η ϕ , shown in Fig. 6. To compute the error bars, we used a bootstrap method, where the selected data points were re-sampled 10000 times, allowing for duplicates.
The symmetry of the full spatial dependence is seen from Fig. 7 for v = v 5 . In Supplementary Fig. 3 we provide similar zoomed-in plots for the other v i , showing the C 4 symmetry at small momenta. In Supplementary Fig. 4, we also show that the contour at the same small momenta only get more 'square-like' with increasing L, indicating that the emergent C 4 symmetry is not a finite-size effect but a true long-wavelength feature of the theory in the thermodynamic limit.

SUPPLEMENTARY NOTE 2. SUPERFLUID DENSITY
In order to test for superconductivity in our system, we apply the method of Ref. [1] to compute the superfluid density ρ s . The system is a superconductor if ρ s surpasses the universal BKT value of ∆ρ s = 2T /π.
The current density operator in thex direction is given by The static current-current correlator in momentum space is defined as where we understand r = (x, y) as usual.
The superfluid density in the thermodynamic limit is then given by For a finite L, we use depends weakly on L, as is seen by others [2]. For all nesting parameters, at the lowest temperature (β = 80), we measure ρ In order to recover the Green's function as an expected diagonal of the form of Eq. (36) we simply take where F is the twisted discrete Fourier transform defined in the Methods Section. Note that by rewriting where G 0 ϕ := D(ϕ) −1 , and the diagonal matrices A and B are defined by and Observe that B = B q implicitly depends on the choice q, but for any specific choice of q, the angle-bracketed expression in Supplementary Eq. (7) can be obtained in computational time dominated by the cost of a constant number of linear solves of the form of Eq. (22). Note that for any term in which diagonals appear within a product, independent diagonal estimators of the form of Eq. (38) must be used for each in order to obtain a consistent estimator for the expression in Supplementary Eq. (7).

SUPPLEMENTARY NOTE 4. ORDER OF THE TRANSITION
Although the effective action of Eq. (1) is modeling a continuous phase transition, the transition can be first-order, especially for u = 0 [3]. In order to check that the transition is indeed second-order for all the parameter values we use in the main text, we plot the inverse spin susceptibility χ −1 = χ −1 (0, 0) as a function of r for v = v 1 , β = 40, L = 20. We can see that even at low temperatures, the susceptibility turns on in a continuous fashion. This is in contrast to The error bars denote the one sigma confidence interval.
the first-order behavior observed in Ref. [4] for certain parameter values, where χ −1 develops a pronounced jump at lower temperatures.

SUPPLEMENTARY NOTE 5. PRECONDITIONER PERFORMANCE
In this section we numerically validate our choice of preconditioner as described in the Methods Section. All experiments are performed for our standard choices of model parameters u = 0.0, g = 0.7 √ 2, c = 3.0, as well as the choice v = v 2 ≈ 0.072 for the nesting parameter and the critical parameter value r = 1.26.
In Supplementary Fig. 6, we fix N τ = 100 and plot the average speedup provided by the preconditioner over an HMC run, measured both in terms of wall clock time and CG iterations. To bolster the practical conclusions of Supplementary Fig. 6 and see more concretely how the preconditioner affects the conditioning of the linear solve, we plot in Supplementary Fig. 7  Finally, in Supplementary Fig. 8 we validate the claim from Methods Section that the advantage of the preconditioner becomes more significant in the limit of small Trotter step ∆ τ . Although in the rest of this work we only consider ∆ τ = 0.1, this limit nonetheless provides one way of understanding the advantage of the preconditioner besides the motivation via exactness for a the choice of zero coupling g = 0. Here we provide more details on the numerical performance of our algorithm, expounding on the summary main in the the Methods Section. As a concrete observable, we track the growth of the integrated autocorrelation time τ int of the total SDW susceptibility χ ≡ χ(0, 0) at criticality. To test the performance, we ran the HMC algorithm in two studies: one in which the lattice volume is fixed to sizes L = 20, 30 and inverse temperature scaled from N τ = 100 to N τ = 300, and one in which the number of imaginary time steps is fixed to N τ = 100, 150 while lattice extent L is scaled from 20 to 60. For each scenario, empirical averages and error bars were computed using 5 Markov chains. In Supplementary Fig. 9, we see that our algorithm exhibits constant scaling of τ int with respect to both lattice volume and inverse temperature.
The number of leapfrog integration steps per effective sample is plotted in Supplementary Fig. 10. This number roughly tracks the total algorithmic cost in terms of the number of linear solves of the form Eq. (22). A power law fit y = ax b is applied to the measurements to extract the scaling exponents z 1 , z 2 . With respect to N τ , we observe that z 1 = b − 1/4 ≈ 0.5. With respect to V , we observe that z 2 = b − 1/4 ≈ 0. As noted in the main text, this implies an absence of critical slowing down with respect to the lattice volume V for the auto-tuned HMC algorithm presented in this work.
Finally, in Supplementary Fig. 11, we show the wall clock time per effective sample with respect to V and N τ , where we can see that the scaling is roughly linear with respect to V and superlinear with respect to N τ .