# Critical behavior of the QED-Gross-Neveu-Yukawa model at four loops

## Abstract

We study the universal critical properties of the QED-Gross-Neveu-Yukawa model with flavors of four-component Dirac fermions coupled to a real scalar order parameter at four-loop order in the expansion below four dimensions. For , the model is conjectured to be infrared dual to the -symmetric noncompact P model, which describes the deconfined quantum critical point of the Néel-valence-bond-solid transition of spin-1/2 quantum antiferromagnets on the two-dimensional square lattice. For , the model describes a quantum phase transition between an algebraic spin liquid and a chiral spin liquid in the spin-1/2 kagomé antiferromagnet. For general we determine the order parameter anomalous dimension, the correlation length exponent, the stability critical exponent, as well as the scaling dimensions of singlet and adjoint fermion mass bilinears at the critical point. We use Padé approximants to obtain estimates of critical properties in 2+1 dimensions.

^{1}

## I Introduction

The study of critical phenomena represents one of the cornerstones of modern condensed matter physics, and the systematic understanding of such phenomena by renormalization group (RG) methods is widely acknowledged as one of the great triumphs of theoretical physics in the twentieth century Fisher (1974). The best-known examples of critical phenomena are phase transitions in systems with an -component vector order parameter, such as Ising (), XY (), or Heisenberg () magnets, which are typically described by the Wilson-Fisher fixed point of the bosonic vector model Wilson and Fisher (1972). Critical exponents for this model have been determined in successive refinements over the years, culminating in the recent tour-de-force calculation of critical exponents at six-loop order Batkovich et al. (2016); Kompaniets and Panzer (2017); Schnetz (2018). Combined with Padé or Borel resummation techniques, the expansion below the upper critical dimension of four is known to yield precise values for the critical exponents Guida and Zinn-Justin (1998).

While the vector model provides a satisfactory description of phase transitions obeying the Landau-Ginzburg-Wilson (LGW) paradigm, that is, transitions at which bosonic order parameter fluctuations are the only relevant long-wavelength degrees of freedom, much attention has been drawn in recent years to two classes of continuous quantum phase transitions at which the purely bosonic LGW approach fails. The first, fermionic quantum critical points, comprises phase transitions at which gapless fermionic degrees of freedom couple to order parameter fluctuations via a Yukawa-type coupling. In cases of interest this fermion-boson coupling is relevant at the purely bosonic (e.g., Wilson-Fisher) critical point and drives the system towards a new universality class with coexisting, and in many cases strongly coupled, bosonic and fermionic degrees of freedom. The prime example in this category is systems with Dirac fermion excitations at low energies, such as graphene or the surface of three-dimensional (3D) topological insulators, coupled with real Herbut (2006); Grover et al. (2014), complex Lee (2007); Roy et al. (2013); Nandkishore et al. (2013); Grover et al. (2014); Ponte and Lee (2014); Zhou et al. (2016); Li et al. (2017); Classen et al. (2017), or vectorial Sorella and Tosatti (1992); Herbut (2006); Honerkamp (2008); Sorella et al. (2012) order parameters. The corresponding critical points are described by the Gross-Neveu-Yukawa (GNY) model Zinn-Justin (1991), which can be studied by perturbative RG in dimensions, or its purely fermionic equivalent, the Gross-Neveu (GN) model Gross and Neveu (1974), which can be studied in dimensions. The critical exponents of the GNY model have been calculated at three-loop Mihaila et al. (2017) and four-loop Zerf et al. (2017) order recently, and those of the GN model have been determined at four-loop order Gracey et al. (2016). Interesting critical phenomena outside the reach of the purely bosonic vector model include the emergence of Grover et al. (2014); Sonoda (2011) and Lee (2007); Grover et al. (2014); Ponte and Lee (2014); Zerf et al. (2016); Fei et al. (2016) spacetime supersymmetry in the real (chiral Ising) and complex (chiral XY) GNY universality classes, respectively.

