Ising nematic quantum critical point in a metal: a Monte Carlo study

The Ising nematic quantum critical point (QCP) associated with the zero temperature transition from a symmetric to a nematic {\it metal} is an exemplar of metallic quantum criticality. We have carried out a minus sign-free quantum Monte Carlo study of this QCP for a two dimensional lattice model with sizes up to $24\times 24$ sites. The system remains non-superconducting down to the lowest accessible temperatures. The results exhibit critical scaling behavior over the accessible ranges of temperature, (imaginary) time, and distance. This scaling behavior has remarkable similarities with recently measured properties of the Fe-based superconductors proximate to their putative nematic QCP.


I. INTRODUCTION
A hallmark of strongly correlated electron systems is the competition of ground states with different kinds of order. [1] In this context, a central set of unsettled theoretical issues concerns the character of quantum critical points (QCPs) in metals [2][3][4]. Such metallic QCPs have been identified in several heavy fermion compounds [5]; evidence for quantum critical behavior has also been found in the ruthenate Sr 3 Ru 2 O 7 [6,7], and the cuprate and iron-based superconductors [8]. Quantum criticality is also often invoked as a possible explanation of the "strange" or "bad" metal behavior seen in a variety of such materials [9][10][11][12][13][14][15][16].
In this paper, we report a DQMC study of a twodimensional sign problem-free lattice model that exhibits an "Ising nematic" QCP in a metal at finite fermion density; in the nematic phase, the discrete lattice rotational symmetry is spontaneously broken from C 4 to * These authors have contributed equally to this work. C 2 . This is a particularly relevant QCP given that nematic order [56][57][58][59][60][61] and nematic quantum critical fluctuations [62][63][64][65][66][67] have been observed in many of the materials mentioned above.
Our simulations are limited to finite system sizes -up to 24 × 24 lattice sites. Within our numerical accuracy, we find evidence for a continuous nematic quantum phase transition with critical behavior that differs significantly from that of a nematic QCP in an insulator. Specifically, in the disordered phase near the QCP, where the dimensionless quantum control parameter h ≥ h c , the thermodynamic (zero frequency) nematic correlation function (defined in Eq. 4) is consistent with the following func- tional description: where T is the temperature, q is the wave-vector, ω n is the Matsubara frequency, and A, b, and κ are positive dimensionful constants. We find that the exponents in this expression take values λ = 1.0±0.1 and γ = 1.0±0.1. (See Fig. 8.) These values are in agreement with the predictions of Refs. [3,16,68,69] for the behavior of the thermodynamic susceptibility near criticality.
This implies that the uniform nematic susceptibility χ(h, T ) ≡ D(h, T, 0, 0) has a Curie-Weiss form with an effective Weiss temperature which varies linearly with (h c −h), as can be seen in Figs. 1 & 5. Eq. 1 is consistent with the scaling relation, D(h, T, q, 0) = ξ γ/ν Φ (x, y), where ξ ∼ |h − h c | −ν , x = T ξz, y = |q|ξ. Here, γ is the conventionally defined susceptibility exponent, the correlation length exponent ν ≡ γ/(2 − η), where η is the anomalous dimension of the nematic field, and we have introduced an "apparent dynamical exponent," z = (νλ) −1 . The values of these exponents derived from our Monte Carlo data are given in Table I. We also include for comparison the exponents of the standard 2dimensional Ising QCP, which apply when the coupling to fermions is set to zero.
Despite the accuracy of the scaling analysis, it is important to test whether the observed behavior is characteristic of an asymptotic quantum critical scaling regime. To this end, we define a "quadrupolar" correlator, Q, in terms of fermion bilinears, which has the same symmetry as the nematic correlator, D, and therefore should have the same asymptotic critical behavior. As can be seen in Figs. 9, 10, 11, 14, Q does not obey a scaling law of the form of Eq. (1) over a comparable range of (h, T, q) as does D. Sufficiently close to criticality, Q may be consistent with Eq. (1), but the dynamic range is relatively small and the error bars on the apparent exponents are substantial. It could be that corrections to scaling are smaller in D than in Q, or it could be that our data simply do not reflect asymptotic quantum critical scaling.
In the fermionic sector we find results consistent with strongly renormalized Fermi liquid, "marginal Fermi liquid" [70], or weakly non-Fermi liquid behavior down to our lowest temperatures (T min ≈ 0.02E F , where E F is the Fermi energy). In particular, at the QCP, the effective quasiparticle weight, Z k F (T ), defined in terms of the single-fermion Green's function in Eq. (15), remains substantial. However, it monotonically decreases on cooling, with downward curvature. In a Fermi liquid, Z k F (T ) would approach a positive limit as T → 0, while in a weak non-Fermi liquid, it would vanish in proportion to a small power of T . Fig. 2 shows maps of the low-energy spectral weight as a function of momentum in the disordered phase, at the QCP, and in the ordered phase. The existence [71] of "cold-spots" along the zone diagonals, where the quasiparticles are relatively weakly scattered, is apparent. Finally, no superconducting transition Fermion Green's function, G(k, τ = β/2), as a function of momentum k across the Brillouin zone for three different values of the transverse field: h = 3.35 (in the nematic phase), h = 3.525 (near the QCP), and h = 3.7 (in the disordered phase). In the nematically ordered phase, the data is averaged over both orientations of the order parameter. G(k, τ = β/2) is proportional to the integral of the spectral function A(k, ω) over an energy window of width T [see Eq. (12)]. The temperature is T = 1/8, and the system size is L = 20. The data were taken from systems with either periodic or anti-periodic boundary conditions. The Fermi surfaces are seen as peaks in G(k, τ = β/2). The lower right panel shows the Fermi surface of the bare band structure.
is found down to our lowest temperatures, although the superconducting susceptibility in the s-wave channel is peaked about h ≈ h c . (See Fig. 16.) The remainder of the paper is organized as follows: in Sec. I we define the lattice Hamiltonian and describe its phase diagram; in the next section, we provide evidence for the statements above regarding thermodynamic nematic and quadrupolar correlations (III A), dynamical nematic and quadrupolar correlations (III B), superconductivity (III C), and single-fermion correlations (III D); finally, in Sec. IV, we discuss various caveats concerning the interpretation of our results, and their bearing on both prior theoretical work and the interpretation of experiments.

