Unconventional non-Fermi liquid state caused by nematic criticality in cuprates

At the nematic quantum critical point that exists in the $d_{x^2-y^2}$-wave superconducting dome of cuprates, the massless nodal fermions interact strongly with the quantum critical fluctuation of nematic order. We study this problem by means of renormalization group approach and show that, the fermion damping rate $\left|\mathrm{Im}\Sigma^R(\omega)\right|$ vanishes more rapidly than the energy $\omega$ and the quasiparticle residue $Z_f\rightarrow 0$ in the limit $\omega \rightarrow 0$. The nodal fermions thus constitute an unconventional non-Fermi liquid that represents an even weaker violation of Fermi liquid theory than a marginal Fermi liquid. We also investigate the interplay of quantum nematic critical fluctuation and gauge-potential-like disorder, and find that the effective disorder strength flows to the strong coupling regime at low energies. Therefore, even an arbitrarily weak disorder can drive the system to become a disorder controlled diffusive state. Based on these theoretical results, we are able to understand a number of interesting experimental facts observed in curpate superconductors.


Introduction
A large amount of experimental and theoretical studies have been devoted to studying the unusual properties of high temperature cuprate superconductors in the past thirty years [1,2,3,4,5,6,7,8,9]. Although some consensuses have been reached, many fundamental problems are still in debate, including the microscopic pairing mechanism [1,2,7,8,9], the origin of pseudogap [2,6], and the description of non-Fermi liquid behaviors of the normal state [1,2,3]. In the past decade, there have been accumulating experimental evidences for the existence of a strong anisotropy in many of the physical properties of YBa 2 Cu 3 O 6+δ (YBCO) [10,11,12] and Bi 2 Sr 2 CaCu 2 O 8+δ (BSCCO) [13,14]. Such an anisotropy is widely believed to be driven by the formation of a novel electronic nematic order [15,16,17,18,19], which spontaneously breaks the C 4 symmetry down to a C 2 symmetry. In case the nematic transition line goes across the superconducting transition line and penetrates into the superconducting dome, there exists a zero-temperature nematic quantum critical point (QCP). The nematic quantum phase transition and the associated quantum critical behaviors have been investigated extensively in recent years [20,21,22,23,24,25,26,27,28,29,30,31].
From a theoretical perspective, there are two widely studied scenarios to induce an electronic nematic order. First, the nematic order can be generated by melting a stripe order that spontaneously breaks both translational and rotational symmetry [15,16,17,18,19]. The other way is related to Pomeranchuk instability which refers to the deformation of the shape of the Fermi surface of a metal due to Coulomb interaction [15,16,17,18,19,32,33,34,35]. In the simplest case, Pomeranchuk instability occurs when the circular Fermi surface of a two-dimensional metal becomes ellipselike via quadrupolar distortion. [15,16,17,18,19]. The Hubbard model defined on square lattices [35,39,40,41,42,43,44] provides a pertinent platform to investigate the electronic nematic order. Halboth and Metzner [35] studied a two-dimensional Hubbard model using functional renormalization group (RG) method, and revealed Pomeranchuk instability and nematic order. More recently, the square lattice Hubbard model is studied by variational cluster approximation [45,46] and found to display a local nematic phase under certain circumstances, which might be applied to understand the intra-unit-cell electronic nematicity observed in the scanning tunneling spectroscopy (STS) measurements by Lawler et al [13].
Another context of studying nematic order is provided by the superconducting dome of cuprate superconductors. For a pure d x 2 −y 2 -wave superconductor, the discrete C 4 symmetry is preserved. However, when superconductivity coexists with a nematic order, the gap nodes are shifted from their original positions and the C 4 symmetry is broken down to C 2 [20,22,23,24]. Such an anisotropic superconducting state is physically equivalent to a d x 2 −y 2 + s-wave superdconducting state [20,22]. At the nematic QCP, the massless fermions excited from the d x 2 −y 2 -wave gap nodes couple strongly to the quantum critical fluctuation of nematic order parameter, which can be effectively described by a (2+1)-dimensional field theory [20,21,22,23,24,25,26,27,28,29,30,31]. This model was first analyzed by Vojta et al. [20,21,22], who made an ǫ-expansion and argued that the nematic phase transition is turned to first order. Kim et al. [23] later tackled the same model by means of 1/N-expansion, where N is a large fermion flavor, and concluded that the transition remains continuous. Huh and Sachdev [24] performed a renormalization group (RG) analysis, and showed that the fermion velocity ratio v ∆ /v F → 0 in the lowest energy limit, where v F is the Fermi velocity of nodal fermions and v ∆ the gap velocity [20,21,22,23,24,25,26,27,28,29,30,31]. The unusual velocity renormalization leads to significant changes of a number of quantities, including density of states (DOS) [25], specific heat [25], low-T thermal conductivity [26], superfluid density [28], and London penetration depth [31].
In this article, we revisit the issue of unusual physical properties caused by the strong interaction between massless nodal fermions and critical nematic fluctuation. We will apply the powerful RG approach to calculate the fermion damping rate ImΣ R (ω) , where ImΣ R (ω) is the imaginary part of retarded fermion self-energy, and the corresponding quasiparticle residue Z f . Kim et al. [23] have previously computed the fermion self-energy and spectral function. Interestingly, it will become clear below that the RG analysis give rise to different results once the singular renormalization of fermion velocities is taken into account. In particular, we will show that the fermion damping rate vanishes upon approaching the Fermi level more rapidly than the energy ω, namely lim ω→0 ImΣ R (ω) /ω → 0. According to the conventional notion of quantum many-particle physics, one would expect the system to behave as a normal Fermi liquid. However, by analyzing the RG results, we find that the quasiparticle residue vanishes, i.e., Z f → 0, in the limit ω → 0. Therefore, the system is actually a non-Fermi liquid that represents an even weaker violation of Fermi liquid theory comparing to a marginal Fermi liquid (MFL). To the best of our knowledge, this type of unconventional non-Fermi liquid behavior has not been reported previously.
In realistic materials, there are always certain amount of disorders, which may play an important role. As demonstrated earlier by Nersesyan et al. [47,48], the most important disorder in cuprates behaves like a randomly distributed gauge potential. Thus, we will mainly study the influence of random gauge potential on the physical properties of nodal fermions and also the stability of nematic QCP. In the absence of nematic order, the coupling between nodal fermions and random gauge potential has attracted considerable interest [47,48,49]. Here, we consider the case in which nodal fermions couple to both the nematic order and random gauge potential, and then study the interplay of these two interactions by means of RG method. We find that the effective strength of gauge potential disorders tends to diverge at the lowest energy. This behavior signals the emergence of a finite zero-energy DOS ρ(0) and the happening of quantum phase transition from an unconventional non-Fermi liquid to a disorder dominated diffusive state. The nodal fermions acquire a finite scattering rate γ, which in turn affects the thermodynamic and spectral behaviors of nodal fermions. The RG results for the self-energy of nodal fermions can be used to understand a number of experimental facts observed in cuprate superconductors. We will show that the RG results are qualitatively consistent with some recent measurements of specific heat, fermion damping rate, and temperature dependence of fermion velocities.
The rest of the paper will be organized as follows. We first present the effective low-energy field theory for the interaction between nematic order and nodal fermions in section 2. The random gauge potential is also introduced in this section. In section 3, we make detailed theoretical analysis for the self-consistent RG equations of the fermion velocities in the clean case. Based on the RG solutions, we proceed to compute the fermion damping rate, quasiparticle residue, and other physical quantities. We will show that the quantum critical fluctuation of nematic order leads to unconventional non-Fermi liquid behaviors of nodal fermions. We consider the influence of random gauge potential in section 4 and find that the effective disorder strength flows to the strong coupling regime at low energies. Therefore, even weak disorders play a significant role and drive the system to enter into a disorder controlled diffusive state. In section 5, we discuss the possible application of the RG results to some experimental findings of cuprates. In section 6, we briefly summarize our main results.

