{'fraction_non_alphanumeric': 0.03893551728477526, 'fraction_numerical': 0.01279780055954075, 'mean_word_length': 4.24471974004874, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pseudogap state from quantum criticality

Upon application of an external tuning parameter, a magnetic state can be driven to a normal metal state at zero temperature. This phenomenon is known as quantum criticality and leads to fascinating responses in thermodynamics and transport of the compound. In the standard picture, a single quantum critical point occurs at zero temperature, which results in a nontrivial critical behaviour in its vicinity. Here we show that in two dimensions the scenario is considerably more complex due to the enormous amount of quantum fluctuations. Instead of the single point separating the antiferromagnet from the normal metal, we have discovered a broad region between these two phases where the magnetic order is destroyed but certain areas of the Fermi surface are closed by a large gap. This gap reflects the formation of a novel quantum state characterised by a superposition of d-wave superconductivity and a quadrupole-density wave, id est a state in which an electron quadrupole density spatially oscillates with a period drastically different from the one of the original spin-density wave. At moderate temperatures both orders co-exist at short distances but thermal fluctuations destroy the long-range order. Below a critical temperature the fluctuations are less essential and superconductivity becomes stable. This new phenomenon may shed some light on the origin of the mysterious pseudogap state and of the high-temperature transition into the superconducting state in cuprates. Our results demonstrate that quantum phase transitions between antiferromagnets and normal metals in layered materials may be the proper playground for search of new high temperature superconductors.

Upon application of an external tuning parameter, a magnetic state can be driven to a normal metal state at zero temperature. This phenomenon is known as quantum criticality and leads to fascinating responses in thermodynamics and transport of the compound. In the standard picture, a single quantum critical point occurs at zero temperature, which results in a nontrivial critical behaviour in its vicinity. Here we show that in two dimensions the scenario is considerably more complex due to the enormous amount of quantum fluctuations. Instead of the single point separating the antiferromagnet from the normal metal, we have discovered a broad region between these two phases where the magnetic order is destroyed but certain areas of the Fermi surface are closed by a large gap. This gap reflects the formation of a novel quantum state characterised by a superposition of d-wave superconductivity and a quadrupole-density wave, id est a state in which an electron quadrupole density spatially oscillates with a period drastically different from the one of the original spin-density wave. At moderate temperatures both orders co-exist at short distances but thermal fluctuations destroy the long-range order. Below a critical temperature the fluctuations are less essential and superconductivity becomes stable. This new phenomenon may shed some light on the origin of the mysterious pseudogap state and of the hightemperature transition into the superconducting state in cuprates. In particular we show that a chequerboard pattern reveals itself upon performing spectroscopies on the oxygen and copper atoms.
Quantum criticality has become one of the most fascinating subjects in theoretical physics for the last two decades. At a certain point in the phase diagram, when the system becomes critical, quantum fluctuations are so strong that the metallic state is broken. Electrons do not behave as stable quasi-particles anymore but become heavy and short living, hence affecting profoundly the properties of those materials in both thermodynamics and transport. Quantum critical points (QCP) are among the strongest disturbances that can be exerted upon the metallic state; as such they are top candidates to explain the mysterious behavior of high-temperature superconducting cuprates 1-3 , heavy fermions 4,5 , or doped ferromagnets 6 .
The beginning of the study of the quantum critical metals dates back to mid-seventies 7,8 . The main issue in dealing with quantum criticality in metals is that the fermions around the Fermi surface create dangerous massless modes which generate singularities in the computation of physical properties. The solution in these early works was sought in integrating the fermionic degrees of freedom out of the partition function from the very beginning. As a result, an originating low-energy Landau-like quantum action describing magnetic waves in the normal phase of the metal -called paramagnonswas used to determine the upper critical dimension and derive renormalisation group (RG) equations. This theoretical model is satisfactory for three dimensional (3D) systems. In two dimensions (2D), however, strong quantum fluctuations constitute a possibly major obstacle, capable of even invalidating the whole approach.
An important step in the study of the quantum criticality was the formulation of a spin-fermion (SF) model designed to picture 2D metals near the antiferromagnetic 9,10 and ferromagnetic 11 instability. In this model, low-energy fermions interact with a collective spin-fluctuation mode φ that close to the magnetic transition is driven into criticality.
The authors of Refs. 9,10 developed a technique similar to the one well known in the context of strongly coupled superconductivity. It consists of renormalising all the fermion propagators while neglecting the vertex corrections. Proceeding in this way, one reproduces the results of the Hertz theory 8 , though the physical picture is enriched by intriguing renormalisations [9][10][11] of Fermi velocities and vertices in 2D. As such, for a quite long time of more than ten years, the (anti)ferromagnet-normal metal transition was believed to have been understood with the critical behavior characterised by scaling relations that resemble those for second-order phase transitions.
However, as recently unveiled [12][13][14][15] , this ansatz is incomplete as the conventional Eliashberg-like theory misses some important contributions in both perturbation theory and renormalisation group thus invalidating the picture. As a result, researchers could not come to a definite conclusion about the character of the quantum phase transition and the problem has remained open.

Model and method of calculations
In this Article, we revisit the issue of quantum antiferromagnet-normal metal transitions in 2D models of itinerant electrons from the perspective of a novel general theory that shall lead us to the conclusion that the physics of the transition is considerably richer and more interesting than it has been thought so far. We demon-strate that within a slightly modified version of the SF model of Refs. 9,10 , the coupling of the bosonic spin mode to the electronic spins generates at the QCP a pseudogap in the spectrum that corresponds to an order completely different from the original spin-density wave (SDW). This new state may be understood as a superposition of dwave superconductivity and a quadrupole-density wave (QDW), and its emergence around the QCP constitutes an unexpected outcome of our theory derived from the SF model. From the high-T c cuprate perspective, the reader should remember that the SF model corresponds to the Cu sites in the CuO 2 lattice only, while the oxygen sites have been "integrated out". Actually, the quadrupole order should be pictured as induced by a corresponding modulated charge order of the four O atoms surrounding a Cu atom [ Fig. 1(a)] that, in turn, leads to an energy modulation on Cu atoms. Altogether, the modulation forms a chequerboard structure.
Interestingly, our formalism shows a certain analogy with the theory of Anderson localisation by disorder 16 . The structure of both theories relies on a summation of ladder diagrams of crossing and non-crossing subtypes, that characterise the emerging effective collective modes and their interaction. In both cases, the low-energy physics is finally captured in terms of a non-linear σmodel. Below, we present a sketch of the derivation and refer the reader to the Supplementary Material for details.
Within the SF model, the physics of electrons interacting via critical bosonic modes is described 9,10 by the Lagrangian L = L ψ + L φ with Herein, L ψ is the Lagrangian of electrons with spectrum ε(p) that propagate in the fluctuating field φ representing the bosonic spin excitations modeled by the Lagrangian L φ . The Lagrangian L φ is a quantum version of the Landau expansion in the vicinity of a phase transition.
We define the spin-wave boson mode D −1 entering Eq. (2) around the QCP through its Fourier transform where v s is the spin-wave velocity, a is a "mass" characterising the distance to the QCP [At QCP, a = 0, while a > 0 on the metallic side.], and Q is the ordering wave vector in the SDW phase. Pursuing the application of the SF model to cuprates, it is usually assumed that the spectrum ε(p) in Eq. (1) leads to a Fermi surface of the shape represented in Fig. 1(b). In order to facilitate a controlled theoretical analysis (see below), we instead consider in the following a slightly deformed Fermi surface as in Fig. 1(c). The points on the Fermi surface connected by the vector Q are the hot spots in the model and, close to criticality, the most interesting physics is formed in their vicinity. Figures 1(b) and (c) illustrate that there are eight hot spots on the Fermi surface. It is implied that the coupling constant λ is small, λ 2 ≪ vp 0 , where v is the Fermi velocity and p 0 the radius of curvature at the hot spots. 17 The quartic φ 4 -term in L φ is usually neglected for a ≥ 0. We keep the theory under control assuming that the Fermi surface has such a shape that the Fermi velocities v 1,2 of two hot spots connected by vector Q are close to being parallel to each other, with the angle δ defined in the inset of Fig. 1(c). The limit (4) favours the Fermi surface of Fig. 1(c) over the one from Fig. 1(b) and shall allow us to have all necessary approximations in our analysis under control. In the more realistic situation of Fig. 1(b), we expect qualitatively similar results but, due to the lack of a small parameter, the theory developed here would formally not be justified. The Landau damping modifies 8 the form of D −1 (ω, q), Eq. (3), adding to the latter the term γ|ω| with γ = (2λ) 2 /(πv 2 sin δ) and v = |v 1,2 |. In the limit (4), the Landau damping is strong, leading to a "weak coupling" limit of our theory, and dominates over the ω 2 -term in the bare propagator D −1 .
In the spirit of the approach of Ref. 16 , we first integrate the partition function Z = exp{− L}DφDψ over the field φ neglecting the quartic term in L φ . As a result, we obtain a model of electrons with an interaction described by the function D (ω, q) including the Landau damping term. Next, we single out those slow pairs that correspond to the mean field order parameters and derive the mean field equations. The order parameter found from these equations is strongly degenerate against SU(2) rotations, which gives rise to gapless excitations that finally are effectively described in terms of a non-linear σ-model. The properties of the latter are studied using the renormalisation group (RG) technique.

Mean field equations and pseudogap state
The mean field approximation leads to a superposition of particle-particle and particle-hole pairings, with the momentum p located at hot spots opposite to each other on the Fermi surface. In Eq. (5), σ 2 is the Pauli matrix for the electron spin and c pp p (c ph 2p ) the amplitude of the particle-particle (particle-hole) pairing. The pairings of the type in Eq. (5) are purely singlet and thus do not lead to any spin order. The signs of coefficients c pp/ph p at neighbouring hot spots on each connected piece of the Fermi surface are opposite, indicating a d-wave-like structure of the gap in the electron spectrum. Therefore, neither local charge nor current density modulations arise. The particle-hole correlations in Eq. (5) have been discussed in Ref. 14 and classified by its authors as corresponding to modulated correlations in a "valence bond solid". In that work, the authors have concluded that the superconducting correlations were stronger due to curvature effects while, at the same time, they have already noticed an emerging SU(2) symmetry in the SF model linearised near the hot spots. It is our central finding that the particle-hole and superconducting particle-particle pairings have to be considered on equal footing even for a finite curvature of the Fermi surface, and that they together form a composite SU(2) order parameter.
Let us have a closer look at the analytical form of the mean-field pseudogap in the SF model. At the end of this section, we discuss its implications for the CuO 2 lattice in the high-T c cuprates. The general solution O of the mean field equations for the order parameter at a given hot spot may be represented in the form O(ε) = b(ε)u with u being an arbitrary SU(2) unitary matrix, u + u = 1, det u = 1, and b(ε) a real positive function of the fermionic Matsubara frequency ε. After a rescaling to dimensionless quantities, ε →εΓ, b →bΓ, and T →T Γ with the characteristic energy Γ = (3λ/8) 2 π sin δ, we obtain at criticality (a = 0) a set of remarkably universal self-consistency equations that are independent of the parameters of the In these equations, sin Θ(ε) =f (ε) b2 (ε) +f 2 (ε) −1/2 .
A quick glance at Eqs. (6) reveals the trivial solutionb(ε) = 0, leading toΩ(ω) = |ω| andf (ε) = sign(ε)(|ε| + 2 π |ε|). This solution is well known as it corresponds to the one-loop self-energy corrections 9,10,19 to the bosonic and fermionic propagators. Here, of a greater interest is the existence of a nontrivial and so far unanticipated energy-dependent solutionb(ε). It can be computed numerically and its dependence onε andT is shown in Fig. 2(a). We have checked that the free energy corresponding to the nontrivial solution is lower than the one in the case of the trivial scenario withb(ε) = 0.
The characteristic value of b(ε) is of order Γ, implying that it scales linearly with the interaction constant λ 2 . This is in a sharp contrast with the exponentially small values of the gap encountered in conventional superconductors 20 . Therefore, one can expect at QCP much higher values of the gap than those obtained for non-singular interaction functions used in BCS theory. Eqs. (6) have been obtained linearising the electron spectrum near the Fermi surface and their solution formally does not depend on the position on the Fermi surface.
In fact, the order parameter O depends not only on the frequency ε but also on the distance from the hot spots, decaying at momenta of order γ/v. The characteristic length of the arc of the Fermi surface under the gap, however, should be calculated taking into account the curvature of the Fermi surface at the hot spots. As a result, the gap is large only in the vicinity of the hot spots. Schematically, b(0, p) on the Fermi surface is depicted in Fig. 2(b). At the same time, increasing the effective coupling constant -which corresponds to the, from the cuprate perspective, more realistic Fermi surface in Fig. 1(b) with the angle δ inevitably of order of unity -leads to spreading the region with a large gap over an essential part of the Fermi surface 21 . In this situation two hot spots located across the boundary of the Brillouin zone can effectively merge, such that gap b(0, p) reaches eventually the maximum value at the antinodal points, see Fig. 2(c). This effect can be enforced by in- (b) In our weak-coupling model, the gap b(0, p) is essentially non-zero only in the vicinity of hot spots. The order parameter has opposite signs at the hot spots located on the same arcs within the Brillouin zone of the Fermi surface, corresponding to a d-wave-like symmetry. (c) The gap function for the SF model on the square lattice appropriate to the cuprates [cf. Fig. 1(b)] as obtained numerically in Ref. 21 . Note that the gaps of two hot-spots adjacent to the same antinode are smeared. Beyond the weak-coupling limit λ 2 ≪ vp0, we may expect them to merge into one single gap situated at the antinode.
homogeneities and, as a result, lead to a d-wave-like dependence on the position on the Fermi surface. We refer to the state that emerges from the non-trivial solution of Eqs. (6) as a pseudogap state.
The SU(2) matrix u reflects the degeneracy of the order parameter O(ε) = b(ε)u and may be parametrised as The complex numbers ∆ + and ∆ − should be interpreted as order parameters for the superconducting and particle-hole order, respectively. In contrast to the conventional superconductivity where electron-electron pairs are formed, we have here quartets consisting of two particles and two holes, see Fig. 3 (a) and (b). Depending on the relation between the horizontal and vertical coupling, one of the pairings is more favourable but one should deal with the entire quartet when considering fluctuations. The nature of the particle-hole pairing in our theory is different from those conjectured in SU(2) theories on the basis of symmetries of t − J models 1 . Studying the symmetries of this order, we find that the d-wave structure does not lead to local charge or current modulations. However, as a consequence of the electron-hole pairing, the rotational symmetry of the electron gas is broken, giving rise to finite modulated quadrupole density 17 with ϕ + and ϕ − denoting phases. This formula describes a spatial oscillation of the off-diagonal elements of the quadrupole moment with the wave vectors Q + = (Q 1 + Q 2 )/2 and Q − = (Q 1 − Q 2 )/2, where Q 1 and Q 2 denote the vectors connecting two hot spots at ±p, cf. Fig. 1(b) and (c). Note that the vectors Q ± are considerably smaller than the SDW wave vector Q and that the resulting chequerboard structure is incommensurate with the original lattice. This new type of the particle-hole order discovered within the SF model shall in the following be referred to as a quadrupole density wave (QDW).
The appearance of a quadrupole structure becomes es-pecially intriguing if we consider it from the perspective of the CuO 2 lattice that is the playground of the physics in the high-T c cuprates. It turns out that the QDW order of the one-band SF model corresponds to a rather unique charge density wave order on the oxygen O(2p) orbitals in the CuO 2 lattice. In contrast, the average charge is the same on all Cu(3d x 2 +y 2 ) orbitals, while, however, we predict a modulation with the same wave vectors Q ± of the transition energies between Cu(3d x 2 +y 2 ) and Cu(2p) states. This effect is due to the interaction of the electrons on 2p orbitals of the Cu atoms with the modulated charges on the four neighboring O atoms.
In the Supplementary Information, we explicitly demonstrate that, as a result of the reduction of the original model for cuprates to the SF model defined on the Cu sublattice, the local electron (hole) density ρ O (r, a 0 ) on an O atom located between the two neighbouring Cu sites at r and r + a 0 transforms into the bond correlation function ψ † rσ ψ r+a0,σ . In the presence of the pseudogap obtained from Eqs. (6), we explicitly find The overall sign in this formula is + or − depending on whether a 0 is directed along an O(2p y ) or an O(2p x ) orbital, respectively. The elementary quadrupole plaquette corresponding to Eq. (9) is represented in Fig. 1(a). Thus, the order obtained within the SF model reflects a charge modulation on the oxygen sites with the same wave vectors Q ± as in Eq. (8), which can clearly be visualised as QDW. In Fig. 4 illustrating the resulting chequerboard structure on the CuO 2 lattice, the vectors R ± = 2πQ ± /|Q ± | 2 reciprocal to Q ± are drawn.
Remarkably, although the vectors Q ± are close to those connecting the antinodes, see Fig. 1(b), we did not assumed any nesting in the model under consideration. Instead, the emergence of the particle-hole QDW order is a consequence of the proximity to the QCP. Note that the period 2π/|Q ± | of the superlattice can considerably exceed the period a 0 of the original Cu lattice.
Concluding this section, we emphasize that the QDW structure is a pure charge order with no spin correlations involved. This contrasts the stripe structure proposed in Refs. 22,23 that involve both charge and spin order. While the latter does not appear in the SF model for a > 0, we do not, however, exclude its existence in the region a < 0. Another important difference between these two types of the modulations is related to the dependence of the period of the modulation on doping. The characteristic length of the stripe structures of Refs. [22][23][24] decays with the doping (as more holes means more stripes) while in our case the period of the QDW grows with the doping as the modulus of the wave vectors Q ± , see Fig. 1(b), decays.

Fluctuations and phase transition into superconducting state
The degeneracy with respect to rotations of the matrix u leads to gapless excitations. These destroy the long-range order and smear the transitions. This phenomenon is well-known in the context of the condensation of a Higgs boson and of the description of the fluctuations around ground states. Following a similar route here but taking into account the more complex symmetries of the spin-fermion model, we study the contributions of these excitations to thermodynamics by deriving the proper non-linear σ-model in a way developed in the localisation theory 16 . At T = 0, the model is 2 + 1 dimensional and the contributions of fluctuations of u as well as of the fluctuations of the amplitude of the order parameter are converging and small in the limit (4). As a result, the mean field equations (6) are applicable. Fluctuations at finite temperatures are more dangerous because the effective dimension is reduced to d = 2. The derivation of the effective free energy functional F [u] that describes the fluctuations yields In Eq. (10), t = c(T )T /λ 2 , where c(T ) is a monotonous function of the order 1 at low temperatures. The latter vanishes at the mean field transition temperature T * ∼ Γ of the pseudogap state. Further, κ ∼ γ 2 v/m with m being the effective electron mass, and R = {x/ sin(δ/2), y/ cos(δ/2)} where x and y are coordinates along the diagonals of the Brillouin zone in Fig. 1(c). Note that the coupling constant t does not contain sin δ and is small for all temperatures that are below but not too close to T * . The matrix τ 3 stands for the Pauli matrix in the space of the matrix u, Eq. (7). Its presence in the σ-model breaks the symmetry between the superconducting and QDW states. As a result, a superconducting order with a large gap ∼ Γ is more favourable at the minimum of F [u]. It is well known that fluctuations in the twodimensional non-linear σ-model at finite t and κ = 0 produce logarithmically divergent contributions destroying the long range order. This means that in our case they destroy both the superconducting and QDW orders. A finite κ serves as the infrared cutoff in the logarithms, so that the superconducting order can be stabilised below a critical temperature T c . The RG method is a standard tool for studying properties of the σ-model 25 . Integrating out step-by-step all momenta exceeding κ in the first-loop approximation, we come to the following expressions for the effective coupling constants t(κ 0 ) and κ(κ 0 ) where t 0 and κ 0 are the bare values written below Eq. (10) and Γ plays the role of the upper energy cutoff in the σ-model. Formally, Eq. (11) is applicable as long as t(κ) ≪ 1 where superconductivity is stabilised. Upon increasing the temperature, the bare coupling t and, hence, t(κ) grow, while κ 2 0 decays. These tendencies imply that fluctuations of the order parameter u become strong but, at the same time, the anisotropy between the superconducting and QDW components vanishes. This should lead for t(κ) ∼ 1 to a phase transition into a "disordered" (pseudogap) phase where no long-range order exists anymore but the gap b(ε) is still finite and remains large up to considerably high temperatures. The pseudogap state should exhibit combined properties of d-wave superconductivity and QDW, although sharp dependencies are necessarily smeared by thermal fluctuations. An estimate for the transition temperature T c between the superconducting and pseudogap phase is provided by the equation t(κ) = 1. The pseudogap can certainly be seen using STM and ARPES techniques and should still be present above the superconducting transition temperature T c roughly up to a temperature T * when the nontrivial solution of Eq. (6) is no longer possible. Of course, there is no sharp transition at T * because it is smeared by the fluctuations. It is evident that there remain strong superconducting fluctuations above T c . Note that, in contrast to the BCS theory of superconductivity, T c is related to the value of the gap Γ in quite a nontrivial way.

Phase diagram
We summarise our findings with the help of the phase diagram in Fig. 3(c). The red dashed line denotes the antiferromagnet-normal metal phase transition in the absence of interaction between the spin-wave modes and electron spins. It is the region to its right that we have studied in the present Article. We identify the pseudogap region and represent the crossover to the normal metal state by the blue line. The green superconducting region is a part of a more broad pseudogap region. Moving to the right by increasing the bosonic mass a in Eq. (3) makes the spin-fermion interaction less singular and effects of curvature of the Fermi surface more pronounced. As a result, QDW pairing is suppressed but one still obtains the d-wave superconductivity, although with a lower transition temperature. At large a, the pseudogap and superconducting states merge.
Although our results cannot be used to the left of the dashed line (a < 0), by assumption, the antiferromagnetic region exists within the SF model, Eqs. (1,2), for not very small |a|. At the same time, the pseudogap remains finite for a while after crossing the dashed line, because it would cost a large energy for it to vanish. Therefore, there must be a region of a < 0 in which the pseudogap is still finite and this region extends to the left until the phase transition (crossover) to the antiferromagnetic state sets in. At low temperatures, there is the possibility for an interesting scenario that, when going to the left, the coefficient κ 2 , Eq. (10), changes its sign. This would correspond to a transition from the superconductivity to QDW followed by a transition into the antiferromagnetic state. At the same time, we cannot exclude a direct transition from the superconductivity to the antiferromagnet and the situation will only be clarified by an explicit study of this region. We have marked this region by the question mark in Fig. 3(c).
The phase diagram represented in Fig. 3(c) looks similar to that of the cuprates (see, e.g. Refs. 1,2 ). The large values of the pseudogap can clearly be understood from the linear dependence of the pseudogap on the effective electron-electron interaction obtained in the present work.

Comparison with experiments on high-Tc cuprates
The number of the experimental works on cuprates is huge and we shall therefore restrict ourselves to the discussion of the most direct evidence from recent experiments to confirm our theoretical predictions. While the existence of the pseudogap at temperatures considerably exceeding the superconducting transition temperature T c has been discovered long ago in NMR experiments 26 followed by a variety of experimental observations including all kinds of transport and thermodynamic probes 2 , more detailed information about the structure of the pseudogap has come from ARPES measurements 27,28 . They demonstrate that the pseudogap has a d-wave-like form and competes with superconductivity. Referring to Ref. 21 , we have argued that the maximum of the pseudogap is indeed expected to move to the antinodal points, thus acquiring the d-wave symmetry.
The unique feature of the pseudogap state resulting from our analysis of the spin-fermion model is the formation of the QDW order specified by Eq. (6) and represented in Fig. 4. The wave vectors Q + and Q − of this modulated order are close to the vector connecting the antinodes.
Using STM technique, the authors of Ref. 29 have come to the conclusion about charge-modulated chequerboard patterns in several high-T c compounds with the wave vectors close to the vectors connecting the antinodes. On the other hand, a quadrupole-like order reducing the C 4 to C 2 symmetry of the lattice has been discussed in a different STM experiment 30 on BSCCO compounds. Periodic structures have been observed with the STM technique also in Ref. 31 .
Competition between superconductivity and charge modulation with and without magnetic field was studied in a very recent experiment on YBa 2 Cu 3 O 6.67 compounds using high-energy X-ray diffraction 32 . As in Ref. 29 , vectors close to those connecting the antinodes determined the periods of the charge chequerboard struc-ture. The charge structure was reported to develop above the superconducting transition temperature T c . Below T c , this structure appeared at magnetic fields suppressing the superconductivity, which is in a perfect agreement with our theory.
Resonant soft X-ray scattering technique with the photon energies near the copper L 3 absorption edge was used in Ref. 33 for studying YBa 2 Cu 3 O 6+x with hole concentrations of 0.09 to 0.13 per planar Cu atom. Again, a charge modulation has been uncovered with approximately the same wave vector. Remarkably, the intensity of the charge density wave signal is maximal at the superconducting transition temperature T c . This effect is well reflected in our theory: above T c , fluctuations destroy the charge correlations whereas below T c , the superconductivity suppresses them.
More detailed information about the modulation has been obtained in a similar experiment in Ref. 34 . The period of the superstructure corresponds to those of Ref. 33 , yet the peaks in the X-ray scattering were attributed to a spatial modulation of energies of the intra-atomic Cu(2p)-Cu(3d x 2 −y 2 ) transition. This feature is also expected from our theory. X-ray measurements at O Kabsorption edges may further clarify the features of the chequerboard order.
A magnetic-field-induced charge order in YBa 2 Cu 3 O y without any spin one has been identified in Ref. 35 with the help of NMR measurements. Although the authors interpreted their results in terms of a stripe order, the data do not seem to exclude a chequerboard one. The fact that the charge modulation appears under a magnetic field destroying the superconductivity agrees with our picture of the composite order parameter O.
We complete our discussion of experimental evidence with a very recent measurement of the sound velocity in underdoped YBa 2 Cu 3 O y in a magnetic field 36 . The authors of this work attributed their results to the existence of a two-dimensional charge order in magnetic fields exceeding the critical field B co ∼ 18 T. Interestingly, the field B co is almost temperature independent up to T 0 ∼ 40 K and then grows very sharply. This behaviour also agrees very well with our picture of the composite SU(2) order parameter with an anisotropy be-tween the particle-hole and superconducting orders that can be controlled by the magnetic field.

Outlook
Solving the spin-fermion (SF) model we have shown that the antiferromagnet-normal metal quantum phase transition in two dimensions is a complex phenomenon belonging to a unique universality class different from those known in the theory of phase transitions. Our most remarkable finding is the emergence of a pseudogap state at and near the critical line. It is characterised by an SU(2) matrix order parameter that physically corresponds to a superposition of superconducting and quadrupole density wave (QDW) orders. We emphasize that the model presented here can be fully studied in a controlled way within the framework of quantum field theory, a feature quite unique in the field of cuprate superconductivity and quantum criticality. At the same time, we note that the spin-fermion model can also efficiently be studied numerically using the Monte-Carlo method due to the possibility to avoid the sign problem here 37 .
Translating the obtained results from the SF model to the CuO 2 lattice, we have come to the conclusion that QDW state is characterised by an energy and charge modulation on the Cu and O sites, respectively. This unique feature can be directly probed with various experimental methods. The charge modulation order is evident from a number of experiments on cuprates that have been carried out recently.
The charge modulation structures observed in the recent experiments have so far not found any theoretical explanation and our theory derived from the microscopic SF is the first one providing a clear picture of what happens and how the QDW state competes with the superconductivity. A good agreement of our findings with the experimental results and their controlled derivation suggest that the theory presented here may be a good candidate for the explanation of phenomena related to high temperature superconductivity. In the earlier considerations 8,10,14 , the spin-fermion model was analysed assuming an arbitrary number N of artificial electron species or "flavours". The general idea was that the limit of large N should help to justify diagrammatic expansions. Only recently 14 , contributions have been found that are, contrarily to what one would expect from conventional 1/N -expansions, not small in the limit N → ∞. In this work, we confirm the conclusion that the large-N limit does not help to solve the model. Instead, it turns out that one can use in the calculations a different small parameter to control the theory and can thus safely put N = 1 corresponding to physical reality. Nevertheless, we keep an arbitrary N in this Supplementary Information (SI) in order to demonstrate explicitly what happens as N becomes large.
The Lagrangian L of the spin-fermion model has been written in Ref. 14 in a form that allows one to see explicitly symmetries of the system and we adopt it here, The Lagrangian (13-15) describes the fermions near the hot spots of the Fermi surface, see Fig. 1(c) of the main text, that interact via a collective bosonic mode. The 2N -component anticommuting vectors χ l α represent the electron fields with 2 being due to spin and N the number of the flavours, and the vectors v l α are Fermi velocities at the hot spots. Herein, the superscript l numerates the pairs of the hot spots connected by the vector Q =(π, π) of the spindensity wave, while the subscript α distinguishes the single hot spots inside each pair. We choose coordinates along the diagonals of the Brillouin zone. The three-component vector field φ describes the bosonic spin-wave modes with velocity v s that couple to the electronic spin, represented in Eq. (14) by the operator σ which is the vector of Pauli matrices. The parameter a(T ), the bosonic "mass", determines the distance from the critical line in the phase diagram and vanishes at criticality. In the spin-fermion model, we consider fields χ and φ depending on (1+2)-dimensional coordinates X = (τ, r) combining imaginary time τ , which varies between 0 and the inverse temperature β = 1/T , and coordinates r for the effectively two spatial degrees of freedom. The parameter a(T ) determines the distance from the critical line of the phase diagram.
The spectrum of the fermions written in Lagrangian L f , Eq. (13), has been linearised. From the Fermi velocities v l=1 in the hot-spot pair l = 1, we obtain all other v l α by inversion and rotation by π/2. The linearised form is not always sufficient for the computation of physical quantities and a slightly more general Lagrangian taking into account curvature effects may become necessary. We shall include the curvature corrections shortly.
For the purpose of clear and compact notations, we introduce 32N -component fields Ψ, Herein, the label τ stands for the particle-hole space and χ is a 16N -component field, The superscript L = 1, 2 distinguishes the pairs of hot spots on the Fermi surface that are connected by Q (for L = 1) from those pairs connected by Q ′ = (−π, π) (for L = 2), cf. Fig. 1(c). Finally, the components of are the spinor fields χ l α appearing in the original formulation (14). In these equations, * maps the anticommuting variable χ onto its conjugate and t denotes transposition. We use the convention that (η i η j ) * = η * i η * j and (η * i ) * = −η i for arbitrary anticommuting variables η i . The symbols Σ and Λ refer to the pseudo-spin spaces corresponding to the subscript α (distinguishing hot spots within one critical pair) and the superscript l (distinguishing the two pairs of hot spots), respectively. Thus, for N = 1, the field Ψ is a tensor in the 32-dimensional space τ ⊗ L ⊗ Λ ⊗ Σ ⊗ σ with σ being the physical spin. In the following, we shall use besides σ the Pauli matrices τ , Λ, and Σ in the corresponding factor spaces. At the same time, all relevant physical quantities are diagonal in L-space.
Besides the hermitian conjugate Ψ + of the vector Ψ, it is convenient to introduce in addition a "charge" conjugate fieldΨ defined asΨ The matrix C satisfies the relations C t C = 1 and C = C t . It is clear that The notion of charge conjugation is naturally extended to arbitrary matrices M (X, X ′ ) as It is easy to see thatM We shall refer to matrices satisfying the relationM as anti-selfconjugated. Such matrices frequently appear when dealing with quadratic forms for anti-commuting vectors.
All the notations (17)(18)(19)(20)(21)(22)(23) are very similar to those used in Ref 16 . Using the multicomponent field Ψ and including also the curvature corrections into the spectrum of the fermions, we write the action S f that corresponds to the Lagrangian L f , Eq. (14), as with Herein, m is the effective fermion mass at the hot spots and the velocity operatorV, written as a matrix acting in the space of χ, Eq. (18), has the form with In H curv only the momentum component directed along the Fermi surface is kept, which is sufficient in the weak coupling limit where v and p 0 are respectively, the modulus of the velocity and the radius of the curvature of the Fermi surface at the hot spots. The velocitiesṽ ± are mapped onto v ± by rotation over the angle π/2. All parts of the Hamiltonian, H 0 , H curv , and λΣ 1 σ t φ (X), are anti-selfconjugate, an immediate consequence of Generally, the action is invariant under transformations where U is a 32N × 32N unitary matrix with unity blocks in spin and particle-hole spaces, while the diagonal blocks in the L ⊗ Λ ⊗ Σ-space, the space in which the matrixV, Eq. (27), acts nontrivially, are arbitrary. In other words, U consists of eight independent U(N ) matrices, reflecting the independent particle number conservation at each single hot spot. If we neglect the term H curv containing τ 3 in the action S f , the latter becomes invariant under arbitrary particle-hole rotations, thus enlarging the symmetry group for transformations of the type (31) to U(2N ). This U(2N ) symmetry of the linearised action is crucial for physical properties of the model and we shall see that its spontaneous breaking in the vicinity of the QCP leads to a rather intricate order parameter.

A. Self-consistency equation
We computate the partition function first by integrating over the bosonic field φ. This is possible in a closed form only if we neglect the quartic interaction in L φ , Eq. (15), which constitutes our basic approximation. In conventional theories of phase transitions it is not legitimate to neglect the quartic term since perturbation theory in g breaks down at the critical point at and below the upper critical dimension. In the present model, one encounters logarithmic divergencies in the perturbation theory at zero temperature 8 . However, renormalisation group (RG) studies of the φ 4 theory show a logarithmic decay of the coupling constant g and therefore following Ref. 8 the quartic term is usually neglected. The situation is more tricky at finite temperatures because the φ 4 model does not have a phase transition in two dimensions, while neglecting the quartic term the system would undergo the phase transition. This problem is discussed later with a conclusion that the quartic term can still be neglected at the cost of assuming that effectively the parameter a (T ) should always remain finite at finite T but vanishing in the limit T → 0. Such a behaviour of a (T ) should mimic the smearing of the phase transition in the absence of the quartic term.
Then, integrating out the bosonic field φ yields where the noninteracting terms for the linearised kinetic spectrum and curvature corrections read The fermion-fermion interaction includes as interaction potential the original bare bosonic propagator ω = 2πT m, m = 0, ±1, ±2, . . ., are Matsubara bosonic frequencies. Note that in SI we measure the electron momenta from the hot spots. As a result, the propagator D (ω, q) in Eq. (37) has formally been shifted by the vector Q with respect to the one in Eq.
In principle, one can study the model defined by Eqs. (33)(34)(35)(36)(37) using an expansion in S int . The bare Green function G 0 for the Hamiltonian H 0 used in this type of perturbation theory is written in Fourier space as with ε = π (2n + 1) T , n = 0, ±1, ±2, . . ., denoting a fermionic Matsubara frequency. This approach has been used in the previous publications 9, 10,14 . It has been found in Ref. 14 that anomalous (in terms of the expansion in 1/N ) contributions come from particle-particle and particle-hole ladder diagrams of the type represented in Fig. 5 These diagrams describe superconducting and some kind of insulating fluctuations demonstrating a tendency to a corresponding particle-particle and particle-hole pairing. Formally, these diagrams resemble "cooperons" and "diffusons" appearing in the localisation theory 16 . More complicated diagrams correspond to an interaction between these effective modes.
A very efficient way to sum the contributions of all these diagrams is to derive a non-linear σ-model 16 and study fluctuations of an "order parameter" using this effective field theory. Below, we follow a similar route. We investigate the model by first writing self-consistent mean field equations, then solving them, and finally studying fluctuations. Of course, one can describe the fluctuations within a perturbation scheme again. However, now the expansion will be performed near another minimum. This is a standard situation in models where a symmetry of the original Hamiltonian is broken in a certain region of parameters.
The first step consists of replacing the Ψ 4 -interaction by a quadratic term with coefficients to be determined in a self-consistent way: The effective action for the fermion-fermion interaction shall contain a bosonic propagator D eff (X − X ′ ) -physically the spin susceptibility -that is renormalised by particle-hole bubbles, see Fig. 6, that in turn are composed by renormalised fermions. In order words, we develop a self-consistency scheme for both fermions and bosons.
The renormalised propagator D eff (ω, q) can be written as In Eq. (40), G (ε, p) is the Green function of the renormalised fermions and Tr is the trace of the 32N × 32N matrices. In order to write equations for G (ε, p), we replace the interaction term S int [Ψ] by a quadratic form, In the spirit of the mean field approximation, the matrix M (X, X ′ ) should be found in a self-consistent way from the equation where the symbol . . . eff is the mean field average, Herein, the effective action S eff [Ψ] reads with the operatorM (X) acting on the field Ψ (X) aŝ The factor of 2 in Eq. (43) arises due to the two possible pairings of "Fock-" and "Cooper"-type that are in fact equal to each other. Now, we are in a position to define the Green function G(X, X ′ ) as the propagator for the action S eff , Eq. (44), It satisfies the equation Using G(X, X ′ ), the self-consistency equation (43) takes the form Equation (48) is the basic equation that determines all the unusual physics at and near QCP. It is an analog of the BCS equation in the conventional theory of superconductivity 20 . However, the singular form of the interaction D eff (X − X ′ ) leads to a more complex mean field solution and, as a result, to richer physics than that given by the BCS theory.

B. Solution of the self-consistency equation
We look for the solution M (X, X ′ ) to Eqs. (47,48) using the ansatz which is Fourier transformed as The self-consistency equation (48) thus becomes This Green function should also be used in Eqs. (40) for the propagator D eff so that the system of the self-consistency equations is closed. Note thatM as a consequence of the anti-selfconjugation property of M , Eq. (49). Eqs. (52, 53) can also be obtained by summation of all non-intersecting diagrams for the self energy like those depicted in Fig. 7.
The presence of the spin matrices σ in Eq. (48) prevents the solution from being invariant with respect to rotations in the spin space. Parametrising the spin structure of M 0 (ε, p) as where Q(ε, p) is the singlet and S (ε, p) the triplet component, we can try solutions of a pure singlet form, Q(ε, p) = 0 while S(ε, p) = 0, or of a pure triplet form, Q(ε, p) = 0 while S (ε, p) = 0. At the same time, the symmetry of M implies thatQ In the first case of a pure spin singlet, the self-consistency reduces to while in the triplet case, we obtain Both the equations (57) and (58) This implies that the Q L (ε, p) are independent of each other, reflecting the independence of the hot-spot pairings in perpendicular directions in the Brillouin zone in Fig. 1(c). For each direction L, we seek for the singlet solution Q L (ε, p) using the representation where 1 1 denotes the unity matrix. The matrixB L (ε, p) is unity in all spaces except Λ-and particle-hole τ -space, a(ε, p) and w(ε, p) are functions without any matrix structure. By construction, these two functions and the matrix B L (ε, p) satisfy the parity relations Further, we assume that the matricesB L (ε, p) anticommute with the matrix Λ 3 , The properties (61, 62) for the matrixB(ε, p) allow us to write it in the form Herein, u L denotes a 2N × 2N unitary matrix rotating particle-hole and flavour components and b (ε, p) is a real even function, b(ε, p) = b(−ε, −p). The matrix U L is unitary, U L+ U L = 1, and also satisfiesŪ L U L = 1, implying for the matrix u the symmetryū As a result, det u L = (det u L ) * so that u L belongs effectively to the SU(2N ) group. For N = 1, we may apply the parametrisation [Eq. (7) of the Article] Recalling the structure of the vectors Ψ, Eq. (17), we have concluded in the Article that ∆ L + plays the role of the superconducting order parameter, whereas ∆ L − is the order parameter for an insulating electron-hole pairing order. Equivalently, we may employ a parametrisation in terms of a 4-component unit vector, The angles θ L determine to what extent the ordering is superconducting and to what extent insulating. At θ L = 0, only the insulating order exists in the L-direction, while at θ L = π/2 the subsystem should manifest superconducting properties. The angles χ L and ϕ L are the phases of the insulating and superconducting order parameters, respectively. In the linear spectrum approximation used so far, the directions of the vectors ∆ L are arbitrary. Using Eqs. (60, 63), we reduce the self-consistency equation (57) to a system of three equations for a (ε, p), w (ε, p), and b (ε, p), where f (ε, p) = ε + a (ε, p) signε and v stands for the Fermi velocity at a hot spot. In the approximation used below the solution weakly depends on the momentum and we omit in the following the momentum in the arguments writing f (ε), w(ε), and b(ε).
Including the mean field in the fermion propagators, the polarisation bubble Π (ω, q), Eq. (40) renormalising the bosonic propagator D eff (ε − ε ′ , p − p ′ ), Eq. (39), takes the form The velocities v 1 and v 2 may be chosen to be those given by Eq. (16) but the result does not depend on the choice of the hot-spot pair l. Also, anticipating the limit to be discussed below, we have neglected the function w(ε). Since v 1 and v 2 are linearly independent, the projections of p on v 1 and v 2 are a convenient choice for the two independent variables. As a result, we immediately see that Π depends only on ω and, neglecting the irrelevant ω 2 -term of the bare bosonic propagator, we find Herein, the dynamic part is characterised by the function Ω(ω) where the coupling constant is given by and δ is the angle between the two Fermi velocities in one hot-spot pair, v l=1 see also the inset of Fig. 1(c). All the approximations made throughout our calculations are justified assuming a small but nonzero angle, while N can be arbitrary, including the physical value of N = 1.
In the absence of the gap b(ε), the function Ω(ω) has the form of the well-known Landau damping 8 , The existence of the gap reduces the damping which in turn pushes the interaction D eff to become stronger. Equations (67-69, 71-73) give us the closed system of equations for the physical quantities we are interested in. We can further simplify Eqs. (67-69) assuming that the important contributions in the integrals come from the region of momenta |p ′ ⊥ | ≪ |p ′ |, where p ′ and p ′ ⊥ are the components of the vector p parallel and perpendicular to the Fermi surface. The inequality (75) guarantees that this is indeed the case. Moreover, this limits allows to neglect the function w(ε, p).
Neglecting |p ′ ⊥ | with respect to |p ′ | in the function D eff (ε − ε ′ , p − p ′ ), we integrate this function separately over p ′ and the rest of the integrand over p ′ ⊥ . As a result, the solutions f and b no longer depend on the momenta and we simplify Eqs. (68, 69) to We see that the typical values of the functions f (ε) and b(ε) are of order Γ and thus important momenta p ⊥ in Eqs.  8: (a) The dimensionless quantitiesb andf as functions of the reduced frequencyε at temperatureT = 0.001 at the QCP (a = 0). Under the peak of the gap functionb, the dynamic term in the fermionic propagator is linear, while the interaction with the Landau-damped bosons leads to the characteristic square root law only at larger frequenciesε ∼ 1. Inset: The same but far away from the criticality (ā = 1). The gapb is considerably suppressed and the Landau damping is ineffective, resulting in the linear frequency dependence as in the case of free fermions. (b) The dynamic termΩ of the boson propagator. At zero bosonic mass, the fermionic gap b suppresses the Landau damping and Ω ∝ ω 2 , while at a large massā = 1 the fermion gapb is small and the |ω|-law is restored.
neglect of the quartic term in the Lagrangian L φ , Eq. (15), may be justified at finite temperatures by assuming that the function a (T ) remains finite at the QCP but tends to zero in the limit T → 0.
One can understand this statement considering the first order correction to the coupling constant coming from the zero Matsubara frequency. In the first order, the zero frequency renormalisation takes the form where c is a numerical coefficient. We can neglect this contribution only if a ≫ T g. At the same time, we should keep in mind that at finite temperatures, the antiferromagnetic transition in 2D is smeared because the thermal fluctuations destroy the antiferromagnetic order. The smearing of the transition means that effectively the "mass" a (T ) cannot turn to zero at critical point at T = 0 because the latter does not exist. In order words, considering fluctuations of the static component of φ, we have a cutoff at k 0 ∼ √ T . In order to avoid all these complications, we simply drop the term withε ′ =ε from the sum overε ′ in Eqs. (81, 82) when studying these equations numerically. This is how our result in Fig. 2(a) has been obtained. Neglecting this term can lead to a somewhat lower mean field transition temperature into the pseudogap state T 0 but its precise value is not very important for our present discussion. Anyway, at low temperatures this should be a good approximation.
The numerical solution of Eqs. (81-83) is also represented atT = 0.001 in Fig. 8. In Fig. 8(a) one can see a finite functionb (ε) leading to a gap in the fermionic spectrum. At a finite a, Fig. 8(b) shows the Landau damping [linear dependence ofΩ (ω) on |ω|] whereas at a = 0 the dependence ofΩ (ω) on ω is quadratic. The latter is a consequence of the existence of the gap in the spectrum of the fermions. A more general picture is given by Fig. 2(a).
Actually, theω 2 asymptotics ofΩ (ω) at smallω leads to a logarithmic divergence of b (ε) in the limit ε → 0, which is clear for Eqs. (81,82). However, this divergency is not important for our consideration. Moreover, this logarithmic divergence may be cut off by taking into account the renormalisation of the coupling constant g and fermionic Green functions by non-zero Matsubara frequencies. This fact can be understood writing at low temperatures instead of Eq. (84) the following renormalisation As Ω (ω) ∝ ω 2 at small ω, one can understand from Eq. (84) that the renormalisation of g is equivalent to the one of a classical three dimensional φ 4 -model. Scaling relations and exponents are very well known for this model. The bosonic propagator D (ω, k, a) can be written at a = 0 in a scaling form 25 where α (x) is a function of x, with a finite α (0) and α (x) ∝ x 2−η , x → ∞ Substituting Eq. (86) into Eqs. (77-79) we can see that the logarithmic divergence of the solution b (ε) is cut because the integral over ε ′ converges. As the exponent η is numerically small, the solution and all subsequent quantities are very close to those obtained for η = 0 when the quartic term in the Lagrangian L φ , Eq. (15), was neglected. As concerns the triplet order parameter S, Eq. (55), we seek a solution S in the form whereē i (ε) = e i (−ε), z + z = 1, andz = z + . This leads us to an equation for e i , i = 1, 2, 3, In Eq. (88) the propagatorD eff is a 3 × 3 matrix. This equation differs from Eq. (68) by the absence of the factor 3 in front of λ 2 , which implies that the triplet solution is energetically less favorable. Therefore, we do not consider the triplet pairing anymore.
In conclusion, the mean field treatment of the interaction between fermions and critical spin modes leads to a spontaneous breaking of the symmetry and formation of a gap. The emerging order is a mixture of a singlet dwave-like superconductivity and a particle-hole insulating order. Of course, we cannot fully trust mean field theories unless we have demonstrated that the mean field state corresponds to a the minimum and we have estimated the contributions coming from fluctuations. We will do this in the next Section, while in the remaining of this one we discuss symmetry properties of the new state and study the effect due to the actually nonzero curvature of the Fermi surface.

C. Symmetry of the order parameter
The off-diagonal part Q off (in Λ-space) of the matrix Q, Eqs. (60, 63), is the matrix order parameter containing the superconducting (∆ + ) and an insulating particle-hole pairing (∆ − ), cf. Eq. (65). This order arises as a result of the pairing of particles and/or holes from hot spots that are located opposite to each other on the Fermi surface. At the same time, the presence of Σ 3 in Q off means that the order parameter changes sign when passing from a hot spot to its partner connected by the vector Q, see Fig. 1(c). Such an oscillation of the sign corresponds to a d-wave-like structure of the wave functions. Recalling that the order parameter Q originates from the spin singlet pairing, we thus conclude that the superconducting part of the order parameter describes d-wave superconductivity. A full body of experimental methods for its identification exists and superconductivity is a phenomenon that can hardly be overlooked in transport or thermodynamic measurements. In contrast, the particle-hole part of this order is somehow "hidden". As it is singlet with respect to spin, no spin order appears. In the following, we are going to find what observable correlations arise due to the the component ∆ − of Q.
Absence of charge density order. Let us begin this study by directly computing the local density modulationρ (r) due to the electron-hole ordering. This quantity can be written in terms of the averagẽ where ψ * p,σ and ψ p+Q L ,σ are fermionic fields at equal times in the conventional formulation [p measured from the Γ-point] and σ is the spin index with summation implied for repeated indices. The vectors Q L with L = 1, 2 connect opposite hot spots ±p, cf. Fig. 1(c) and also Ref. 14 , and are directed along the diagonals of the Brillouin zone. Using the fields χ l a and the definition of the vectors Ψ, Eqs. (17)(18)(19), we rewrite Eq. (89) as where . . . eff stands for the averaging with the action S eff , Eq. (45), the momentum k is counted from the hot spots, and, by definition, Q l = −Q −l .
Using the Green function G, Eq. (46), we obtaiñ and U L is given by Eq. (63). The momentum S is the length of the arc on the Fermi surface covered by the gap. At first glance, the integral determining S diverges at large p . This is a consequence of neglecting the curvature of the Fermi surface. In the next subsection devoted to the effects of curvature we shall estimate the integral S as A quick glance at the second line of Eq. (91) reveals however that the trace yields zero anyway due to the presence of the matrix Σ 3 . Therefore, the electron-hole pairing considered here does not lead to any local charge density modulationρ (r) = 0 (94) In fact, this reflects the d-wave character of the pairing, leading to an oscillation of the order parameter when moving along the Fermi surface.
Absence of orbital currents. As concerns an average orbital current j, we start with a formula analogous to Eq. (89), where v = ∂ε(p)/∂p =(v x (p), v(p) y ) is the velocity. Expressing the current in terms of the vector field Ψ as in Eq. (90) and then in terms of the order parameter U , we find The trace in Eq. (96) vanishes due to the presence of the matrix Λ 3 inV, cf. Eq. (27). The presence of the matrix Σ 3 under the trace would not be sufficient for the trace to vanish becauseV also contains this matrix.
Absence of dipole order. Now let us check whether an average dipole moment d arises. We write this quantity as or so that no average dipole moments exist. Since no macroscopic scalar and vector order parameters can be constructed for the electron-hole pairing, we now try to compose a tensor order parameter. The natural choice is to consider the 2 × 2 quadrupole tensor.
Quadrupole order. In the model considered here, the 2 × 2 quadrupole moment tensorD (r) can be composed aŝ where the traceless 2 × 2 matrix operatorÎ p equalŝ and e is the electron charge. Proceeding in the same way as previously, we writê Now we can calculate the derivatives with respect to momentum p. It is clear that when neglecting the curvature only the derivatives along the vectorV are not equal to zero. This allows us to write the quadrupole tensor D (r) in the form where i, j ∈ {x, y}. The evaluation of the integral yields where Computing the trace in Eq. (104) and using Eqs. (80,93), we find the quadrupole moment tensor in the form D (r) =D (r) 0 1 1 0 with (106) This formula shows that the electron-hole pairing order corresponds to a quadrupole moment oscillating in space in the directions Q 1 and Q 2 , cf. Fig. 1(c). This results in a chequerboard structure as illustrated in Fig. 4, which is a quadrupole density wave (QDW). The vectors Q 1,2 are directed along the diagonals of the Brillouin zone but their length is considerably smaller than that of Q = (π, π). The tensorD is diagonal for coordinate axes parallel to the boundaries of the Brillouin zone. Remarkably, the quadrupole density in Eq. (107) does not depend on the strength of the spin-fermion interaction. The quadrupole moment essentially depends on the angle θ, Eq. (66), and vanishes at θ = π/2, which corresponds to a pure d-wave superconducting state.
Valence bond order. Nonzero off-site correlations have already been discussed in Ref. 14 . Here, we explicitly compute the real part of the correlation function ψ * rσ ψ r+a0,σ in the presence of the pseudogap. Seeking the application to cuprates, we think of r as a site in the Cu sublattice and of a 0 as a primitive vector of this lattice pointing either in horizontal or vertical direction in Fig. 1(a). In line with Eq. (89), we write ψ * rσ ψ r+a0,σ + c.c. = ψ * p−Q l ,σ ψ pσ + ψ * pσ ψ p+Q l ,σ exp ipa 0 + iQ l r dp (2π) 2 . (108) In the linearised SF model valid around the hot spots, we can re-express the right-hand side in terms of the 32component field Ψ, Eq. (17), and the momentum k is from now on counted from the hot spots, whose location is determined by ± (Q L /2 + (π, 0)) and ± (Q L /2 + (0, π)). This is the origine of the shift of the coordinate r in the exponent in Eq. (109). In the mean-field approximation for the pseudogap state, the average in Eq. (109) is to be replaced by the Green function of Eq. (46), In the coordinates of Fig. 1(c), a 0 = (a 0 / √ 2)(1, ∓1) with the sign in the second vector component depending on whether a 0 points along a horizontal Cu-Cu bond (−) or a vertical one (+). The parameter a 0 = |a 0 | is the Cu lattice constant. For the projection of the momentum vector k = k ⊥ + k onto the subspace perpendicular to the Fermi surface, we may thus write We observe that Σ 3 also appears in the exponent, and as a result the trace in Σ-space does not vanish in the presence of the pseudogap. The integration over the momentum k extends over the length S, Eq. (93), along the Fermi surface. As this length is small, Sa 0 ≪ 1, we may neglect k a 0 in the exponent. The remaining integration over k ⊥ is straightforward, yielding The function Y 3 (T ) is defined as The contributions that are traceless in Σ-space have already been omitted and the energy Γ has been introduced in Eq. (80). Calculating the trace in Eq. (112) we finally obtain The modulation of the off-site correlation function corresponds to the modulation of the on-site quadrupolar order, Eq. (107), and is as the latter proportional to cos θ = |∆ − |, cf. Eq. (66). Note that the different overall sign for horizontal (−) and vertical (+) bonds reflects the quadrupolar character of the particle-hole component in the pseudogap order.

D. Effect of curvature
The mean field equations (81-83) have been derived neglecting the curvature of the Fermi surface. As a result, the solution Q (ε) does not depend on the position on the Fermi surface and one might come to the conclusion that the entire Fermi surface is covered by the gap. Of course, this is not the case and the Fermi surface is covered by the gap only near the hot spots where the spectrum has been linearised. In this Section, we estimate the length of the gapped part of the Fermi surface around a hot spot. Figure 9 illustrates the geometry for one of the four pairs of hot spots in the Brillouin zone [ Fig. 1(c)]. Let us reconsider Eq. (68) focussing on this pair. The two hot spots are connected by the vector Q = (π, π). In terms of the original fermionic operators ψ p,α with p close to the Fermi surface, we write the equation for the gap in the form How large is the region near the hot spots where it is justified to consider the gap b (ε, p) as independent of the momenta p along the Fermi surface? To answer this question, let us assume that b (ε ′ , p ′ ) in the integrand on the right-hand side does not depend on p ′ . Furthermore, assuming that typical values of the momentum components p ⊥ perpendicular to the Fermi surface are much smaller than parallel components p (to be verified a posteriori) we may neglect p ⊥ in the bosonic propagator D eff and perform the p ⊥ -integral. Finally, we have to investigate the integral and find out when it does not depend on p.
Taking into account the curvature close to the hot spots, we approximate the Fermi surface in their vicinity by circular arcs as shown in Fig. 9. The radius p 0 of these arcs has to be chosen in such a way that the curvature of the arc equals the radius of the curvature of the actual Fermi surface at the hot spot. When integrating over p ′ ⊥ we neglect any non-linear dependence of the fermionic spectrum on this variable. This is justified in the weak coupling limit (30). The centers of the circles for the two hot spots of the pair in Fig. 9 are located at a distance 2P 0 from each other. Let us choose coordinates according to Fig. 9 with the center of the coordinate system in the center of the Brillouin zone. Then, we express the vector Q and the two points p and p ′ on the Fermi surface as Q = 2p 0 P 0 /p 0 − sin(δ/2) , 0 , p = p 0 P 0 /p 0 − sin(δ/2 + α) , − cos(δ/2 + α) , The angle δ is the angle between the Fermi velocities v 1 and v 2 introduced previously in Fig. 1(c). With these notations, we have close to the hot spots (|α| , |α| ′ ≪ 1) Since p ≃ p 0 α and p ′ ≃ p 0 α ′ , the integral (116) takes the form . two. It is clear that the solution b is independent of p so long as this last term in the denominator can be neglected. Hence, we estimate the length S of the gapped region as As we have seen in the previous subsections, typical values of Ω(ε − ε ′ ) are of order Γ, Eq. (80). Thus, at the critical point, a = 0, we come to the estimate already given in Eq. (93), We have assumed throughout the mean field scheme that the essential |p ⊥ | are much smaller than p . As |p ⊥ | ∼ Γ/v and p ∼ S, we obtain which justifies the approximation. Making the parameter a finite increases the area S under the gap until the entire Fermi surface gets covered by the gap at sufficiently large a. For the length of the gapped region S (a) at finite a, we find where S (0) is given by Eq. (93). This is clearly not more than a qualitative estimate but it shows the increase of the gapped region. The value of the gap b itself decreases when increasing a, though. All formulas obtained in the mean field approximations could also be justified in the limit of large N for any sin δ. However, this limit is not useful in establishing of the validity of a σ-model considered in the next Section and cannot substitute the inequality (75).
The bottom line of this subsection is that, whenever integrals dp appear, they should be simply replaced by the value S, Eq. (93).

III. FREE ENERGY AND NON-LINEAR σ-MODEL
A. Free energy in the mean field approximation All the results of the previous Section have been obtained in the mean field approximation. Now, we would like to calculate the energy of the ground state and consider fluctuations. This is certainly important for checking the applicability of the mean field scheme but, at the same time, fluctuations enrich the physics and lead to new phase transitions. Therefore, after having established the mean field picture, the investigation of their role is the next step in the study of the spin-fermion model. However, let us first concentrate on calculating the free energy.
An important peculiarity of the mean field approximation is that we not only renormalise self-consistently the fermionic Green functions, Fig. 7, but also the bosonic propagator D is renormalised by including polarisation bubbles, Fig. 6, which in turn contains the renormalised fermionic Green functions. In principle, we would encounter the pseudogap state also just using the Landau damped form 8,9,10,14 of the propagator D, i.e. inserting the bare Green functions into Eq. (40). However, this approximation would be incorrect now and, even though a finite solution for b (ε) exists in this case, the values of the gap b (ε) would be considerably smaller. As we have seen, including the self-consistent renormalisation of the fermion and boson modes into the mean field scheme is not difficult but the calculation of the ground state energy and the derivation of a field theory for the fluctuations requires a certain care.
The most useful feature of the diagrammatic technique is the possibility to replace parts of diagrams by "blocks" (see, e.g. Ref. 25 ). In other words, one can calculate vertices or self-energies within a certain approximation and eventually work with these quantities rather than with sums of diagram classes. While such a scheme is possible for the calculation of Green functions, it does not apply as directly for computation of the free energy. This problem can be overcome by considering a derivative of the free energy or of ln Z with Z being the partition function, Eq. (33). Here, let us therefore write its dependence on the coupling constant λ explicitly. We have Then, the average of S int [Ψ] entering Eq. (123) is calculated using the block structure of our theory. The partition function Z λ is recovered at the end as Within the mean field approach of the previous Section, we have to insert the Green function G λ , Eq. (53), and the propagator D λ eff , Eq. (39, 40), both depending on λ, into ln Z λ . In the lowest order in the perturbation theory, we have for the derivative entering Eq. (124) where the bare λ-independent bosonic and fermionic propagators D and G 0 have been defined in Eq. (37, 38). The self-consistency scheme is equivalent to the summation of non-intersecting diagrams for the self-energy of fermions, Fig. 7, and polarisation loops consisting of the dressed fermion lines for the self-energy of the bosons, Fig. 6. As a result of such a summation, one comes instead of Eq. (125) to where the Green functions D λ eff (ω, q) and G (ε, p) are given by Eqs. (39, 40) and (53, 60), F λ is the free energy per unit volume, and V is the volume of the system.
Using the fact that the spin blocks of the Green functions G are proportional to unit matrices we obtain Using Eqs. (53, 60, 63) and subtracting the derivative of the free energy for the solution f 0 (ε) of Eq. (81) with b(ε) ≡ 0, we obtain for the difference of the free energies ∆F λ where E (ε) = f 2 (ε) + b 2 (ε) and sin Θ = f (ε) /E (ε). Now, we are in the position to write down a formula for the difference ∆F between the free energy in the case of the solution with a gap b (ε) = 0 of Eqs. (81, 82) and the free energy for the gapless solution [b (ε) = 0]. Measuring all quantities with the dimension of energy that enter this formula in units of Γ, Eq. (80), we find Figure 10 displays ∆F as a function of the reduced temperatureT and the bosonic mass a. This function is clearly negative, which shows that the solution with b = 0 is energetically more favourable than the one without the gap [b(ε) ≡ 0]. Although the function ∆F vanishes in the limit T → 0 the function F/T logarithmically diverges. As the latter enters the weight of the states, we conclude that the pseudogap does exist below the critical temperature T 0 .

B. Non-linear σ-model
The mean field solution Q (ε), Eq. (60, 67, 68), does not depend on the coordinates. Once we know Q (ε), we can write the free energy F in this approximation. At the same time, the solution Q (ε) is infinitely degenerate because of the arbitrariness of the SU (2N )-matrix U in Eq. (63). As a result of the spontaneous breaking of the symmetry, gapless Goldstone modes exist in the system and strongly affect properties of the 2D model under study. The mean field theory developed so far does not take into account these fluctuations and is thus insufficient for describing important effects. Therefore, we have to go beyond the mean field theory and derive a field theory that finally describes the gapless excitations.
A similar problem is encountered in the Anderson localisation theory and can be successfully treated with the help of a supermatrix non-linear σ-model. The mean field theory for the disorder problem corresponds to the self-consistent Born approximation but the latter is not sufficient for studying the localisation and this is why one should study fluctuations with the help of the σ-model. The supersymmetry is needed in order to be able to average over the quenched disorder in the very beginning of the calculations before doing any approximations.
For the spin-fermion model considered here the integration over the field φ can be done exactly [neglecting the φ 4 term in L φ , Eq. (15)], leading to the Ψ 4 interaction, Eq. (35). Following the scheme developed in the localisation theory, we would single out slow pairs of the the fields Ψ and decouple them by integration over an auxiliary matrix field Q (Hubbard-Stratonovich transformation). Then, integrating over the field Ψ, we would obtain a functional of the field Q, for which the next step would be to find the saddle point equations. These equations coincide with the mean field equations. In the final step, one should consider fluctuations near the saddle point and derive the non-linear σ-model.
In principle, we could follow this route if we did not have to take into account the gap in the spectrum of the fermionic Green functions when calculating the polarisation loops that renormalise the bosonic propagator D eff , Eq. (39, 40). The problem is that the gap should be found self-consistently, so that we cannot choose the self-consistent propagator D eff as the interaction in the original model. Since there is no simple way to include the self-consistent renormalisation of D eff into the scheme based on the Hubbard-Stratonovich transformation, we have to proceed in a somewhat different way.
We assume here that the entire volume in {X}-space is subdivided into subspaces {X f } and write the free energy of {X} space as where F {X f } is the free energy calculated in the mean field approximation but the integration over X in this formula is extended over the region X f . At the energy minimum, the gap b (ε) is the same in each region but the matrices U f , Eq. (63), can be completely different provided the regions {X f } are isolated from each other. An interaction between the regions will couple U f with one another and the free energy counted from the mean field energies will be proportional to the sum of differences between the matrices U f . Following this idea we assume that the free energy functionalF for the variation of U is given by the following expressionF where δu (X) dX = 0 .
However, in order to proceed with the derivation ofF [u] we should find first the proper form of the mean field energy F . The derivative of this energy has been written in Eq. (127) and we have to reconstruct the energy itself. However, writing the free energy in a form an integral over λ 2 of the expression is not helpful because we would have to derive the σ-model for all λ, assuming at the same time that the energy of the excitations are smaller that the gaps for all λ including the smallest values of λ. Unfortunately, writing the explicit form of the function whose derivative is given by Eq. (127) is not simple. In order to avoid complicated derivations, we simply replace the function F λ in Eq. (127) by another function F λ , such that where the propagator D λ is introduced as This replacement is not justified by a small parameter but we do not expect that using Eq. (133) instead of Eq.
(127) may be dangerous. First, the propagators D λ eff and D λ have the same asymptotics in the limit of large ω, vk much exceeding the energy scale Γ, because in this limit the polarisation loops can be neglected. Moreover, they have the same asymptotics in the opposite limit of ω, vk much smaller than Γ because the function Ω (ω), Eq. (72), can be estimated in this limit as Ω (ω) ∝ ω 2 /Γ and the propagator D λ eff , Eq. (71), does not depend on λ being equal according to Eq. (134) to D λ . Therefore, the function D λ can serve as interpolation formula for the propagator D λ eff . As the gapless fluctuations arise due to the degeneracy of the order parameter, the replacement of D λ eff by D λ can numerically change coefficients in front of the gradient terms in the σ-model but not the form of the σ-model iteself. The new approximate coefficients will have the same order of magnitudes as the exact ones.
First, we guess the form of the function F writing it in the form and the matrix P (X, X ′ ) is related to the solution Q (X, X ′ ) of the self-consistency equations (60, 67, 69) as In Eq. (135) the symbol D λ eff P λ means and trace Tr includes both the trace of matrices and integration over X. So, the free energy F λ is composed as a functional of the solution Q of the self-consistency equations and the propagator D λ eff (X − X ′ ) containing in the polarisation loops the Green functions G, Eq. (53), with the same Q. In order to check that F λ is really the correct free energy we calculate its derivative Next, we calculate the energy of a slowly in space and time varying matrix U for each configuration of these matrices. We represent the matrix Q in the form with the matrix U defined in Eq. (63). Alternatively, the free energy functional F [u], Eq. (139), could be obtained by decoupling the interaction term, Eq. (35), singling out slow pairs of the fields and decoupling them with the help of the Hubbard-Stratonovich transformation. It should be kept in mind though that following such a procedure one would have to use as the bare interaction the propagator λ 2 D λ eff (X − X ′ ) with the preliminarily found gap function b (ε). Then, looking for the saddle point and expanding in fluctuations one should not differentiate D λ eff (X − X ′ ) with respect to the gap function.
The first term in Eq. (139) does not give any contribution for Q (ε, X) from Eq. (140) because U (X) drops out. As concerns the second one, we write and expand the second term in δU (X) up to second order terms in space and time derivatives. We express F [u] in the form where and f (ε) = (|ε| + a (ε)) signε.
Then, we expand the integrand in ω and q up to the second order and finally obtain the free energy functionalF that describes the gapless excitations in the form of two decoupled non-linear 2 + 1 dimensional σ-models, with u (X) belonging to the SU (2N ) group. The velocities v L x and v L y are defined as and the superscript L in Eqs. (146-148) specifies the two different pairs of the opposite pieces of the Fermi surface. We see that the fluctuations of the order parameters corresponding to these pairs of the pieces are independent. The dimensionless constant J 1 is written using the rescaled dimensionless functions, cf. Eq. (80), as This equation is valid at any temperature for T < T 0 . As concerns the constant J 0 , we write it in the limit of low temperatures, The precise numerical values are not important for the present study because we are interested mainly in universal properties of the model. The contribution of the fluctuations to the partition function takes the form As the σ-model, Eq. (147), is effectively three dimensional at T = 0, all contributions originating from fluctuations are convergent at T = 0. Moreover, they are small provided the inequality (75) is fulfilled. We can estimate the typical fluctuation of the unitary matrix u writing it as and expanding F L in δu up to quadratic terms. Then, we can easily calculate the average square of a small fluctuation in the Gaussian approximation and obtain using Eqs. (80, 93) that The free energy functional, Eq. (147) has been derived assuming that the relevant energies are lower than Γ, Eq.
(80), which thus is the upper limit of the integration over the energies. As a result, Eq. (154) shows that the average fluctuation is small at T = 0 as sin δ is the small parameter of our theory. At the same time, the large-N limit does not help because the fluctuations are not sensitive to N in this limit. This estimate demonstrates that the mean field theory is a good approximation in the limit (75). However, the situation changes at finite T because the contribution of the fluctuations diverges. We can see this from Eq. (153) replacing the integral over ω by the Matsubara sum T ω over the bosonic frequencies. Then, the integral over k in Eq. (153) diverges for ω = 0 in the infrared and more sophisticated methods are required for the investigation of the σ-model.
In order to finally confirm the applicability of the σ-model Eq. (147), we estimate the fluctuations of the eigenvalues of the matrix order parameter B (ε, X). For this purpose, we represent the order parameter as cf. Eq. (63). The h-excitations are gapped as seen in the corresponding free energy functional Herein, the constants K g , K 0 , and K 1 are of order unity. The estimate for the second moment of the fluctuations, [h 2 ] ii , is obtained similarly to Eqs. (153, 154), Even at finite temperatures, we do not need to struggle with the Matsubara sum as the estimate is given in a sufficiently accurate way by the integral in Eq. (156) as long as the temperature is not close to T 0 , where the constants K g , K 0 , and K 1 are small. Then, our estimate for the fluctuations of the eigenvalues of the matrix B is which is the same order of magnitude as in Eq. (154). So, the fluctuations of the modulus of the order parameter are small at all temperatures not close to T 0 in the limit (75). Note that, again, the large-N limit does not help. In order to complete the derivation of the non-linear σ-model, we now take into account the curvature term violating the symmetry between the superconducting and QDW states. The degeneracy between these states discussed so far is a consequence of the linearisation of the fermion spectrum. This degeneracy is broken as soon as the curvature of the Fermi surface is included. We shall see that due to the curvature effects the moduli of the order parameters become different, favouring superconductivity. This is because the curvature of the Fermi surface affects mainly the particle-hole channel and while the particle-particle channel is left unchanged. Formally, this difference splits the values of h. However, the most essential effect of curvature is the breaking of the SU (2N ) symmetry of the gapless excitations. This is because low energy excitations are very sensitive to violations of the symmetry. In the σ-model, the finite curvature of the Fermi surface gives rise to an additional term.
In order to derive this term, we write the second term in the free energy functional in Eq. (139) as with H curv introduced in Eq. (26) and expand it in H curv . The first order of the expansion does not contribute because the trace tr [U τ 3 ] vanishes as U is off-diagonal in Λ-space. The second order of the expansion yields an additional term F a in the σ-model, The integral over the components of the momentum parallel to the Fermi surface decouples from the integral over the perpendicular one. It formally diverges but we have to restrict the region of the integration by the length S closed by the gap at the hot spots. The rest is straightforwardly calculated and we come to the anisotropy term where the constant J 1 is given by Eq. (149). The full free energy functional of the σ-model takes the form The anisotropy term F a contains τ 3 , thus breaks the symmetry between the superconductivity and QDW and introduces a gap in the spectrum of the QDW excitations. However, this gap is small and fluctuations remain important. Their study is presented in the next Section.

IV. RENORMALISATION GROUP EQUATIONS AND SUPERCONDUCTING PHASE TRANSITION
As we have understood, the mean field approximation works well at T = 0 provided the inequality (75) is fulfilled. Due to the anisotropy term F a , Eq. (159), the system is at T = 0 in the superconducting state. However, fluctuations are strong at finite T and can eventually destroy the superconductivity. Therefore we consider the static limit of the σ-model neglecting time-dependent fluctuations and fluctuations of the modulus of the order parameter. This approximation is enforced by the inequality (75) but we believe that the σ-model can be used also for an arbitrary shape of the Fermi surface assuming that the coefficients of the σ-model are renormalised by those fluctuations.
It is convenient to introduce new coordinates R, and we consider here the σ-model only for one of the values L since they behave identically after rotation of the coordinates by π/2. The free energy functional F =F/T to be studied now takes the form In Eq. (163), S (a) [given by Eq. (122)] is the length of the Fermi surface under the gap for each hot spot near the critical line. We see that with increasing a, which corresponds going away from the critical line, the anisotropy parameter κ grows making superconductivity more favorable. The minimum of the energy is achieved for u anticommuting with τ 3 , {u, τ 3 } = 0. In considering the fluctuations near the minimum, it may be convenient to rotate matrix u by π/2, i.e. to replace u → τ 1 u. Then, the second term changes sign and the minimum is reached at u = 1. We assume that this rotation has been done. The standard method of studying 2D σ-model is to use the renormalisation group (RG) technique (see, e.g. Refs. 16,25 ) and, proceeding in this way, all imaginable types for the symmetries of the matrices u have been studied in the literature. However, for completeness, we present shortly the scheme of the derivation of the RG equations for the parameters of the model t and κ.
Following the RG method we represent the matrix u as where u 0 is a "fast" varying part andũ is a "slow" one. Integrating over the fast part, we obtain a renormalised functional. Due to the symmetries of the σ-model, the form of Eq. (53) is reproduced after the renormalisation but the constants t and κ are modified. This allows us to write RG equations for the dependence of t and κ on the "running" momentum cutoff. Proceeding in this way, we have to choose a suitable parametrisation for the unitary matrix u, The explicit form of w is where s = s + is a N × N hermitian matrix and p = p t is complex symmetric. Substituting Eq. (164) into Eq. (162) and expanding u in w, we write the functional F as where Φ = ∇ũũ + andF is the slow unrenormalised σ-model Integrating over fast variables one obtains where The integral in Eq. (170) gives contributions toF, thus renormalising the latter. The calculation of the averages in Eq. (170) reduces to Gaussian integrals. With the help of Wick's theorem, all the integrals can be reduced to the averages of the type Evaluating the averages, we obtain renormalised valuest andκ, which leads us to the RG equations with ξ = − ln k. At small k the RG flow for t (ξ) and κ (ξ) should stop at k ∼ κ, and these parameters have to be equal to the bare parameters t 0 , κ 0 in Eq. (162) at ξ = ln Λ c , where Λ −1 c is the lowest length in the theory. Then, the solutions for the effective temperature t (κ) and the anisotropy parameter κ 2 read where t 0 and κ 0 are the bare values of the parameters t and κ in Eq. (162). We see from Eqs. (177) that in the absence of anisotropy between the superconductivity and QDW (κ 0 = 0) the logarithm diverges and the effective temperature t goes to infinity when approaching the pole. This means that the long range order in the pseudogap state does not exist and correlation functions decay exponentially. A finite κ can stop the growth of t in the region when this value is not very large, t ≪ N . If this happens, which corresponds to large κ, the superconductivity is stabilised (though with a power law decay of correlation functions of the order parameter). Decreasing κ 0 or increasing t 0 , one approaches the region when t becomes of the order N and we can expect here a transition into the disordered pseudogap phase. The critical temperature T c of this transition has a completely different origin with respect to the one in the BCS theory. Its ratio to the value of the gap depends on many parameters and definitely is not universal.
Remarkably, the anisotropy parameter κ 2 vanishes simultaneously with the divergence of t. This means that at the transition not only the coherence of the superconducting state is lost but the system passes into the disordered mixture of the superconducting and QDW states.
Eqs. (176) demonstrate again that the limit N → ∞ is not special and the divergences in the perturbation theory would arise anyway without being suppressed by the parameter 1/N . An interesting possibility of the superconductor-QDW transition could exist if one could change the sign of the anisotropy parameter κ 2 in Eq. (162). Certainly, it is not possible to the right of the antiferromagnet-normal metal transition (where a > 0). However, the pseudogap cannot sharply disappear after crossing the line a = 0 to negative values of a. We thus do not exclude the possibility of changing the sign of κ 2 in the region of negative, not very large a. If this happened, the phase diagram in Fig. 3(c) would contain an interesting new phase of the QDW state with a transition into the disordered pseudogap state at a temperature T QDW and a quantum phase transition into the superconducting state at κ = 0.

V. SMALL PARAMETER OF THE THEORY
We have performed the calculations exercising care in keeping under control all the approximations. Let us summarize what we have assumed and what we should assume in addition in order to justify the results obtained.
First, we have assumed that the interaction is weak, Apparently, the very existence of the spin-fermion model can be justified only in this limit. Integrating over momenta, we assumed that the main contribution comes from the vicinity of the Fermi surface and this is only possible provided the inequality (178) is fulfilled. This is the usual weak coupling limit of a non-ideal Fermi gas.
The singular form of the propagator D (ω, q), Eq. (37), at the QCP makes the situation more complicated and the inequality (178) is not sufficient to control a perturbation theory. In order to make use of the diagrammatic technique, it was suggested 9,10 to introduce an artificially large number N of fermionic species. This became the basis of the subsequent considerations of the critical behavior near the QCP in the spin-fermion model with the conclusion that the earlier Hertz theory 8 properly described the system. Recently, it has been realized 12,13,14 that there is a class of diagrams that are, against the first intuition, not suppressed by 1/N in the large-N limit.
In our work, we confirm the conclusion of Refs. 12,13,14 that a large-N limit does not help in summing ladder-like diagrams that can interact with each other. This is analogous to the situation arising in the theory of the Anderson localisation where one is able to justify the use of diagrams without intersections leading to the conventional Drude transport formulas but corrections require the resummation of smaller diagrams leading to effective Goldstone modes called "diffusons" and "cooperons" 16 . We encounter such diverging contributions due to the effect of the spontaneous breaking of the symmetry, which generates gapless excitations. Estimating the contribution of the massive and massless modes we have seen that, indeed, the large-N limit does not help to keep the approximations under control. Fortunately, we have realised that a small parameter may exist in the theory, namely δ, where δ is the angle between the Fermi velocities of two hot spots connected by the vector Q.
Indeed, assuming allows us to distinguish between the scales of the momenta parallel and perpendicular to the Fermi surface when deriving the basic equations for the order parameter and neglect the massive modes when deriving the σ-model. Moreover, it justifies to neglect the quantum fluctuations at zero temperature and thus confirms the validity of the mean field approximation at T = 0. The inequality (179) corresponds to a tendency to nesting of the parts of the Fermi surface at the hot spots with either almost antiparallel, δ = π, or almost parallel, δ = 0, Fermi velocities v 1 and v 2 . Proceeding in this way, we used the random phase approximation to calculate the effective bosonic propagator D eff , Eqs. (39, 40). In doing so, we did not take into account logarithmic contributions to the Fermi velocities at the hot spots, however. These effects have been discussed in details in previous publications 9,10,14 but neglecting them has to be justified. All these contributions vanish in the limit N → ∞ resulting in an Eliashberg-type theory, i.e. a renormalisation of propagators while neglecting corrections to vertices. In principle, we could also assume the N → ∞ limit. However, as the latter does not help to justify the use of the σ-model, we would have to complement the N → ∞ limit with the inequality (179). Of course, this would be sufficient to prove the validity of our results but such a model may seem too artificial.
Actually, we are in a position to justify our results for any N including the physically relevant value of N = 1 if we assume instead of (179). The inequality (179) implies almost parallel Fermi velocities at the hot spots. The case of the almost antiparallel velocities does not help to keep the approximations under control. The inequality (180) is considerably stronger than the smallness in Eq. (179) but it allows to control the approximations employed during the derivation of the σ-model. What remains to be done is estimating the logarithmic contributions to the vertices and Fermi velocities as discussed in Refs. 9,10,14 . Explicit formulas are available in these publications and, for simplicity, we adopt them for our study. First of all, the logarithmic contributions calculated in Refs. 9,10,14 in the absence of the pseudogap diverge in the limit T → 0. As the pseudogap exists and is of the order Γ, Eq. (80), all the logarithmic contributions are cut off at low energies by Γ.
The renormalisation of the velocities v x and v y has been studied 9,10,14 introducing the ratio α α = v y /v x = cot(δ/2) An RG procedure yields the equation where ξ is a logarithmic running variable. Eq. (182) has a stable fixed point α = 0 and an unstable one at α = ∞. The general solution of Eq. (182) is where the constant C should be determined from the condition α (ln Λ) = α 0 .
Herein, α 0 is the bare value of α and Λ is the ultraviolet cutoff of the theory so that The flow determined by Eq. (183) should stop at ξ = ln Γ. The case of almost antiparallel Fermi velocities corresponds to small bare values α 0 ≪ 1 and one obtains Equation (185) shows that starting from an "incomplete" nesting with antiparallel Fermi velocities, α 0 ≃ (π − δ) /2 ≪ 1, one flows towards to a stronger one. However, one cannot move far if which holds unless Γ/Λ becomes exponentially small. The case of almost parallel Fermi velocities is characterised by large α 0 ≃ 2δ −1 ≫ 1 In this limit, the solution takes the form The flow moves away from the bare value but, again, cannot go far due to the low-energy cutoff Γ. So, for not exponentially small Γ one can neglect the renormalisation of α in both the cases of almost parallel and antiparallel Fermi velocities regardless of the value of N . The renormalisation of vertices can be expressed in terms of the parameter α. Using results of Refs. 9,10,14 , the vertex Γ φψ1ψ + 2 for the boson-fermion interaction in the lowest order in ln (Λ/Γ) is Of course, the correction is always small for N ≫ 1 but it remains small for an arbitrary N provided α ≫ 1. We see from Eq. (188) that the nesting with almost antiparallel Fermi velocities does not help to neglect the renormalisation of the vertex for N = 1 but the proximity to the nesting with the almost parallel velocities does this job very well. The same is true for other coupling constants of the theory. For example, the bosonic field strength renormalisation denoted in Ref. 14 by Z φ reads in the leading order in ln (Λ/Γ) In the limit α ≫ 1, we may expand arctan(1/α) in 1/α up to (1/α) 3 term, finding Again, we recover a small pre-factor of order (αN ) −1 in front of the logarithm. The correction to the inverse propagator D −1 (ω, q) has also been calculated in Ref. 14  The limit α ≫ 1 allows us to expand arctan(1/α) and we see that all the pre-factors in front of the logarithms are of the order (αN ) −1 , which justifies the neglect of these logarithms. (In addition, we should assume that the renormalised g is small, g ≪ 1, which is the case.) Thus, the condition (180) is really sufficient for proving the applicability of our theory for any number of flavours N including the actual physical value N = 1. We emphasize, though, that we do not expect a qualitatively different behaviour of the system for an arbitrary shape of the Fermi surface including scenarios in which the inequality (180) is not fulfilled, cf. Fig. 1(b). We expect that going from small angles δ to angles of arbitrary size is analogous to going from the weak to the strong coupling limit in the theory of (conventional) superconductivity. Although it is more difficult to control approximations in the latter limit, the superconducting gap is even larger.
The above inequalities allow one to derive the mean field equations (81-83). At the same time, the most interesting regime is where the QDW state and the superconductivity are close to each other and where the SU(2) fluctuations described by the σ-model are strong. This regime is possible when the anisotropy parameter κ in Eq. (162) is sufficiently small. As the shortest length in the σ-model, Eq. (162), is of order v/Γ, we assume that the inequality is fulfilled. Using Eqs. (80,93, 163) we come at a = 0 to the inequality If (192) is not fulfilled, fluctuations are small, the QDW state is supressed, and only superconductivity is possible. If (192) is fulfilled at a = 0, one can have the transition from the superconductivity to the pseudogap phase. However, going away from QCP by increasing a, one increases the anisotropy κ as the latter contains S (a), and for sufficiently large a the inequality (191) is no longer valid. In that region, one can have only the superconductivity.
We have seen that the major results correlate very well with properties of the high-T c cuprates. The present discussion leads to the conclusion that the ultimate quantitative theory of high-T c does not have to be a strong coupling theory (where the fixed point goes to strong coupling and thus escapes approximate treatment) and the great variety of physical effects can actually be obtained from the spin-fermion model in the limit of weak to moderate coupling.

VI. EXTENDED SPIN-FERMION MODEL AND PHYSICS OF CUPRATES
Until now all the calculations have been carried out for the spin-fermion (SF) model. Having in mind applications of our theory based on this model to the high-T c cuprates, we need to be able to translate physical quantities defined on the CuO 2 into corresponding ones defined in terms of the SF model and vice versa. Thus in this Section, let us discuss the connection between the one-band SF model considered throughout the paper and the original model of interacting electrons propagating in the CuO 2 lattice. This is particularly interesting because, as we are going to demonstrate, the particle-hole "valence bond order" of the SF model corresponds to a nontrivial charge modulation on the O atoms as well as an energy modulation on the Cu atoms within the CuO 2 lattice.
The one-band SF model written for the Cu lattice is partly phenomenological and thought of as a low energy limit of a more general model of interacting electrons moving on the CuO 2 lattice. At the same time, it is sufficiently general in the sense that the only information input about the system is the existence of an antiferromagnetic QCP. Therefore, in order to describe the physics on the CuO 2 lattice, we follow the same approach separating the degrees of freedom into paramagnon and fermion parts. We assume as previously that the formation of antiferromagnetic correlations takes place on the Cu atoms, whereas the fermions can be located both on the Cu and O atoms.
We introduce the Lagrangian L ext of such an extended spin-fermion model on the CuO 2 lattice, see Fig. 1(a), in a form that is a direct generalisation of Eqs. (1, 2) of the main text, In this Section, we restrict the number of fermionic flavours to its physical number N = 1. The first three terms in Eq. (193) describe the fermion motion while the last one is the Lagrangian of the critical antiferromagnetic field. As previously, L φ has the form L φ = 1 2 r φ r (τ ) −v −2 s ∂ 2 ∂τ 2 + a + 2 φ r (τ ) − 1 2 r,a 0 φ r (τ ) φ r+a 0 (τ ) + φ r (τ ) φ r−a 0 (τ ) , where r denotes a site of a Cu atom and a 0 a primitive vector in the Cu sublattice, i.e. it connects two neighbouring Cu sites in the directions either "upward" or "rightward". According to Fig. 1(a), there is an O(2p y ) orbital on the bond between the Cu sites r and r + a 0 if a 0 points in an upward direction and an O(2p x ) orbital if a 0 points to the right. As the field φ r (τ ) is critical, we assume that essential momenta q in its Fourier transform are close to the vector Q, which allows one to reduce Eq. (194) to Eq. (3) of the main text.
In the fermionic part of the Lagrangian L ext , the terms L Cu and L O are the Lagrangians for fermions on isolated Cu and O atoms, respectively, while L T introduces tunneling between them. Explicitely, L T = −t 0 r,a 0 ψ * r+a 0 ,σ (τ ) + ψ * r,σ (τ ) ψ r+ 1 2 a0,σ (τ ) + ψ * r+ 1 2 a0,σ (τ ) (ψ r+a 0 ,σ (τ ) + ψ r,σ (τ )) , where µ is the chemical potential and summation over repeated spin indices is implied. Within the SF model, we may assume that high energy degrees of freedom have been integrated out. While the energies ε Cu and ε O possibly differ from the energies of isolated Cu(3d x 2 −y 2 ) and O(2p x,y ) atomic orbitals, it is reasonable to assume that ε O − ε Cu > 0 and that this energy difference has become the highest energy in the low energy theory under consideration. Note that in the cuprates, ε Cu is close to the chemical potential µ.
where the wave functions extend over the region |R| a 0 . The function |Φ 3d x 2 −y 2 (R)| 2 is invariant under the exchange of the coordinates x ⇆ y, wheareas V O (r + R) changes sign. Thus, the integral in Eq. (208) vanishes and δE 3d x 2 −y 2 = 0. In contrast, the function |Φ 2px,y (R)| 2 is not symmetric under x ⇆ y and therefore δE 2px = −δE 2py is nonzero. As a result, the first-order correction to the transition energy ∆E x,y between the Cu(3d x 2 −y 2 ) and Cu(2p x/y ) states, δ∆E x,y = δE 3d x 2 −y 2 − δE 2px,y is nonzero as well, modulated with the same wave vectors Q ± as the QDW order, and having opposite signs for Cu(2p x ) and Cu(2p y ) states.