II. MODEL AND PHASE DIAGRAM
Our model is illustrated in Fig. 3, and is defined on a two-dimensional square lattice. Every lattice site hosts a single (spinful) fermionic degree of freedom. Each link has a pseudospin-1/2 degree of freedom that couples to the fermion bond-density. The system is described by the following Hamiltonian: where Here, c † jσ creates a fermion on site j with spin σ =↑, ↓, i, j denotes a pair of nearest-neighbor sites on the square lattice, t and µ are the hopping strength and chemical potential, respectively, τ a i,j (a = x, y, z) denote Pauli matrices that act on the pseudospin that lives on the bond connecting the neighboring sites i and j, V > 0 is the nearest-neighbor Ising interaction between neighboring pseudospins (here, i, j ; k, l denotes a pair of nearest-neighbor bonds), h is the strength of a transverse field that acts on the pseudospins, and α is the dimensionless coupling strength between the pseudospin and the fermion bond density.
The pseudospins are not related to the physical spins of the electrons; their ordering corresponds to a nematic Illustration of the lattice model (2). Spinful electrons reside on the sites, while the Ising pseudospins live on the bonds. The pseudospins interact with their neighbors antiferromagnetically, and are coupled to the fermion bond density. transition. The model (2) should be viewed as an effective lattice model, designed to give a nematic QCP. Microscopically, the nematic degrees of freedom could represent (via Hubbard-Stratonovich transformation) interactions involving the same electrons that form the Fermi surface, or another, independent degree of freedom (such as a phonon mode that becomes soft at a structural transition). So long as the properties of the QCP are universal, the low-energy behavior does not depend on its microscopic origin.
For α = 0, the system is composed of two decoupled sets of degrees of freedom: free fermions, and pseudospins governed by H b , which has the form of a d = 2 transverse field Ising model. At zero temperature, the pseudospins undergo a second-order quantum phase transition from a paramagnet to an "antiferromagnet" that breaks 90 • rotational symmetry at h = h c0 . This transition is in the 3 dimensional classical Ising universality class. At T = 0, the pseudospin degrees of freedom are gapped in both the nematic and isotropic phases. At finite temperatures, a line of second-order classical d = 2 Ising transitions extends from the QCP in the h − T plane.
For non-zero α, the phase diagram is similar, but exhibits quantitative and qualitative modifications. Fig. 1 shows the phase diagram, obtained by DQMC, for both α = 0 and α = 0. The α = 0 transition between the nematic and isotropic phases remains second order, and extrapolates to a new QCP, shifted relative to h c0 . More striking is the change in the slope with which the phase boundary, T N (h), approaches the QCP. For α = 0, the slope diverges at low temperature, consistent with the expectation for the transverse field Ising transition, where T N ∝ |h − h c0 | νz with ν = 0.63 and z = 1. In contrast, for α > 0, we find that T N ∝ (h c − h). On the disordered side of the transition, we define a crossover line by identifying h cross (T ) as the value of h at which the nematic susceptibility at fixed T has fallen to half of its value at h = h c . T cross (h), the inverse of h cross , also vanishes linearly with h upon approaching the QCP, although its slope is steeper than that of T N (h).
The linear behavior of T N (h) for small h c − h is also seen for other model parameters. In Fig. 4 we show the phase diagram for two values of the fermion density, controlled by the chemical potential µ. As the fermion density is reduced, both h c − h c0 and the range over which T N (h) is linear become smaller, indicating that the effect of the coupling between electrons near the Fermi surface and the nematic modes becomes weaker. The fact that T N (h) appears linear at low temperature for both values of fermion density is consistent with this being a universal property of the metallic QCP.