Effective field theory
We start from an effective action S = S Ψ + S φ + S Ψφ . The free action for the nodal fermions is [20,21,22,23,24] where τ (x,y,z) denote the Pauli matrices. The two component Nambu spinors are defined f 3a and f 4a represent fermions excited from the nodal points (K, K), (−K, K), (−K, −K), and (K, −K) respectively [20,21,22,23,24]. The repeated spin index a is summed from 1 to N, where N is the number of fermion spin components with the physical value being 2. The action S φ describes the nematic order parameter, which is expanded in real space as where τ is imaginary time and c velocity of field φ. It is convenient to choose c = 1. Mass parameter r tunes a quantum phase transition from d x 2 −y 2 -wave superconducting state to a state where the superconducting and nematic orders coexist, with r = 0 defining the QCP. Moreover, u 0 is the quartic self-interaction strength. The nematic order couples to fermions via a Yukawa-coupling term where λ 0 is the coupling constant. In addition to the nematic order, the nodal fermions also couple to gauge-potential type disorders, which is described by with Γ = (τ z , τ x ) and v Γ = (v Γ1 , v Γ2 ). Here, the random gauge potential A(x) is assumed to be a Gaussian white noise, defined by where g is the impurity concentration and v Γ measures the strength of a single impurity. We will follow Huh and Sachdev [24] and perform RG analysis by employing 1/Nexpansion. The inverse of free propagator of φ behaves as q 2 + r. After including the polarization, there will be an additional linear-in-q term. The q-term dominates at small q over q 2 -term, so the q 2 -term can be neglected. Near the nematic QCP, we keep only the mass term and rescale φ −→ φ/λ 0 and r −→ N f rλ 2 0 , leading to The free propagators for fermions Ψ 1a and Ψ 2a are written as respectively. At the nematic QCP, r = 0, the propagator for the nematic order field φ is where the polarization Π(Ω, q), to the leading order of 1/N-expansion, is given by After carrying out RG calculations [24,27], we find that all the physical parameters flow with a varying length scale l as follows where we have introduced a quantity The expressions of C 1,2,3 can be found in the Appendix A. We emphasize that the impact of disorders is characterized by the quantity C g , rather than v Γi with i = 1, 2. It will be shown below that C g flows to strong coupling at large l, leading to remarkable physical properties.

Unconventional non-Fermi liquid behaviors of nodal fermions
It is well established that conventional metals can be described by the standard Fermi liquid theory, which states that the fermionic excitations of a normal Fermi liquid must have a sufficiently long lifetime and exhibit a sharp quasiparticle peak in their spectral peak despite the existence of Coulomb interaction [50]. The conventional notion is that the fermionic quasiparticles constitute a normal Fermi liquid if their zero-T damping rate ImΣ R (ω) vanishes more rapidly than ω in the limit ω → 0. This criterion can be mathematically expressed as Another criterion to identify normal Fermi liquid is to define an important quantity: the quasiparticle residue, also called renormalization factor, Z f . The residue Z f is usually calculated through the definition where the real part of retarded fermion self-energy ReΣ R (ω) is related to ImΣ R (ω) via the Kramers-Kronig (K.-K.) relationship. The residue Z f is finite in a normal Fermi liquid, but vanishes in a non-Fermi liquid. Generically, the fermion damping rate in an interacting fermion system can be formally written as where C F is a constant. With the help of K.-K. relation, the corresponding real part of the fermion self-energy, in the low energy limit, is given by where ω 0 is a cutoff, and I(x) is a function that depends only on x.
It can be easily checked that Z f = 0 for 0 < x ≤ 1 and Z f = 0 for x > 1. Therefore, the above two criteria are actually equivalent because Z f automatically takes a finite value whenever equation (17) is fulfilled.
However, in the present nodal fermion system, the above two criteria are no longer equivalent. We will show by explicit calculations in this section that the damping rate    Figure 1. of nodal fermions vanishes more rapidly than the energy, but the residue Z f vanishes, i.e., Z f → 0.