The second category of critical phenomena not captured by the standard universality classes are phase transitions involving dynamical gauge fields. In the condensed matter context these occur as a result of the fractionalization of microscopic degrees of freedom, under the influence of strong correlations, into slave particles with fractional quantum numbers. The paradigmatic example is the fractionalization of bosonic local moments into neutral fermionic spinons Wen (2004). One may further distinguish two subclasses of critical points in this category: those for which the gauge field deconfines only at the critical point itself, dubbed deconfined quantum critical points Senthil et al. (2004a, b), and those for which the gauge field deconfines in at least one of the two phases separated by the critical point. While the former subclass corresponds to (LGW-forbidden) transitions between conventional ordered phases, the latter describes transitions involving at least one fractionalized phase. When fermionic spinons acquire a Dirac dispersion Affleck and Marston (1988), one obtains a theory of Dirac fermions interacting with a gauge field as well as with a bosonic order parameter. The transitions of interest taking place in 2+1 dimensions, we will refer to such models as QED-GNY models, since the fermion-gauge sector is described by massless quantum electrodynamics (QED). Examples of transitions recently studied in this way include transitions between an algebraic spin liquid and a chiral spin liquid He et al. (2015); Janssen and He (2017) and between an algebraic spin liquid and a spin liquid Boyack et al. (2018), which are described by the chiral Ising (real order parameter) and chiral XY (complex order parameter) QED-GNY models, respectively. Interestingly, it was recently conjectured Wang et al. (2017) that a critical point in the first subclass, the deconfined quantum critical point between a Néel antiferromagnet and a valence bond solid on the 2D square lattice Senthil et al. (2004a, b), is dual to a critical point in the second subclass, that of the chiral Ising QED-GNY model with two flavors of two-component Dirac fermions Janssen and He (2017). By varying the number of flavors of Dirac fermions in the theory, one may obtain an infinite family of new universality classes, distinct from both the purely bosonic and GN/GNY universality classes.

Motivated by these recent developments, in this paper we study the critical properties of the chiral Ising QED-GNY model as a function of the number of flavors of four-component Dirac fermions. The study of the critical properties of this model in dimensions via the expansion was initiated in Ref. Janssen and He (2017), where calculations at leading (one-loop) order were performed; here we study this model up to four-loop order. The paper is structured as follows. In Sec. II we define the model. In Sec. III we discuss basic aspects of the RG procedure and give our results for the beta functions and anomalous dimensions. Results up to three-loop order are given in the main text; four-loop contributions are given separately in Appendix B and C. In Sec. IV and V we present our -expansion results up to for the usual thermodynamic critical exponents as well as the scaling dimensions of certain fermion bilinear operators; Padé approximants are then used to obtain rough estimates in . The procedure for the calculation of the stability critical exponent is briefly explained in Appendix E. In addition to our results for the chiral Ising QED-GNY model, we also compute the scaling dimension of fermion bilinears at the pure QED (see Appendix D) and GNY fixed points. In Sec. VI we discuss the application of our results to the duality mentioned above. In Sec. VII we discuss technical aspects of the automated procedure employed for the determination of renormalization constants at four-loop order. A brief conclusion is found in Sec. VIII.

## Ii Model

We study the chiral Ising QED-GNY model with flavors of four-component Dirac fermions , ,

(1) |

where the are Euclidean gamma matrices, is the field strength tensor, is a gauge-fixing parameter, and is a real scalar field. In the rest of the paper we will simply refer to this model as the QED-GNY model. The model has a global discrete chiral symmetry,

(2) |

where and , under which the fermion mass bilinear changes sign. The scalar mass squared tunes a quantum phase transition from a symmetric phase () with massive scalars and massless fermions, described for momenta by pure massless QED, to a phase with spontaneously broken symmetry () where the scalar acquires a vacuum expectation value and a fermion mass is dynamically generated.

As mentioned in Sec. I, when extrapolated to dimensions this model has been argued to be relevant to two problems of current interest in quantum magnetism. For , it has been suggested as a possible fermionic dual Wang et al. (2017) to the bosonic -symmetric noncompact P (NCP) model, which describes a deconfined quantum critical point between a Néel antiferromagnet and a valence bond solid on the 2D square lattice Senthil et al. (2004a, b). For , it describes a putative time-reversal breaking quantum phase transition between an algebraic spin liquid and a chiral spin liquid in a spin-1/2 kagomé antiferromagnet He et al. (2015); He and Chen (2015).

## Iii Renormalization group analysis

To perform an RG analysis of the model (II) we use the standard field-theoretic approach with dimensional regularization and modified minimal subtraction () Zinn-Justin (2002). Comparing the bare Lagrangian