III. DQMC RESULTS AND ANALYSIS
We have simulated the Hamiltonian (2) with system sizes between L = 8 and L = 24, and temperatures between 0.025t and 5t. These simulations do not suffer   from the minus sign problem: the fermion determinants obtained from integrating out the spin up and spin down fermions are identical and real, so the product is nonnegative. Global updates of the pseudospin space-time configurations are introduced in order to overcome critical slowing down.
Some technical details related to the DQMC simulations are discussed in Appendix A. Here, we shall focus on the results and data analysis. Unless stated otherwise, the model parameters used in the simulations were: V = t, α = 0.5, µ = −0.5t, t = 1.

A. Thermodynamic correlations
The two-point function for the nematic order parameter field is defined as   Here, r ij = r i − r j and the nematic order parameter is defined as N i = j η ij τ z ij , with η ij = 1/4 for r ij = ±x, η ij = −1/4 for r ij = ±ŷ, and η ij = 0 otherwise. The phase diagram discussed above was determined using the nematic susceptibility and standard finite size scaling techniques, described in Appendix B. We also define a quadrupolar order parameter, , the susceptibility is seen to follow a power law, χ(h = h c , T ) = A/T λ with λ ≈ 1, over more than a decade of temperature between 0.1V < T 2V . The slope of log(χ) vs. log(T ) is independent of the system size for L ≥ 12.
To establish the scaling behavior of χ away from h c , we plot the susceptibility for different temperatures and system sizes, as a function of h − h c [ Fig. 6(a)]. For the system sizes displayed (16 ≤ L ≤ 24), χ is only weakly dependent on L. In Fig. 6(b) we plot T λγ χ(h, T ) for h > h c as a function of (h − h c )/T λ , with λ = 1.0 and γ = 1.0. The different curves collapse onto each other, as expected for a scaling function.
Next, we examine the momentum dependence of the static correlator D(h, T, q, iω n = 0) near the QCP. Fig. 7 shows D −1 (h, T, q, iω n = 0) as a function of |q| 2 for various values of h − h c and T . For small momentum, the momentum dependence consists of an essentially isotropic term, approximately proportional to |q| 2 . Combined with the h and T dependences discussed above, we deduce that the static nematic correlator is In each case the momentum dependence remains quadratic. Momenta from multiple system sizes (represented by squares for L = 16, circles for L = 20, and diamonds for L = 24) fall on the same curve, indicating that the properties shown are in the thermodynamic limit. well-described by Eq. (??), with the aforementioned exponents, and with A ≈ 2.1, b ≈ 2.9t, κ ≈ 2.6t, and the lattice spacing set to unity. This is explicitly demonstrated in Fig. 8. (In the figure, A(h, T, q, ω n ) is the functional approximant to D(h, T, q, ω n ) , defined in Eq. 6.) In order to evaluate the consistency of the scaling form, Eq. (1), we have also measured the quadrupolar correlator, Q, over a similar range of parameters as D. In Fig.  9, we compare the nematic and the quadrupolar suscep- FIG. 9. Nematic and quadrupolar susceptibilities vs inverse temperature at h = hc. Multiple system sizes are shown with squares for L = 16, circles for L = 20, and diamonds for L = 24. The error bars reflect uncertainties due to finite system sizes, estimated from the differences between the measured susceptibilities of systems with different boundary conditions. 10. Inverse nematic and quadrupolar susceptibilities vs h − hc at various temperatures. Multiple system sizes are shown with squares for L = 16, circles for L = 20, and diamonds for L = 24. The error bars reflect uncertainties due to finite system sizes, estimated from the differences between the measured susceptibilities of systems with different boundary conditions.  tibilities at h = h c as a function of temperature. Like D, Q varies approximately linearly as a function of 1/T , although there are large error bars at low temperatures. In Fig. 10 we show the h dependence of D −1 and Q −1 at various temperatures. In both cases, the h dependence is approximately linear for small h − h c . However, the range of h − h c in which D −1 is linear is larger than the corresponding range of linearity of Q −1 , and the slope of The momentum dependence of the two correlators is shown in Fig. 11. Both appear isotropic, depending only on |q| 2 . However, Q −1 , in contrast to D −1 , has noticeable downward curvature. D −1 depends on temperature through an essentially momentum-independent shift, while the temperature dependence of Q −1 is more complicated. For more details see Appendix E.