Fermion damping rate and quasiparticle residue
Now we calculate the fermion damping rate and the residue Z f utilizing the solutions of the RG equations (11)- (15). The unusual renormalization of fermion velocities need to be taken into account in an appropriate manner. We only present the results obtained in the clean limit C g = 0 in this section, and include the effect of random gauge potential in the next section. After solving RG equations (11)-(13) self-consistently, we show the l-dependence of fermion velocities v F , v ∆ and ratio v ∆ /v F obtained at different initial values of ratio v ∆0 /v F 0 in figures 1(a), (b) and (c) respectively. All the quantities v F , v ∆ , and v ∆ /v F decrease with growing l and flow eventually to zero as l → +∞, but apparently v F decreases much more slowly than v ∆ . If we use the ratio v ∆ /v F to characterize the velocity anisotropy, it is clear that the nematic order drives an extreme velocity anisotropy v ∆ /v F → 0. For later use, it is helpful to extract an approximate analytical expressions for the velocities. Considering the leading and sub-leading terms, the with c s ≈ 0.3809/N, which is consistent with the expression in reference [24]. In the long wavelength limit, l → ∞, ln(c s )/ ln(l) → 0, so the sub-leading term can be ignored. Retaining only the leading term, the asymptotic behavior of velocity ratio can be well approximated by Substituting the asymptotic form of v ∆ /v F into equation (11), we obtain with c 1 ≈ 0.078 for large l. Solving this equation, we express the renormalized with c 2 = π 2 c 1 8 . It can also be found that v ∆ at large l behaves approximately as According to equations (24) and (25), both v F and v ∆ vanish as l → +∞.
To examine the impact of singular fermion velocity renormalization and extreme velocity anisotropy on the properties of nodal fermions, we will compute a number of physical quantities. Coming first is the quasiparticle residue Z f . Apart from the widely used definition given by equation (18), the residue Z f can also be calculated within the RG framework. The interaction induced renormalization of fermion field Ψ is encoded in the residue Z f , which exhibits the following l-dependence or alternatively The l-dependence of Z f can be easily obtained from the above equation, and is presented in figure 2(a). We find that Z f flows to zero in the limit l → +∞, which indicates the breakdown of normal Fermi liquid and the absence of well-defined Landau quasiparticles. However, Z f decreases very slowly with growing l, thus the deviation of the system from a normal Fermi liquid ought to be quite weak. To see this point, we plot the l-dependence of Z f l in figure 2(b), and find that Therefore, the residue Z f obtained at nematic QCP vanishes with growing l more slowly than that of a MFL [51], where Z f ∼ 1/l. According to equation (18), one can speculate since a large length scale l corresponds to a low energy ω. This speculation in turn implies that the imaginary part of the retarded fermion self-energy ImΣ R nem (ω) should display the following low-energy behavior, To make the above discussions more quantitative, we are now going to make a detailed analysis of the asymptotic behavior of Z f . At large running scale l, the equation of Z f can be approximated by where c 4 ≈ 0.426 and c 5 = π 2 c 4 8 ≈ 0.523. Solving this equation leads to the asymptotic behavior in the limit l → +∞. To proceed, it proves convenient to utilize the relationship between ω and l: where ω 0 is a cutoff. Then the real part of retarded self-energy is approximated by as ω → 0. Using the K.-K. relation, we obtain the imaginary part of self-energy From equation (35), it is easy to verify that this fermion damping rate is smaller than that in a MFL and manifests the asymptotic behavior lim ω→0 ImΣ R nem (ω)/ω → 0, confirming the above analysis based on numerical results. According to the conventional notion of quantum many-body physics, one would expect the system to behave like a normal Fermi liquid. However, in the low energy limit ω → 0, the residue Z f actually vanishes: which clearly implies that the system under consideration is actually a non-Fermi liquid and the nodal fermions do not have a well-defined quasiparticle peak in the spectral function. We see that the above residue Z f flows to zero at a lower speed than that of a MFL, i.e., Z f ∼ 1 ln(ω 0 /ω) . Therefore, the quantum critical fluctuation of nematic order gives rise to an even weaker violation of ordinary Fermi liquid theory than a MFL. To the best of our knowledge, this sort of unconventional non-Fermi liquid state has not been reported previously.
It is now interesting to compare this unconventional non-Fermi liquid state with graphene. In graphene, early perturbative calculations revealed that Dirac fermions exhibit MFL behavior with damping rate ImΣ R (ω) ∼ ω and vanishing Z f [53]. Nevertheless, subsequent careful RG studies [53,54,55,56] have found that, the fermion damping rate actually depends on energy as ImΣ R (ω) ∝ ω/ ln 2 ω at zero temperature due to the long-range Coulomb interaction, whereas the corresponding Z f flows to a finite value as ω → 0. Therefore, graphene is a normal Fermi liquid. In contrast, the massless nodal fermions constitute a non-Fermi liquid at the nematic QCP, because Z f vanishes in the lowest energy limit. The crucial difference between these two cases can be understood as follows. At nematic QCP, the fermion velocities are driven to vanish by the nematic order, so the effective fermion-nematic interaction is significantly pronounced. In graphene, however, the fermion velocity is dramatically enhanced at low energies by Coulomb interaction, which then weakens the effective Coulomb interaction and guarantees the validity of the Fermi liquid description. These two interacting Dirac fermion systems provide interesting new insight on the effects of strong electron correlations and also on the criterion of non-Fermi liquid states.
Moreover, an important lesson one can learn from the research experience of graphene is that, RG may lead to qualitatively different spectral properties of fermions from that obtained by ordinary perturbative expansion approach. Indeed, this is the main motivation that has promoted us to make an extensive RG analysis of the spectral properties of nodal fermions at the nematic QCP.