(3) |

and the renormalized Lagrangian

(4) |

where is a renormalization scale, we find that the bare and renormalized fields are related by

(5) |

implying that the bare and (dimensionless) renormalized couplings are related by

(6) | ||||

(7) | ||||

(8) |

The quantum critical point is found by tuning the renormalized scalar mass squared to zero. To calculate the correlation length exponent however, one must determine the RG eigenvalue of the scalar mass squared at the critical point, for which we need the relation between bare and renormalized masses,

(9) |

Finally, one must also track the flow of the gauge-fixing parameter , which for the choice of gauge-fixing term in Eq. (II) is controlled by the relation

(10) |

The calculation of the renormalization constants , at four-loop order is done using an automated setup; the main technical steps of the computation are explained in Sec. VII.

### iii.1 Beta functions

The beta functions are defined as

(11) |

and we work with rescaled couplings for . Using Eq. (6) and the fact that the bare couplings are independent of , we have

(12) | ||||

(13) | ||||

(14) |

where we define the anomalous dimension

(15) |

associated to the renormalization constant . We express the four-loop beta functions as

(16) | ||||

(17) | ||||

(18) |

Here we display our result up to and including three-loop order; the four-loop contributions are lengthy and given in Appendix B and also in Ref. Sup (). The beta function for the gauge coupling is given by

(19) | ||||

(20) | ||||

(21) |

Likewise, the beta function for the Yukawa coupling is specified by

(22) | ||||

(23) | ||||

(24) |

where is the Riemann zeta function, with Apéry’s constant. Finally, the contributions to the beta function for the quartic scalar coupling are

(25) | ||||

(26) | ||||

(27) |

The beta functions (16)-(18) can be checked against known results in various limits. Setting and , the model reduces to the bosonic Ising universality class; our result for in that limit agrees with the known four-loop result vla (). Setting and , reproduces the four-loop QED beta function Gorishny et al. (1991). Setting only, our expressions for and agree with those for the pure GNY model in the chiral Ising class, which were computed at four-loop order recently Zerf et al. (2017). Finally, for the full QED-GNY theory with all three couplings nonzero we recover the one-loop beta functions recently obtained in Ref. Janssen and He (2017).

### iii.2 Anomalous dimensions

The anomalous dimensions of the fields , , , evaluated at the quantum critical point,

(28) | ||||

(29) | ||||

(30) |

are universal, gauge-invariant quantities. Considering first, from Eq. (12) we see that at a fixed point with nonzero gauge coupling one must have the exact relation Herbut (2007). This is a consequence of gauge invariance, since Eq. (12) follows from the fact that the gauge coupling and gauge field wave function renormalizations are related by a Ward identity. As shown at one-loop order in Ref. Janssen and He (2017), and confirmed at four-loop order in Sec. IV, the QED-GNY critical point indeed has , implying that exactly in dimensions. In this section we give expressions at four-loop order for and , which will then be evaluated at the quantum critical point in Sec. IV to yield the universal exponents and . As for the beta functions, we express the anomalous dimensions as a sum of contributions at fixed loop order,

(31) | ||||

(32) |

To calculate the anomalous dimensions, we make use of the chain rule when taking the derivative with respect to , as well as of the fact that the renormalization constants and have no dependence since the associated fields are gauge invariant Zinn-Justin (2002),

(33) |

for . As for the beta functions, here we only list the contributions up to three-loop order and provide the four-loop contributions in Appendix C and Ref. Sup (). The anomalous dimension of the scalar field is given by

(34) | ||||

(35) | ||||

(36) |

while for the scalar mass operator , we find

(37) | ||||

(38) | ||||

(39) |

### iii.3 Fermion bilinears

Besides and , another class of gauge-invariant operators one can consider are fermion bilinears. Restricting ourselves to Lorentz scalars, i.e., mass terms, a generic fermion bilinear can be expressed as a linear combination of an flavor singlet mass , which appears in the Yukawa interaction in Eq. (II), and an flavor adjoint mass

(40) |

where the generators of are linearly independent traceless Hermitian matrices. The scaling dimensions of the singlet and adjoint bilinears are in general different. To calculate the scaling dimension of a fermion bilinear where , we add it to the bare and renormalized Lagrangians,