B. Dynamic correlations
The dependence of D on Matsubara frequency, ω n = 2πT n, is shown in Fig. 12. At intermediate ω n (0.5t ω n 2t) and for q = 0, D −1 is an approximately linear function of |ω n | with a slope that is independent of T , and also independent of both the direction and magnitude of q. However, at the smallest non-zero momenta, a different frequency dependence is visible for |ω n | < 0.5t. This is emphasized in the inset of  where D −1 (q, iω n = 0) has been subtracted. The frequency dependence at q = 0 is different from the behavior seen for non-zero q; for ω n t, D −1 (q = 0, iω n ) appears to saturate to a frequency independent value, followed by a sudden drop at ω n = 0. The frequency dependence of D −1 away from the critical point, shown in Fig. 23 in the Appendix E, is similar to the behavior of D −1 at criticality, apart from a constant positive shift (that corresponds to the fact that we are on the disordered side of the transition, and hence D(q = 0, iω = 0) remains finite in the thermodynamic limit). The bulk of the data at small values of h − h c , T , q, and ω n is well described by the simple "functional approximant" As already shown in Fig. 8, this approximant is extremely accurate for all the data at ω n = 0. As can be seen from Fig. 24 in Appendix E, it does an equally good job for the dynamical response for q not too small, but fails to describe the observed dynamics for the few lowest values of |q|.
The equal time correlator gives complementary information about the dynamics, since it integrates over all frequencies. Moreover, "equal time" is the same in both real and imaginary time. Thus, in Fig. 13 we exhibit the inverse of the equal time correlator,D(τ = 0, r) for h = h c at T = 0.1t, as a function of distance for various values of L. Note that for the largest size system (L = 24),D(0, r) ∼ 1/|r| 2 , consistent with what one would get by Fourier transform of A.D is shown for all orientations of r, so the narrow spread of points shows the high degree of isotropy exhibited by the data. Note also that for relatively small |r|, the data is independent of L, at least for L ≥ 16, but that at longer distances the data for the smaller systems deviate systematically from the |r| −2 behavior.
The frequency dependence of Q is shown in Fig. 14. It is clearly qualitatively different from the behavior of D in Fig. 12, and correspondingly cannot be approximated by any expression of the form in Eq. 6. The uniform quadrupolar density is nearly conserved, as is apparent from the very large values of Q −1 (q = 0) at nonzero frequencies, shown in the inset.

C. Superconductivity
We have probed for superconducting correlations in the vicinity of the QCP by measuring the superfluid density and the superconducting susceptibility, as a function of h and T . The superfluid density is given by [72,73] where  Here, the current density operator is given by j x (r i ) = describes the response of the system to a static orbital magnetic field B(q): where A x (q) = iB(q)/q y is the vector potential in an appropriate gauge.
Extrapolating K xx to the thermodynamic limit is challenging due to the presence of strong finite-size effects associated with the Fermi surface. These finite size effects are dramatically reduced by performing the computations in the presence of a weak a uniform magnetic field [74]. We adopt this procedure, but make the field spin-dependent, with the total flux through the system being 2π (−2π) for spin up (spin down). As shown in Ref. [49], such a procedure does not introduce a sign problem for the Monte Carlo. Moreover, in the limit of large system size, the magnetic field vanishes. (See Appendix A for the implementation of this effective magnetic field in the simulations.) Fig. 15 shows K xx as a function of q y , for different system sizes at temperature T = 0.05 and h = h c . The data collapse onto a single curve that extrapolates to a value near zero in the limit q → 0. The extrapolated value is clearly much below the universal value at the BKT transition, ρ BKT s = 2T /π. This implies that at h = h c , the system is not superconducting down to T = 0.05t. Similarly, the superfluid density for h = h c is found to be below ρ BKT for all T ≥ 0.05t. Note that, from Eq. (8), lim qy→0 K xx /q 2 y is the orbital susceptibility. From Fig. 15, we see that in our system, this quantity is negative, hence the orbital response is paramagnetic [75]. This indicates that there are no strong superconducting fluctuations (which would produce diamagnetism).
We have also computed the superconducting susceptibilities in the s-wave (A 1g ) and d-wave (B 1g ) channels, as a function of h and T : Here, ∆ s (r i ) = c i↑ c i↓ and ∆ d (r i ) = j η ij (c i↑ c j↓ − c i↓ c j↑ ), where η ij is defined below Eq. (4). The superconducting susceptibilities are displayed in Fig. 16. As expected for a non-superconducting state, both P s and P d saturate as a function of system size at fixed T and h [see Fig. 16(a,b)]. P s has a maximum at h ≈ h c , while P d is smaller and more or less independent of h.
The enhancement of P s by nematic fluctuations is further reflected in its behavior as a function of temperature at a fixed h [ Fig. 16(c,d)]. P s rises more rapidly than logarithmically at low temperatures. In contrast for non-interacting electrons P ∼ log(T 0 /T ), a functional form consistent with the T dependence of P d at low T . We interpret this as an indication that the s-wave susceptibility is enhanced by induced attractive interactions mediated by nematic fluctuations; the ground state for h near h c is likely an s-wave superconductor, although with T c < 0.05t.