Density of states and specific heat
We have showed in the last subsection that the nodal fermions exhibit unconventional non-Fermi liquid behaviors at the nematic QCP. In this subsection, we will compute two important quantities, namely DOS and specific heat, on the basis of the RG solutions with the goal to gain a better understanding of the non-Fermi liquid state. The fermion DOS can calculated from the retarded Green functions of nodal fermions via the definition [25] In the absence of interactions, the fermion DOS is well-known to be linear in ω, namely ρ(ω) ∝ ω. This linear behavior will be changed once the interaction effects are considered. After including the fermion velocity renormalization, employing the method presented in [25], the flow equation for ρ(ω) is given by In cuprates, the initial value of velocity ratio v ∆0 /v F 0 is known to be much smaller than 1 [57]. Additionally, due to the quantum fluctuation of nematic order, v ∆ /v F decreases monotonously as the energy scale is lowering. Therefore, for a given initial value v ∆0 /v F 0 < 1, we only need to consider the case v ∆ < v F . In order to show that the conclusion is independent of the condition v ∆0 /v F 0 < 1, we also plot the curves for the case v ∆0 /v F 0 = 2 in figures 3 and 4. The RG equations of DOS and specific heat are also given with a generalized form. In the clean limit, we plot the results for ρ(ω) in figure 3. We see from figure 3(a) that ρ(ω) is apparently not linear in ω, but displays the asymptotic behavior On the other hand, figure 3(b) implies that This DOS ρ(ω) is qualitatively different from the power law function ρ(ω) ∼ ω 1−α , where α is a small finite value, obtained previously in reference [25]. In the non-interacting limit, the fermion specific heat depends on T as C V (T ) ∝ T 2 . As shown in Appendix B, including the influence of the quantum fluctuation of nematic order and random gauge potential leads us to At low T , C V (T ) behaves as which is visualized in figure 4(a). Apparently, the original quadratic T -dependence of C V (T ) obtained in the non-interacting limit is significantly altered. According to figure  4(b), we can express C V (T ) at very small T as which is also distinct from the power law T -dependence C V (T ) ∼ T 2−β , where β is a small finite constant, obtained previously [25]. Simple analysis reveal that the three parameters C 1 , C 2 , and C 3 appearing in the RG equations (38) and (41) all flow to zero in the lowest energy limit [24]. This asymptotic behavior makes it impossible to express ρ(ω) and C V (T ) by power law functions. As given in Appendix C, the low energy behavior of ρ(ω) and the low-T behavior of C V (T ) can be approximately expressed as where a ρ and a C are two negative constants. At N = 2, a ρ ≈ −1.896 and a C ≈ −1.466. We notice that our RG equation for C V (T ) (41) is not exactly the same as that obtained by Xu et al. [25]. This difference does not affect our conclusion that C V (T ) is not a power law function at low T . Indeed, if we start from the RG equation of C V (T ) presented in [25], we reach the same conclusion. A more detailed discussion is presented in Appendix B and Appendix C.

Impact of Random gauge potential
In this section, we investigate the impact of random gauge potential on the quantum critical behaviors near the nematic QCP. We will first show the RG solutions obtained in the presence of random gauge potential, and then re-calculate the fermion DOS and specific heat after considering the influence of disorder scattering. To this end, we need to retain a nonzero C g in the RG equations (11)-(15), (27), (38), and (41), and then solve these RG equations numerically.

Flow of effective disorder strength
In order to make a direct comparison to the clean case in which C g = 0 and to explicitly see the influence of a nonzero C g on the running behaviors of various parameters, we plot the l-dependence of v F , v ∆ , v ∆ /v F , and Z f in figures 5(a)-(d). Moreover, we show the ω-dependence of ρ(ω) in figure 5(e) and the T -dependence of C V in figure 5(f), respectively. We will discuss under what circumstances these RG results are modified in the next subsection.
Comparing figures 3(a) and (b) with figures 5(a) and (b), we see that the detailed l-dependence of v F and v ∆ are both altered dramatically by the random gauge potential. As shown in figure 3(c) and figure 5(c), the velocity ratio v ∆ /v F exhibits exactly the same l-dependence in the clean and disordered cases, which originates from the fact that C g does not enter into the RG equation of velocity ratio v ∆ /v F . According to figure 2(a) and figure 5(d), the renormalization factor Z f flows to zero in the disordered case much more rapidly than the clean case.
As shown in figure 5(e), the DOS ρ(ω) is divergent in the lowest energy limit due to random gauge potential. This is completely different from the behaviors of clean case presented in figure 3. In figure 5(f), we plot the ratio between C Dis V /C Cl V , where C Dis V and C Cl V are the specific heat obtained in disordered and clean cases, respectively. Since the solutions of the RG equations about DOS and specific heat are modified substantially in the disordered case, it turns out that random gauge potential is a relevant perturbation in the present system. To verify the relevance of random gauge potential, we should appeal to the RG analysis of the effective disorder strength.     In the action term S dis given in equation (4), the parameters that characterize the fermion-disorder coupling seem to be v Γ1 and v Γ2 . It can be seen from the RG equations that v Γ1 flows in precisely the same way as v F , and that v Γ2 as v ∆ . Thus, both v Γ1 and v Γ2 are strongly renormalized and driven to vanish as l → +∞. However, this does not mean that the disorders can be simply neglected. Indeed, whether disorders are important is determined by the ratio between the interaction energy given by S dis and the fermion energy We can see that this ratio is defined by C g , which enters into the RG equations for the parameters v F , v ∆ , v Γ1 , and v Γ2 . The effective strength of random gauge potential should be measured by C g , rather than v Γ1 and v Γ2 . Recall that C g is a function of five parameters, i.e., v F , v ∆ , v Γ1 , v Γ2 , and g. Among these parameters, g is assumed to be a dimensionless constant, but the other four parameters flow strongly with the varying l. Detailed RG analysis revealed that C g goes to infinity, namely in the limit l → +∞ as a consequence of the singular renormalization of the fermion velocities ratio v ∆ /v F . More quantitatively, the large scale behavior of C g can be described by where c 3 = 4g It is necessary to explain here why the effective strength of random gauge potential should be represented by C g , rather than solely by the coefficients v Γi with i = 1, 2, and why C g → ∞ in the lowest energy limit. In an interacting fermion system, the effective strength of the interaction is characterized by the ratio between the interaction energy scale and the kinetic energy of fermions. This ratio is widely used in condensed matter physics to judge whether an interacting fermion system can be defined as a strongly correlated system or not. For instance, the normal metal with a high density of itinerant electrons is believed to be a weakly interacting system since the energy scale of Coulomb potential is much smaller than the Fermi energy. In a massless Dirac fermion system (such as graphene), the effective strength of long-range Coulomb interaction is defined as α ∼ e 2 v F , where e 2 appears in the action of the Coulomb interaction as a coupling coefficient and the fermion velocity v F reflects the order of the kinetic energy [54,56]. Another example comes from the effective BCS model of Dirac fermion systems, where the effective strength of pairing interaction is characterized by g ∼ u v F with u being the coupling coefficient of pairing interaction and v F fermion velocity [58,59].
This criterion also applies to disordered systems. When massless Dirac fermions couple to random gauge potential, the effective strength of random gauge potential should be defined by C g ∼ v 2 Γi v F v ∆ , rather than the coefficient of fermion-disorder coupling v Γi [47,48]. When v Γi flows to zero in the lowest energy limit, C g does not necessarily vanish since there is a possibility that v F and v ∆ may vanish more (or equally) rapidly than v Γi . When the nodal fermions couple to both the quantum fluctuation of nematic order and random gauge potential, the four parameters v F , v ∆ , v Γ1 , and v Γ2 all flow to zero in the lowest energy limit, but the effective strength of random gauge potential becomes very large, namely C g → ∞. This originates from the fact that are constants, while at the same time the velocity ratio v ∆ /v F → 0. The behavior C g → ∞ at low energies indicates that random gauge potential is relevant. To see this point, we neglect the nematic order and consider only the coupling of nodal fermions to random gauge potential, which results in simplified RG equations: Using the two relations In which all vanish rapidly in the limit l → ∞. Since C g is a constant, random gauge potential is marginal. From the above calculations, we can conclude that the behavior C g → ∞ obtained in the presence of both nematic fluctuation and random gauge potential is directly related to the extreme velocity anisotropy v ∆ /v F → 0 induced by the quantum fluctuation of nematic order. The flow of C g towards strong coupling regime is a clear signature that random gauge potential should have substantial physical effects on the low-energy behaviors of nodal fermions, which will be discussed in the next subsection.
We notice other interesting correlated electron models in which the coupling coefficient of an interaction vanishes at low energies, but the interaction is not negligible due to the more (or equally) rapid decrease of the kinetic energy of electrons. Recently, Sur and Lee studied the influence of quantum fluctuations of an antiferromagnetic (AF) order at an AF quantum critical point in a metal supporting one-dimensional Fermi surface [60]. In particular, they showed that the coupling coefficient of the interaction flows to zero at low energies. However, the fermion velocity also vanishes, thus the interaction cannot be simply neglected. Actually, Sur and Lee found that the interaction drives the system to become a so-called quasi-local strange metal that is apparently qualitatively different from the free fermion system.