(41) |

where we use the shorthand . This implies the relations

(42) |

and thus the beta functions

(43) |

where . Note that and are not separately gauge invariant, i.e., they depend on the gauge-fixing parameter , but all the gauge-dependent terms must cancel out in Eq. (43), since the fermion bilinears are gauge-invariant operators. Taking into account its gauge dependence, the fermion field anomalous dimension is given by

(44) |

where we have used Eq. (10) to express . We find that up to four-loop order depends on only in the one-loop term, as in pure QED Gracey (2007); Kißler and Kreimer (2017).

To calculate we compute the fermion two-point function at four-loop order with all possible single fermion bilinear and fermion bilinear counterterm insertions. Flavor-adjoint bilinear insertions preserve the chiral symmetry of the massless theory (for a proof of this statement, see Appendix A). In this case the scaling dimension of the bilinear is simply determined by the slope of the beta function (43) evaluated at the fixed point,

(45) |

where and are the anomalous dimensions of the fermion field and adjoint bilinear evaluated at the quantum critical point,

(46) | ||||

(47) |

Accounting for the dependence of , one has

(48) |

At four-loop order, we obtain

(49) |

with

(50) | ||||

(51) | ||||

(52) |

The four-loop contribution is given in Eq. (C.3) of Appendix C. The absence of gauge dependence up to four-loop order is an additional consistency check on the calculation. Furthermore we checked that our result for agrees with the results available in the literature in the pure QED limit Chetyrkin (1997); Vermaseren et al. (1997); Chetyrkin (2005); Czakon (2005), using the fact (discussed in Sec. V.1) that in the pure QED limit the singlet and adjoint bilinear scaling dimensions are identical.

By contrast with adjoint bilinear insertions, singlet bilinear insertions explicitly break the chiral symmetry and thus symmetry-breaking interactions will be radiatively induced. The only relevant (or marginal) such interaction below four dimensions is a interaction, which must be kept to preserve renormalizability of the theory. We must therefore additionally include the bare and renormalized couplings in the Lagrangian, with a new renormalization constant. This implies the additional relation

(53) |

and the corresponding beta function,

(54) |

with . Note that one has to introduce this extra coupling in order to obtain a finite/local result for starting at three-loop order. This is because radiative corrections to the cubic scalar vertex arise for the first time in the fermion two-point function with single mass insertions in three-loop diagrams (first diagram of Fig. 1). To calculate we use an equation analogous to Eq. (48), but the sum over must additionally include the couplings and .

In fact, as soon as a flavor-singlet fermion bilinear insertion is present, the theory is already nonrenormalizable at one-loop order without a cubic scalar vertex, because the scalar three-point function becomes divergent through fermion loop diagrams containing a single flavor-singlet bilinear insertion. Exactly these diagrams reappear as subdiagrams at three-loop order (first diagram in Fig. 1) and carry a subdivergence which renders the corresponding mass renormalization constant nonlocal (containing logarithms of ), if one does not subtract their subdivergence via a corresponding vertex counterterm insertion (second diagram in Fig. 1).

Moreover, even the scalar one-point function becomes divergent at the same loop order, when one allows for a operator insertion. So one has to introduce a tadpole counterterm into the Lagrangian in order to be able to render the one-point function finite. In our case we renormalize the parameter in a full subtraction scheme, meaning we have to all orders. This means effectively we do not need to consider any diagrams containing a -tadpole insertion, because for each such diagram there is a corresponding counterterm diagram which exactly cancels its contribution.

Since we are interested in the scaling dimension of in the massless, symmetric theory, it is sufficient to calculate and up to linear order in and . Note that will contain the singular ratio as a symmetry-breaking mass term can be radiatively generated at two-loop order by the cubic scalar vertex (i.e., the second diagram of Fig. 1 but with the counterterm insertion replaced by the cubic scalar vertex itself). Likewise, will contain the ratio as the cubic scalar vertex can be radiatively generated at one-loop order by a closed fermion loop with a single singlet mass insertion and three external scalar legs (i.e., the fermion loop subdiagram in the first diagram of Fig. 1). These singular ratios lead to terms linear in in and terms linear in in , i.e., mixing between the operators and . To linear order one thus obtains the linear system

(55) | ||||

(56) |