D. Single Fermion Correlations
Having found that the system is not superconducting down to the temperatures reached in this study, we examine the single-fermion Green's function to glean information about the metallic state in the vicinity of the QCP. The metal is naturally characterized via the single particle density of states N (ω), the Fermi velocity v F , and the quasiparticle weight Z k F . However, measuring these quantities within DQMC is not straightforward, since they generally require real-frequency information inaccessible without analytic continuation. In addition, v F and Z k F are well defined only at zero temperature, while our simulations are at finite temperature. We partially resolve these difficulties by defining finite temperature objects N (T ), v F (T ), and Z k F (T ) which can be measured numerically, and whose zero temperature limits are N (ω = 0), v F , and Z k F respectively. Readers interested in the raw frequency dependence of the fermionic Green's function and self-energy are referred to appendix D.
Real frequency information can be extracted from our imaginary-time data using the identity [76]: where A(k, ω) is the fermion spectral function. G(k, β/2) is evidently equal to the spectral function integrated over a range of energies ∼ T . We therefore define We explain in Appendix C how these quantities attain the desired zero temperature limits. We also discuss in the appendix the procedures used to estimate these quantities in finite size systems with discrete momentum and imaginary time. Fig. 17 shows N (T ) for an L = 20 system as a function of h, for different temperatures. N (T ) remains nonzero for all values of h, consistent with a finite density of states in the limit T → 0. This indicates the absence of a substantial gap in the fermionic spectrum, even a gap with nodal points on the Fermi surface, down to T = t/12, corroborating the conclusion of Sec. III C that the system is not superconducting down to these temperatures. Moreover, the fact that N (T ) seems to neither vanish nor diverge at low temperature and h = h c has important implications for the nature of the metallic state of the QCP. Were the fermions to acquire a positive (negative) anomalous dimension, η F , the density of states would go to zero (diverge) at T = 0 [28,77]. We elaborate on this point in Sec. IV. Instead, the saturation of the density of states is consistent with a Fermi liquid, at least down to the temperature scales in our simulations. Note, as well, that the density of states appears to be depressed with increasing nematic order for h < h c The results for Z k F (θ, T ), where θ is the angle between k F and the x axis, are shown in Fig. 18   of h ≥ h c . At θ = 45 • , Z k F is close to unity at all temperatures. This can be understood as a consequence of the fact that at θ = 45 • , the coupling between fermions and the nematic critical modes vanishes by symmetry; hence, at these "cold spots", the effects of scattering off the critical modes are minimal. Z k F decreases systematically as θ approaches 0 • , and as h approaches h c . Z(k F , θ = 0, T ) decreases with decreasing temperature with negative curvature. On the basis of the present data, it is difficult to say whether it extrapolates to a finite value in the T → 0 limit (indicating a FL ground state), or vanishes with a small exponent, Z ∼ T a Z where a Z ≈ 0.1 (indicating a breakdown of FL behavior). It would be difficult to reconcile the data with a Z larger than about 0.15. Fig. 19(a) shows the same data as a function of h ≥ h c for θ = 0 • and θ = 45 • , emphasizing the decrease in Z k F as h approaches h c . Measuring Z k F close to the QCP for h < h c (in the nematic phase) has proven difficult, since the two configurations of the order parameter correspond to two different Fermi surfaces with a small splitting between them. In Fig. 19(b) we show the renomalization of the Fermi velocity, v F /v F,0 , where v F,0 is the Fermi velocity for the non-interacting band structure (α = 0). As with Z k F , v F is hardly renormalized at the cold spot θ = 45 • . In contrast, at θ = 0 • , v F is strongly reduced from its non-interacting value, and the velocity renormalization increases as h approaches h c . However, down to T = 0.1t, there is no obvious indication that v F vanishes (i.e., the effective mass diverges) as T → 0, even at h = h c .
Over the temperature range explored in this study, the single fermion Green's function is consistent with a renormalized Fermi liquid, with an angle-and h-dependent quasi-particle weight and velocity. However, the expected differences between a Fermi liquid and a "marginal Fermi liquid" or a weak (small a Z ) non Fermi liquid are subtle at non-zero T . We can rule out a non-Fermi liquid with a substantial value of a Z and/or η F , but not one with a z and η F sufficiently small.

IV. DISCUSSION
We view our results as a numerical "experiment." For an explicit lattice Hamiltonian, we have obtained detailed and extensive "measurements" of physical quantities. In the neighborhood of the quantum phase transition we find an enhanced tendency to superconductivity, substantial scattering of the fermionic quasiparticles, and a clearly delineated quantum critical fan. Moreover, for a large range of parameters, a simple scaling form provides a good description of the thermodynamic nematic fluctuations. In laboratory experiments on solid state systems, good scaling collapse over a similar dynamical range is often interpreted in terms of universal critical exponents. However, in our system, the fermion bilinear quadrupolar correlations are only consistent with the apparent scaling form of the nematic correlator over a narrow region as a function of (T, h − h c , q). This is reason to question whether the behavior we are seeing is characteristic of the asymptotic approach to a quantum critical point. It is possible that the success of the scaling analysis is fortuitous. Alternatively, the observed scaling may reflect an intermediate asymptotic regime, possibly associated with an unstable multicritical point. It is also possible that corrections to universal scaling behavior are simply smaller for D than for Q, in which case the scaling form of D might reflect the true critical behavior. At the very least, the observation of simple power laws invites a comparison with theories of metallic quantum criticality.