Physical effects of random gauge potential
How should we understand the divergence of C g ? In order to answer this question, we first consider the non-interacting system that contains only nodal fermions and random gauge potential. In this case, the RG equation of fermion DOS becomes d ln ρ d ln ω with If C g0 < 1, we have where α satisfies 0 < α < 1. In this case, the RG equation for specific heat is given by The solution of this equation is given by If C g > 1, it is easy to get ρ(ω) → ∞ (56) in the limit ω → 0. The divergence ρ(ω) indicates the emergence of a disorder dominated diffusive state, in which a finite disorder scattering rate γ and a finite zero-energy DOS ρ(0) are generated [49,61,62]. The value C g = 1 defines a QCP for the quantum phase transition between the ballistic and diffusive states of nodal fermions. Therefore, a weak random gauge potential gives rise to power law behavior of ρ(ω), whereas a sufficiently strong disorder can trigger the quantum phase transition to a diffusive state.
In the presence of both nematic critical fluctuation and random gauge potential, the fact that C g → ∞ in the lowest energy limit signals the development of a disorder dominated diffusive state and the generation of a finite γ and a finite ρ(0) even in the case of weak random gauge potential. Although the perturbative RG method provides a powerful tool to judge whether and when a phase transition takes place, it cannot be used to compute the precise value of γ. To calculate γ, one usually needs to construct a self-consistent equation for the retarded fermion self-energy by properly considering both the fermion-disorder scattering and the influence of quantum critical fluctuation of nematic order [66]. This is an interesting yet complicated issue, which is beyond the scope of the present work and subjected to future research. Here, we use the RG method to make a rough estimation for the energy scale of γ. As shown in Fig. 5(e), for small given values of C g0 and v ∆0 /v F 0 , the solution of the RG equation of DOS should have the following properties: ρ(ω) decreases as the varying energy scale decreases, but tends to increase when the energy scale exceeds a critical value E c (C g0 , v ∆0 /v F 0 ), which is a function of C g0 and v ∆0 /v F 0 . The magnitude of γ should be an increasing function of E c . In the following calculations of ρ(0) and C V (T ), we will regard γ as an undetermined constant. Fortunately, the qualitative behaviors of ρ(0) and C V (T ) in the low energy regime is independent of the precise value of γ.
The imaginary part of retarded fermion self-energy can be generically written as where ImΣ R nem (ω) is the contribution induced solely by the nematic order. The disorderinduced scattering rate γ represents a characteristic energy scale. At energies higher than γ, namely ω > γ, ImΣ R nem (ω) dominates over γ and all the RG results for v F , v ∆ , Z f , ρ(ω), and C V (T ) shown in figures 5(a)-(f) are still applicable. At ω < γ, the diffusive motion of nodal fermions and its interplay with critical nematic fluctuation govern the low-energy properties of the system.
Once a finite γ is generated, the renormalized velocities v F and v ∆ no more vanish at low energies, which is apparently different from the clean case. Instead, as the energy scale decreases, both v F and v ∆ are saturated to certain finite values, denoted by v ′ F 0 and v ′ ∆0 , below the energy scale set by γ. Hence, there is no extreme velocity anisotropy in the diffusive state. The fermion DOS and specific heat also exhibit different behaviors comparing to the clean case. To demonstrate the difference in DOS, we write the retarded propagators of nodal fermions in the forms: Calculations find that the fermion DOS depends on ω as In In order to compute the specific heat, it is convenient to invoke the standard Matsubara formalism of fermion propagators, i.e., where ω n = (2n + 1)πT with n being integer. The fermion free energy is given by Summing over ω n leads to Now the fermion specific heat in the low energy regime can be approximately given by which depends on T linearly. We can see that DOS and specific heat obtained in the diffusive state exhibit entirely different behaviors from the unconventional non-Fermi liquid state below the energy scale γ.
We finally make a brief remark on the behavior of the system staying away from the nematic QCP. Suppose the system stays at a distance r to the QCP, the RG results are still valid and the fermion velocity ratio is still renormalized at energies higher than the scale corresponding to r. However, the renormalization terminates at certain low energy scale. Therefore, the effective strength of random gauge potential is moderately enhanced compared its bare value, though C g does not diverge.