As in Sec. III.1 we find full cancellation of the gauge dependence in those beta functions Sup (), which is a strong check on the calculation. The scaling dimensions of the new eigenoperators are given by where

(57) |

are the eigenvalues of the matrix , evaluated at the quantum critical point. By inspecting the corresponding eigenvectors we find that can be associated with , and likewise .

## Iv Critical exponents

From the knowledge of the beta functions [Eq. (16)-(18)] and anomalous dimensions [Eq. (31)-(32)] one can calculate the usual critical exponents. We begin by searching for fixed points with couplings at one-loop order Janssen and He (2017). At that order one finds eight fixed points: the Gaussian fixed point , the conformal QED Di Pietro et al. (2016) fixed point , the Wilson-Fisher fixed point , a conformal QED Wilson-Fisher fixed point , two GNY-type fixed points with and , , and two fixed points with all three couplings nonzero. In agreement with Ref. Janssen and He (2017), of all those fixed points only one of the latter two is stable, the so-called QED-GNY fixed point:

(58) | ||||

(59) | ||||

(60) |

defining

(61) |

Furthermore, all three couplings are positive for all . In the following we study this fixed point at four-loop order, looking for a zero of the beta functions in the form

(62) |

with the one-loop coefficients given in Eq. (58).

Besides the previously determined exponent , the critical exponents we compute here are the scalar field anomalous dimension , the inverse correlation length exponent , and the stability critical exponent . The exponent is defined as the RG eigenvalue associated with the (relevant) scalar mass term,

(63) |

From Eq. (9) and Eq. (28), one obtains Zinn-Justin (2002)

(64) |

The exponent is defined as the RG eigenvalue associated with the least irrelevant operator in the basin of attraction of the fixed point (i.e., the critical hypersurface ), and controls the leading corrections to scaling. In practical terms, it is given by the smallest eigenvalue of the Jacobian (stability) matrix

(65) |

The approach utilized to diagonalize order by order in is tantamount to ordinary quantum-mechanical perturbation theory, and is briefly summarized in Appendix E.

At one-loop order, we find

(66) | ||||

(67) | ||||

(68) |

with and in agreement with Ref. Janssen and He (2017). At higher loop order, analytical expressions for the critical exponents with general can be obtained but are extremely cumbersome Sup (). As a nontrivial check on the calculation, we have verified that our four-loop result for , when expanded in inverse powers of to , agrees with the corresponding expansion result for the QED-Gross-Neveu (QED-GN) model in dimensions, when expanded to Gracey (1992a); Not (). Here we only give explicit expressions for the critical exponents for the case, relevant for the conjectured duality with the NCP model,

(69) | ||||

(70) | ||||

(71) |

and for case, appropriate for the spin-1/2 kagomé antiferromagnet:

(72) | ||||

(73) | ||||

(74) |

In both cases coefficients are given numerically to four significant digits.

### iv.1 Padé approximants

To obtain approximate values of the critical exponents in physical dimensions, corresponding to , we employ standard one-sided Padé approximants (see, e.g., Ref. Di Pietro and Stamou (2017)), defined as

(75) |

with , where is the desired loop order. They reproduce the -expansion results when expanded to and constitute an extrapolation from down to . We have also attempted to use two-sided approximants (see, e.g., Ref. Fei et al. (2016)) by combining information from the expansion of the QED-GNY model and the expansion of the fermionic QED-GN model,

(76) |

where are two-component spinors and the are gamma matrices. Just like the standard GNY fixed point in dimensions is believed to be in the same universality class as the fixed point of the purely fermionic GN model in dimensions Zinn-Justin (1991) when extrapolating , so the QED-GNY fixed point in dimensions is believed to extrapolate to the same universality class as one of the two charged critical points of Eq. (IV.1) Janssen and He (2017); Janssen (2016). Two-sided Padé extrapolation indeed appears to give accurate results for the pure GNY/GN fixed point Fei et al. (2016). However, since the gauge coupling in Eq. (IV.1) is strongly relevant near two dimensions, the charged critical point of interest, the QED-GN fixed point, is not perturbatively accessible at finite in a strict expansion, by contrast with the neutral fixed point of the pure GN theory Zinn-Justin (1991). Thus one is forced to proceed in a combined and expansion Janssen and He (2017); Janssen (