A. Comparison with theory
The theory of metallic QCPs, along with the closely related theory of metals coupled to fluctuating gauge fields [78][79][80], has traditionally been treated by integrating out the fermions to obtain an effective action for the bosonic modes that is highly non-local in both space and time. In this "Hertz-Millis" approach, a renormalization group (RG) analysis is then undertaken for the effective field theory so defined. While much can be learned from this program, there are at least two general causes for concern: 1) The analyticity of the β-function, which is a core ingredient of RG, is typically ensured by the locality of the action. For a non-local action, all aspects of the analysis, especially the enumeration of potentially relevant interactions that can be generated under RG, need to be tested explicitly. 2) In many cases, including in particular the 2D nematic QCP [71,81], if one considers the back effect of the collective fluctuations on the fermions, one finds that Fermi liquid theory is perturbatively unstable. The question then arises concerning the self-consistency of the effective action which is obtained by integrating out an assumed Fermi liquid.
Nonetheless, there is a general feeling that this approach is reliable in large enough spatial dimension [82] [83], and its main results are similar in structure to those of many subsequent studies. Accordingly, we briefly re-view the main features of the Hertz-Millis treatment of the Ising nematic QCP.
In this context, the quadratic part of the effective action for the nematic mode is of the form where r and κ are (renormalized) functions of T and h, and Q 0 is the suitably defined particle-hole response function (with the zero-frequency piece subtracted) where g(k, q) is a dimensionless form factor of appropriate symmetry, whose form depends on microscopic details ( e.g. g(k, q) = cos(k x ) − cos(k y )). The essential features of this expression relate to its asymptotic (IR) behavior in the limit E F |v F (k q ) × q| |ω n |, wherek q is the point on the Fermi surface (assumed to be unique mod inversion) at which the Fermi velocity, v F (k q ), is perpendicular to q. (The fact that the integral in Eq. (17) for given q is dominated by the neighborhood of this point on the Fermi surface is the essential observation motivating the "patch" theories of this problem.) The Hertz-Millis approach implies that the dynamical exponent z = 3 at tree level. From this, a naive scaling analysis (neglecting the possibility of higher-order nonlocal interactions) concludes that, since d + z = 5 is well above the Ising upper critical dimension, mean-field exponents ν = 1/2, γ = 1, and η = 0 apply, consistent with what we find. (See Table I.) Moreover, one would expect violations of scaling (for instancez = z and λ = 1/νz). Indeed in Refs. [3,16,68,69] it is shown that at criticality, such an analysis leads to r ∼ T log[E F /T ]. As we do not have the dynamical range to detect a logarithm, this, too, is consistent with our results.
However, on the face of it, our results for the dynamical nematic response are inconsistent with a dynamical exponent z = 3 and with the sort of anisotropies one might expect from the patch construction. Moreover, the relative weakness of any deviations of the single-particle properties from Fermi liquid behavior appears inconsistent with the Hertz-Millis predictions [84]. On the other hand, at the smallest momenta, deviations of D from the scaling form, A, Eq. (6), become apparent. This may indicate a crossover [85] from an intermediate z ≈ 2 regime to a different behavior in the deep infra-red [86].
The foregoing discussion of the relation between our data and Hertz-Millis theory can be applied to numerous subsequent attempts to bring the problem under theoretical control. Among the approaches employed by these studies, whose results have similar structure to those of Hertz-Millis, are large N F methods (where N F is the number of fermion flavors) [17,20,21], dimensional [34,87,88] and dynamical [18,19,30] regularization, and higher dimensional bosonization [89,90]. The related problem of orbital loop current criticality in a metal has been argued [11,91] to lead to a critical bosonic correlator with |ω| frequency dependence (as seen in our results), but with an asymptotic absence of momentum dependence (not seen in our results) characteristic of "local quantum criticality" [92]. In the vicinity of d = 3, several varieties of intermediate asymptotics corresponding to different schemes of large N F and N B (the number of boson flavors) have also been explored [33,93], at least one of which has z = 2 upon naive extrapolation to d = 2. We believe our results motivate revisiting existing theories of nematic quantum criticality in metals. In particular, the isotropy of bosonic correlations we observe appears to be at odds with "patch" constructions.
Finally, there has been considerable theoretical interest in the issue of superconductivity in the neighborhood of a nematic quantum critical point. It has been argued [37][38][39]43] that superconductivity with a critical temperature of order 1 times microscopic scales should be expected near the nematic quantum critical point. Indeed, it has been suggested that the superconducting T c at criticality may be so high that it preempts the quantum critical regime in which non-Fermi liquid behavior would be expected. More modestly, in the small α limit, it has been shown [42] that nematic fluctuations do indeed mediate attractive pairing interactions, and that the resulting transition temperature grows singularly as the quantum critical point is approached either from the ordered or the disordered side. In agreement with this latter result, the superconducting susceptibilities shown in Fig. 16 do show a maximum at h = h c . However, if T c at criticality is a number of order 1 times the Fermi energy, that number is apparently sufficiently small that T c < E F /70 for our chosen parameters. At larger values of the boson-fermion coupling constant, preliminary results [94] show that the ground state is an s−wave superconductor whose T c is maximal near the nematic QCP.