Comparison with experiments
In this section, we address the possible connection between the theoretical results obtained in the last sections and the phenomenology of cuprate superconductors. We are particularly interested in three existing experimental findings about some of the unusual properties of the superconducting dome.

Anomalous residual linear-T term of specific heat in cuprates
We first discuss the residual specific heat in cuprates. Due to the line nodes of dwave gap, the specific heat in the superconducting phase of cuprates is expected to exhibit a T 2 behavior, i.e., C V (T ) ∝ T 2 , at low T . This expectation is generically consistent with experiments [67]. In the lowest T limit, experiments have observed a residual linear T term of C V (T ) [67,68,69,70], which is usually attributed to the finite fermion DOS generated by disorder scattering and is also well consistent with the result given by equation (63). There is, however, an unexpected experimental finding [67,70] that the residual specific heat of YBCO is obviously larger in magnitude than that of La 2−x Sr x CuO 4 (LSCO), although the former material is known to be cleaner than the latter. Apparently, disorder scattering alone is not capable of accounting for this experimental fact. Recently, coexistence of d-wave superconductivity with a loop current order was proposed to give a possible explanation [71,72,73] for the large residual linear-T term of specific heat in YBCO. In this scenario, when d-wave superconductivity coexists with a loop current order, two of four nodal points are converted to finite Fermi pockets of Bougoliubov quasiparticles, which then generates a finite ρ(0) and a residual linear-T term of specific heat [71,72,73,74]. The recent ultrasound measurements showed the possible evidence for the existence of loop current order in YBCO [75,76].
Here we propose an alternative explanation for the above anomalous experimental results of residual specific heat. Our RG analysis found that the effective strength of random gauge potential, being the most relevant disorder to cuprates [48], is strongly enhanced by the critical nematic fluctuation, which in turn increases the residual value of the specific heat. To understand the role of quantum nematic fluctuation, we firstly consider only the coupling of nodal fermions to random gauge potential. In this special case, C g does not flow and thus remains a constant, namely C g = C g0 . If C g0 is very small, which means the system is only slightly disordered, the behavior of ρ(ω) and specific heat would be governed by equations (53) and (55), respectively. In this case, the system does not develop a finite ρ(0) and there is no residual liner-T term of specific heat. If the system is quite disordered such that C g0 > 1, it enters a diffusive state with a finite scattering rate γ, which induces a finite ρ(0) and a residual linear-T term of specific heat. Apparently, the residual specific heat is larger in more disordered systems, which is not consistent with the aforementioned experiments of residual specific heat. If we consider both the quantum fluctuation of nematic order and random gauge potential, C g is significantly enhanced and flows to strong coupling at low energies even if C g0 takes an arbitrarily small value. This implies that a cleaner compound might acquire a larger amount of γ and naturally a larger ρ(0).
To make a more careful comparison between theories and experiments, we now briefly discuss the doping dependence of our results. We use δ to denote the doping concentration and δ c the nematic QCP. At zero temperature, the mass of nematic field φ is proportional to the difference between δ and δ c , namely r ∼ |δ − δ c |. When the cuprate is at a distance r away from the nematic QCP, the quantum fluctuation of nematic order is not critical, but remains important for small r. At the energy scales larger than r, the fermion velocity anisotropy is still considerably enhanced, which leads C g to flow to larger values. For a given small value of C g0 , C g can flow to a sufficiently large value to induce a diffusive state and generate a finite scattering rate γ, provided that r is made sufficiently small. In this case, the quantum fluctuation of nematic order can result in a finite ρ(0) and a finite residual linear-T term of specific heat. If the bare value C g0 is large enough, random gauge potential itself suffices to generate a finite γ. Including the quantum fluctuation of nematic order leads to larger values of both C g and γ. In any case, the quantum fluctuation of nematic order tends to amplify γ, which naturally increase ρ(0) and the residual specific heat.
A number of recent experiments provided strong evidence supporting the existence of nematic order in YBCO, whereas there is little evidence for nematicity in LSCO. Although YBCO is cleaner than LSCO, the quantum fluctuation of nematic order in the former superconductor can drive C g to flow to much larger values. As a consequence, the residual linear-in-T term of specific heat of YBCO would be larger than that of LSCO. We emphasize that, disorder itself cannot explain the anomalous behaviors of the residual specific heat observed in YBCO and LSCO, and it is necessary to consider the interplay of quantum nematic fluctuation and random gauge potential.
For the above elaboration, we know that the roles played by the quantum nematic fluctuation and random gauge potential depends on doping δ. To simplify discussion, we assume the bare value C g0 displays only a weak δ-dependence. Since the quantum fluctuation of nematic order is most pronounced at the nematic QCP, the scattering rate γ and consequently the coefficient of the residual liner-T term of specific heat are maximal at the QCP, and decrease as the system moves away from this QCP. This doping dependence is observable, and can be examined by experiments. Within the loop current order scenario, both the zero-energy DOS ρ(0) and the residual linear-T term of specific heat are proportional to the order parameter of the loop current order, which decreases with the growing doping in the underdoped region.
Early experiments [67,68,69] found that the coefficient of linear-T term of specific heat at optimally doped YBCO is roughly 2 mJ·mol −1 ·K −2 . More recent measurements performed in underdoped YBCO [70] revealed that this coefficient is 1.85 ± 0.06 mJ·mol −1 ·K −2 . We feel that the currently available experimental data about the doping dependence of this coefficient are still quite limited. We expect more extensive measurements would be performed in the future to extract a more quantitative doping dependence of the coefficient, which could help to judge whether the scenario proposed in this paper works.

Strong damping rate of nodal fermions in optimally doped BSCCO
We next apply the RG results to understand the observed damping rate of nodal fermions in cuprates. Valla et al. [77] have performed extensive angle resolved photoemission spectroscopy measurements in optimally doped Bi 2 Sr 2 CaCu 2 O 8+δ (BSCCO). Their main discovery is that the nodal fermions exhibit a MFL-type damping rate in the normal state above T c , which is in general consistency with the observed linear resistivity. They further found [77] that the linear damping rate is not sensitive to the onset of superconductivity and persists well below T c . This was out of expectation since previous BCS weak coupling analysis [1,78] had predicted a quite weak damping rate, i.e., ImΣ R (ω, T ) ∝ max(ω 3 , T 3 ), in the superconducting phase. Several scenarios [1,20,79] were proposed to account for the nearly MFL behavior. In particular, Vojta et al. [20] and Khveshchenko and Paaske [79] have argued that the strong fermion damping may arise from a secondary phase transition from a d x 2 −y 2 -wave superconducting state to a d x 2 −y 2 + is or d x 2 −y 2 + id xy -wave superconducting state.
Experimentally, the existence of a nematic phase was observed in BSCCO by Lawler et al. [13]. More recent experimental studies of Fujita et al. provided further evidence pointing towards the existence of a nematic QCP in the vicinity of optimal doping in BSCCO [14,19]. Therefore, it seems natural to account for the experimental finding of Valla et al. by considering the quantum critical fluctuation of nematic order. We have showed through RG analysis that the nodal fermion damping rate caused by critical nematic fluctuation, described by ImΣ R nem (ω), is slightly weaker than that of a MFL. In realistic experiments, it is difficult to distinguish this non-Fermi liquid state from a MFL state. Therefore, presence of a nematic QCP provides an alternative scenario for the nearly MFL behavior observed in optimally doped BSCCO. However, the nearly MFL behavior occurs only at energies higher than the value ω γ set by disorder scattering rate γ. Indeed, at ω < ω γ , γ dominates over ImΣ R nem (ω), thus a finite zero-energy DOS ρ(0) is generated. To summarize, our RG results are qualitatively consistent with the quasiparticle self-energy ImΣ(ω, T ) = γ + max(ω, T ) observed in reference [77].

Temperature dependence of fermion velocity in BSCCO
Another experiment is the observation of Plumb et al. [80] that the fermion velocity v F along the nodal directions increases as T grows in BSCCO. We know from RG results that the critical nematic fluctuations drive the fermion velocities to vanish in the lowest energy limit. This means the velocities must increase as the energy scale is growing. Therefore, the nematic quantum fluctuation will induce increment of fermion velocities when the thermal energy increases with growing T , which is qualitatively well consistent with the observation of Plumb et al. Technically, one can translate the l-dependence of fermion velocity to a T -dependence by making the transformation T = T 0 e −l [31,81]. With the help of this transformation, it is easy to show that the fermion velocity v F is an increasing function of T . Therefore, the singular renormalization of fermion velocities of nodal fermions induced by the quantum fluctuation of namatic order, which is first discovered by Huh and Sachdev [24], agrees with the T -dependence of v F observed in reference [80].