B. Relation to experiment
There is increasingly extensive evidence that electron nematic phases are common in highly correlated electronic fluids. [95] In particular, it has been suggested that a nematic QCP occurs near optimal doping in both the cuprates [96,97] and in the Fe-based superconductors [98,99]. However, the experimental situation in all these cases is complicated, and in some cases controversial.
It is premature to attempt any sort of serious comparison between the present results and the experiments in these materials. However, the experimental case is clearest in certain Fe based superconductors, so it is worth noting a few points of comparison. Most strikingly, elastoresistance [62,67], Raman scattering [63,66], and elastic constant [65] measurements show a large range of T and doping over which a large nematic susceptibility can be documented with a remarkably systematic Currie Weiss T dependence χ ∼ 1/[T − T * (x)] where T * (x) appears to depend roughly linearly on "doping" concentration, x, and to pass through zero at a critical doping concentration, x c , which approximately coincides with the optimal doping for superconductivity. It is impossible not to be encouraged by the similarity between this experimental finding and our numerical results.
Note added: After our paper was posted, a new theoretical analysis of the finite T dynamics near a nematic metallic QCP [100] has found a possible route to understanding the observed behavior of D, and in particular for the apparent dynamical exponent z = 2. wherê is a row vector of spin 1/2 fermionic creation operators for the spatial sites 1, 2, ..., L 2 , where L is the linear dimension. H h , H V , and H µ are the terms in Eq. 3 proportional to h, V, and µ respectively. We use a checkerboard decomposition to describe the kinetic energy matrices K (m) , whose elements are K i.e, m = 1, 2, 3, 4 enumerates the horizontal/vertical bonds originating from a site with an even/odd index. We allow t i,j to depend on spin in order to implement the spin-dependent magnetic field discussed below.
Plugging in unity operators in the τ ij sector at every time slice and taking the trace over the fermions, we bring the partition function to a form which can be sampled using Monte-Carlo techniques: n is a matrix which depends implicitly on the c-numbers τ i,j;n = ±1 through the relation K i,j (σ)(1 − ατ i,j;n )δ σ,σ , and the bosonic part of the action is given by i,j ,n τ i,j;n τ i,j;n+1 We have used ∆τ = 0.05 in our simulations.
Monte-Carlo sampling is most efficient when the determinant in (A3) is positive. Here, the absence of a sign problem is be guaranteed by microscopic time reversal symmetry. To see this, note that the matrices T n are block diagonal in the spin sector, and so the total determinant is the product of the spin up and spin down determinants. If time reversal is preserved , i.e then the product of the determinants is clearly positive. It has long been known [74] that finite size effects in fermionic systems, in particular at low temperatures, can be greatly improved by the application of a weak orbital magnetic field. On a finite size system with periodic boundary conditions, the field must be such that an integer multiple of the flux quantum passes through the system, Φ = nΦ 0 . In order to preserve time reversal symmetry as required in (A5), in this work we have applied the opposite field for both spin species, such that Φ ↑ = −Φ ↓ = Φ 0 . Such a field vanishes in the thermodynamic limit.
Although our simulations do not suffer from the sign problem, close to the quantum critical point our simulations do suffer from critical slowing down. To avert this, we found it was helpful to include global Monte-Carlo updates, in addition to local (Metropolis) ones. Our global updates consisted of a slight modification to the Wolff algorithm [101]: Starting from a space-time configuration of pseudo-spins, a cluster of pseudo-spins is constructed using the usual Wolff algorithm, supposing our action is given by (A4), i.e without taking into account the coupling to the fermions. Next, we propose flipping the pseudo-spins which belong to the cluster, and accept the move with a probability which is the ratio of the determinants in (A3) between the two configurations. It can easily be seen that such a move obeys detailed balance. The finite temperature phase boundary T N (h) (or equivalently h N (T )) can be computed by using classical finite size scaling techniques for two dimensional Ising transitions, characterized by susceptibility exponent γ = 7/4 and and correlation length exponent ν = 1. At a given temperature and for a given system size L, the most singular part of the thermodynamic susceptibility, χ s , satisfies h N can therefore be identified using a single parameter data collapse, as illustrated in Fig. 20. At low temperatures, the region of parameter space corresponding to classical criticality narrows, so that larger system sizes are needed to reliably estimate h N . Here, |n is a many-body eigenstate with energy ε n , and Ξ is the grand partition function. In this representation, the fermion spectral function is If the ground state is a Fermi liquid, we can extract the quasi-particle weight and Fermi velocity renormalization from Eq. (C3). We assume that at low temperatures, the spectral function in the vicinity of the Fermi surface obtains a Fermi liquid form, A(k, ω) = A qp (k, ω) + A reg (k, ω), where A qp is the quasi-particle contribution and A reg (k, ω) is a regular background. In the limit T → 0, A qp becomes a delta-function like peak: is the quasi-particle dispersion near the Fermi surface, Z k F is the quasi-particle weight, v F is the Fermi velocity. At finite but small T and ε k , the quasi-particle peak at k F obtains a width that scales as max(T 2 , ε 2 k ). The regular part satisfies A reg (k, ω) = O[max(T 2 , ω 2 , ε 2 k )].
Inserting the Fermi liquid form of A(k, ω) into Eq. (C3), we obtain that at low temperatures, We can use this relation to estimate Z k F and v k F along the Fermi surface. For a point on the Fermi surface, . The dispersion can be estimated from ε k = −d log[G(k, τ )]/dτ | τ =β/2 , and the Fermi velocity is then evaluated according to v k = ∇ k ε k . Some caution is needed because in our finite size systems, the values of k are quantized, and the k grid does not in general intersect the Fermi surface. One can instead estimate Z k from the points nearest to the Fermi surface. Also, our measurements of G(k, τ ) are taken at finite temperatures, yielding the finite temperature estimators N (T ), v F (T ), and Z k F (T ) discussed in the text (Equations 13, 14, and  15). An extrapolation to the limit T → 0 is necessary in order to deduce the "true" Fermi liquid properties, to the extent they are well defined.