Summary and discussions
In summary, we have found that the nodal fermions of d-wave superconductors constitute an unconventional non-Fermi liquid, which exhibits a weaker violation of Fermi liquid description than a MFL, due to the quantum critical fluctuation of nematic order. This unusual state represents a novel quantum state of matter that cannot be well captured by the traditional classification of (non-) Fermi liquids and thus enriches our knowledge of strong electron correlation effects. We also have calculated the fermion DOS and specific heat after incorporating the unusual renormalization of fermion velocities. When a gauge-potential-type disorder is added to the system, we have analyzed its interplay with the quantum fluctuation of nematic order, and found that the effective disorder strength flows to strong coupling, leading to diffusive motion of nodal fermions. Therefore, even a weak random gauge potential can drive a quantum phase transition from a unconventional non-Fermi liquid state to a disorder dominated diffusive state. However, the unusual fermion damping induced by the nematic order is more important than the disorder scattering at high temperatures, where the nodal fermions still display the unconventional non-Fermi liquid behaviors. We finally have discussed the connection between our theoretic results and a number of interesting experiments in the context of some cuprate superconductors.
We now would like compare our work to the existing extensive works on the non-Fermi liquid behaviors in two-dimensional metals produced by the quantum critical fluctuation of nematic order. At the QCP of Pomeranchuk instability in two-dimensional metals, the quantum fluctuation of nematic order can lead to very strong fermion damping [36,82,83,84,85,86]. To the leading order, it is found [36,82,83,84,85,86] the fermion damping rate behaves as ImΣ(ω) ∼ ω 2/3 and the quasiparticle residue Z f ∼ ω 1/3 . Since Z f vanishes in the limit ω → 0, this QCP exhibits non-Fermi liquid behavior. In this paper, we have considered the interaction between the quantum fluctuation of nematic order and massless nodal fermions in the superconducting dome of cuprate superconductors. It is apparent that the quantum critical fluctuation of nematic order gives rise to a stronger fermion damping effect in metals than in the superconducting dome of cuprates. This difference should be owing to the different forms of the kinetic energies of fermionic excitations. In the context of cuprates, the kinetic energy of the massless Dirac fermions excited from the superconducting gap nodes is E = v 2 F k 2 x + v 2 ∆ k 2 y . In contrast, in a two-dimensional metal with a finite Fermi surface, the kinetic energy of fermions can be written as E = v F k x + k 2 y 2m , where k x is the momentum component perpendicular to the Fermi surface and k y is the momentum component along the tangential direction [36,82,83,84,85,86]. In the low-energy regime, the latter kinetic energy is smaller that the former for the same given values of momenta, which indicates that the interaction plays a more important role in the latter system than the former. To further demonstrate this difference, we consider the different roles of the long-range Coulomb interaction in a two-dimensional Dirac semimetal and a two-dimensional semi-Dirac semimetal. In a Dirac semimetal, the kinetic energy of Dirac fermions is simply E = v F k with k = k 2 x + k 2 y . RG analysis showed that the residue Z f approaches a finite value at low energies, so the system is a normal Fermi liquid [52,53,54,55,56]. In a semi-Dirac semimetal, the kinetic energy of fermions is [87,88]. In this case, the long-range Coulomb interaction drives Z f to vanish in the lowest energy limit, i.e., Z f → 0, which apparently implies the breakdown of Fermi liquid behavior [87]. Once again, we see that the ratio between the interaction energy scale and the kinetic energy is a crucial quantity to determine the low-energy behaviors of an interacting fermion system.
The coupling of nodal fermions to the quantum fluctuation of antiferromagnetic (AF) order is also an interesting problem [7,89]. Uemura [7] considered the coupling of nodal fermions to the (π, π) AF fluctuation, and suggested a possibility that the (π, π) AF fluctuation can connect two nodal charges in different hole pockets and then generate a bound state of two nodal charges. More recently, Pelissetto et al. studied a number of different couplings between nodal fermions and AF fluctuations using RG method [89]. An interesting results is that, though most of these couplings are irrelevant, there emerges a nearly marginal coupling between nodal fermions and an effective, AF-order induced nematic fluctuation. This nearly marginal coupling is found [89] to results in a fermion damping rate that is nearly linear in energy or temperature.
The electronic nematic state has been observed not only in some cuprate superconductors, but also in a number of iron-based superconductors [90]: 122 family, such as hole doped Ba 1−x K x Fe 2 As 2 , electron doped Ba(Fe 1−x Co 2 ) 2 As 2 , and isovalentdoped BaFe 2 (As 1−x P x ) 2 ; 111 family, such as NaFeAs; 1111 family, such as LaFeAsO; 11 family, such as FeSe [90,91]. In most of these compounds, the nematic order emerges in accompany with a spin density wave (SDW) order. However, there are also exceptions. For instance, the nematic order is observed in FeSe without any evidence for SDW order [90,91,92,93,94,95,96,97,98]. Whether the nematic order observed in iron-based superconductors is generated by the fluctuation of SDW order or the orbital degrees of freedom is still in fierce debate [90,91,92,93,94,95,96,97,98,99]. Recent experimental studies have unambiguously showed that there is a QCP in the superconducting dome at the optimal doping of BaFe 2 (As 1−x P x ) 2 . This QCP may correspond to the critical point for a SDW order or nematic QCP [100,101], and is expected to exhibit rich quantum critical phenomena. Moreover, there are clear evidences that the superconducting gap of BaFe 2 (As 1−x P x ) 2 has nodal line points [100]. Since the quantum fluctuation of nematic order is peaked at zero momentum [31,91,102], the nodal fermions excited from the nodal line points might couple strongly to the quantum fluctuation of nematic order at the nematic QCP [31]. This coupling is physically analogous to the model considered in this work, and it would be interesting to study this coupling by means of RG method.
The RG analysis performed in this work could be generalized to study the possible non-Fermi liquid behavior and disorder effects in the context of BaFe 2 (As 1−x P x ) 2 and other iron based superconductors, where the multi-band effects and different gap symmetry need to be seriously taken into account.
which then gives rise to The scaling equation for C V is converted to (B.7) Since d lnp ∼ −dl, we find that (B.8) Substituting equations (11) and (12) into (B.8) leads us to The quantum fluctuation of nematic order drives the velocity ratio v ∆ /v F to monotonously decrease as one goes to lower and lower energies. It is known that the bare ratio v ∆0 /v F 0 < 1 in cuprate superconductors [1], which make it possible to simplify the RG equation to d ln C V d ln T = 2 + 2C 1 − C 2 − C 3 − 2C g 1 − C 1 + C 2 + C g . (B.10) In the clean limit, C g = 0, so the equation becomes This equation is slightly different from the RG equation for specific heat presented in [25], where the equation is The numerical solutions for equation (B.12) are shown in figure B1. We notice that the specific heat obtained from equation (B.12) also satisfies lim T →0 C V (T )/T 2 → ∞ and lim T →0 ln (C V (T )) / ln(T ) → 2, which indicates that the specific heat is enhanced comparing to the case of non-interacting nodal fermion system. However, as will be shown in Appendix C, it cannot be expressed as a power law function. The reason is that the three parameters C 1 , C 2 , and C 3 appearing in equation (B.12) all vanish in the lowest energy limit, rather than approaching certain finite values.

Appendix C. Approximate analytical expressions of DOS and specific heat
In order to show that the DOS and specific heat in clean limit do not exhibit power law behaviors, we now derive their approximate analytical expressions in the low-energy regime from the RG results. The RG equation of DOS is given by d ln ρ d ln ω in the absence of disorders. In the lowest energy limit, we know that the velocity ratio v ∆ /v F → 0. The parameters C 1 , C 2 , and C 3 , being functions of v ∆ /v F , all vanish [25] as v ∆ /v F → 0. Thus the RG equation can be approximately written as d ln ρ d ln ω where we have used three constants a 1 ≈ 0.426, a 2 ≈ 0.348, and a 3 ≈ 0.96 [24]. Substituting the approximate low-energy expression of v ∆ /v F given in equation (22) into (C.2), we obtain d ln ρ d ln ω which is applicable for small ω. The RG equation for specific heat is shown in equation (B.11). At low energies, it can be approximated as where a C = (2a 1 − a 2 − a 3 ) π 2 8 + ln 8 π 2 N . At N = 2, a C ≈ −1.466. Using the transformation T = T 0 e −l , we find that the specific heat at low temperature can be well approximated by the expression (C.6) By applying the same treatment, the RG equation of specific heat presented in [25], shown in equation (B.12), can be approximately expressed as where a ′ C = (4a 1 − 3a 2 − a 3 ) π 2 8 + ln 8 π 2 N . At N = 2, a ′ C ≈ −1.273. The above two functions do not exhibit power law dependence on temperature.