Appendix D: Fermion Green's function and self-energy
In Fig. 21 we show the Matsubara frequency dependence of the inverse fermion Green's function. Linear frequency dependence, consistent with a Fermi liquid, seems to dominate down to the lowest frequencies available. Finer details are visible by considering the self energy (right panels) which shows both a noticeable temperature dependence and a deviation from linearity for θ = 0 • . Appendix E: Additional data Fig. 11 of the main text shows the momentum dependence of D −1 and Q −1 at h ≈ h c and temperatures T ≥ 0.1t. The momentum dependence at low temperatures is shown in Fig. 22, here plotted versus |q| instead of |q| 2 . Whereas D −1 is quadratic in |q|, Q −1 is apparently linear in |q| at low temperatures over a substantial range of momenta, a finding we have not yet understood.
Figs. 12 and 14 of the main text show D −1 and Q −1 versus frequency for a variety of momenta, at low temperature and h ≈ h c . Fig. 23 shows similar data for h ≈ 1.3h c , illustrating that deviation from criticality appears mostly through an additive shift of the inverse correlator. Fig. 8 of the main text shows the validity of the functional approximant A of Eq. (6) in describing thermodynamic nematic correlations D(h, T, q, ω n = 0). Fig. 24 plots D −1 vs A for all frequencies. The agreement is very good except for the smallest momenta, shown in blue and black.  6). On the plot of Q −1 at right, the dotted lines are a linear fit to the data at nonzero momentum, with a temperature-independent slope. The apparent linear momentum dependence of Q −1 suggests a critical power law, though with an apparently different exponent than the one found in D −1 . Fig. 13 of the main text shows the position dependence of the equal time nematic correlator D. The analogous data for the equal time quadrupolar correlator Q are somewhat more messy. This can be seen in Fig. 25 where we plot the dependence ofQonr for r along the symmetry directions, (1, 0) and (11); as is apparent, the error bars (associated with the sensitivity of the results to boundary conditions) are sufficiently large that it is difficult to conclude much from this figure. Here, there is a non-vanishing energy scale associated with nematic fluctuations since we are deep in the disordered phase, but the frequency dependence is similar to that at h = hc (see Fig. 12) . Different colors represent different values of q, and filled symbols mark frequencies |ωn| < vF |q|. The insets show the subset of the data with T = 0.025 and momenta along the (10) direction shifted by their zero frequency values. For the inset of the bottom panel the data at q = 0 are scaled down by a factor of 0.01.
[21] A. Abanov    FIG. 25. The inverse of the equal time quadrupolar correlator at h = hc and T = 0.1, plotted versus the square of the spatial separation r along high symmetry directions for several system sizes. In contrast to the data of Fig. 13, substantial tetragonal anisotropy can be seen by comparing the left panel (horizontal displacements) and the right panel (diagonal displacements). However, the data apparently represent the thermodynamic limit only for |r| 4, and therefore we cannot reliably extract long-distance behavior.