Gross-Neveu-Yukawa model at three loops and Ising critical behavior of Dirac systems

TL;DR

This paper advances the understanding of quantum critical behavior in Dirac systems by calculating three-loop critical exponents of the Gross-Neveu-Yukawa model in $4- ext{epsilon}$ dimensions, with applications to various condensed matter phase transitions.

Contribution

It provides the first three-loop calculations of critical exponents for the Gross-Neveu-Yukawa theory, extending previous two-loop results and applying them to real physical phase transitions.

Findings

Critical exponents computed at three loops agree with known two-loop results.

Emergent super-scaling relations hold up to three loops in single-component fermion limit.

Padé approximants of series provide insights into phase transition behaviors.

Abstract

Dirac and Weyl fermions appear as quasi-particle excitations in many different condensed-matter systems. They display various quantum transitions which represent unconventional universality classes related to the variants of the Gross-Neveu model. In this work we study the bosonized version of the standard Gross-Neveu model -- the Gross-Neveu-Yukawa theory -- at three-loop order, and compute critical exponents in $4-\epsilon$ dimensions for general number of fermion flavors. Our results fully encompass the previously known two-loop calculations, and agree with the known three-loop results in the purely bosonic limit of the theory. We also find the exponents to satisfy the emergent super-scaling relations in the limit of a single-component fermion, order by order up to three loops. Finally, we apply the computed series for the exponents and their Pad\'e approximants to several phase transitions of current interest: metal-insulator transitions of spin-1/2 and spinless fermions on the honeycomb lattice, emergent supersymmetric surface field theory in topological phases, as well as the disorder-induced quantum transition in Weyl semimetals. Comparison with the results of other analytical and numerical methods is discussed.

Tables

Table 1. Table 1: Critical exponents for the chiral Ising universality class to order ϵ 3 superscript italic-ϵ 3 \epsilon^{3} for three different choices of N 𝑁 N .

$N=2:$	${\nu }^{-1}$	$\approx$	$2-0.952ϵ+0.00723{ϵ}^2-0.0949{ϵ}^3$
	${\eta }_{\psi }$	$\approx$	$0.0714ϵ-0.00671{ϵ}^2-0.0243{ϵ}^3$
	${\eta }_{\varphi }$	$\approx$	$0.571ϵ+0.124{ϵ}^2-0.0278{ϵ}^3$
$N=1:$	${\nu }^{-1}$	$\approx$	$2-0.835ϵ-0.00571{ϵ}^2-0.0603{ϵ}^3$
	${\eta }_{\psi }$	$\approx$	$0.1ϵ+0.0102{ϵ}^2-0.033{ϵ}^3$
	${\eta }_{\varphi }$	$\approx$	$0.4ϵ+0.102{ϵ}^2-0.0632{ϵ}^3$
$N=1/4:$	${\nu }^{-1}$	$\approx$	$2-0.571ϵ-0.0204{ϵ}^2+0.024{ϵ}^3$
	${\eta }_{\psi }$	$\approx$	$0.143ϵ+0.0408{ϵ}^2-0.048{ϵ}^3$
	${\eta }_{\varphi }$	$\approx$	$0.143ϵ+0.0408{ϵ}^2-0.048{ϵ}^3$
$N=0:$	${\nu }^{-1}$	$\approx$	$2-0.333ϵ-0.117{ϵ}^2+0.125{ϵ}^3$
	${\eta }_{\psi }$	$\approx$	$0.167ϵ+0.0478{ϵ}^2-0.0469{ϵ}^3$
	${\eta }_{\varphi }$	$\approx$	$0.0185{ϵ}^2+0.0187{ϵ}^3$

Table 2. Table 2: Chiral Ising universality in D = 3 𝐷 3 D=3 : Inverse correlation length exponent 1 / ν 1 𝜈 1/\nu and anomalous dimensions η ϕ subscript 𝜂 italic-ϕ \eta_{\phi} and η ψ subscript 𝜂 𝜓 \eta_{\psi} for bosons and fermions, respectively. In this work , we provide results within the ( 4 − ϵ ) 4 italic-ϵ (4-\epsilon) expansion to order  ϵ 3 superscript italic-ϵ 3 \epsilon^{3} .

$N=2$		$1/\nu$	${\eta }_{\varphi }$	${\eta }_{\psi }$
this work (Padé $\left[2/1\right]$)		1.048	0.672	0.0740
$\left(2+ϵ\right)$, (${ϵ}^4$, Padé)Gracey:2016mio		0.931	0.745	0.082
functional RGPhysRevB.94.245102		0.994(2)	0.7765	0.0276
Monte CarloPhysRevD.88.021701		1.20(1)	0.62(1)	0.38(1)
$N=1$		$1/\nu$	${\eta }_{\varphi }$	${\eta }_{\psi }$
this work (Padé $\left[2/1\right]$)		1.166	0.463	0.102
functional RGPhysRevB.94.245102		1.075(4)	0.5506	0.0645
Monte Carlo1367-2630-17-8-085003		1.30	0.45(3)
$N=1/4$		$1/\nu$	${\eta }_{\varphi }$	${\eta }_{\psi }$
this work (Padé $\left[2/1\right]$)		1.419	0.162	0.162
functional RGHeilmann:2014iga		1.408	0.180	0.180
conformal bootstrapIliesiu:2015akf			0.164	0.164

Equations

$$\mathcal{L}=\bar{\psi }(/\partial +g\varphi )\psi +\frac{1}{2}\varphi (m^2-{\partial }_{\mu }^2)\varphi +\lambda {\varphi }^4,$$

$$\mathcal{L}=$$

$$+Z_{\varphi \bar{\psi }\psi }g{\mu }^{ϵ/2}\varphi \bar{\psi }\psi +Z_{{\varphi }^4}\lambda {\mu }^{ϵ}{\varphi }^4,$$

$$m^2$$

$$y$$

$${\beta }_y=$$

$$\displaystyle+\frac{y}{64}\Big{(}1152(7+5N)y^{2}\lambda+192(91-30N)y\lambda^{2}$$

$$\displaystyle+\big{(}912\zeta_{3}-697+2N(67+112N+432\zeta_{3})\big{)}y^{3}$$

$$\displaystyle-13824\lambda^{3}\Big{)}\,,$$

$${\beta }_{\lambda }=$$

$$\displaystyle-72Ny\lambda^{2}-816\lambda^{3}+\frac{1}{32}\Big{(}6912(145+96\zeta_{3})\lambda^{4}$$

$$+49536Ny{\lambda }^3-48N(72N-361-648{\zeta }_3)y^2{\lambda }^2$$

$$+2N(1736N-4395-1872{\zeta }_3)y^3\lambda$$

$$\displaystyle+N(5-628N-384\zeta_{3})y^{4}\Big{)}\,,$$

$$(y_{\ast },{\lambda }_{\ast }{)}_{\pm }=\left(\frac{1}{3+2N}ϵ,\frac{3-2N\pm s}{72(3+2N)}ϵ\right),$$

$${\gamma }_{\psi }=$$

$$\displaystyle+4N(47-12N))y^{2}+768y\lambda-2112\lambda^{2}\Big{)}\,,$$

$${\gamma }_{\varphi }=$$

$$\displaystyle+48\zeta_{3})+960Ny^{2}\lambda-2880Ny\lambda^{2}-6912\lambda^{3}\Big{)}\,.$$

$${\gamma }_{{\varphi }^2}=-2(6\lambda +N{y}^2-12Ny\lambda -72{\lambda }^2)-72{\lambda }^2(4Ny+87\lambda )$$

$$-4N{y}^3(4N-9+3{\zeta }_3)-\frac{3}{2}N{y}^2\lambda (11-24N+120{\zeta }_3).$$

$$\displaystyle\nu^{-1}=\theta_{1}=-\frac{d\beta_{\tilde{m}^{2}}}{d\tilde{m}^{2}}\Big{|}_{(y^{\ast},\lambda^{\ast})}=2-\eta_{\phi}+\eta_{\phi^{2}}\,,$$

$$\frac{1}{\nu }=$$

$$\displaystyle+\frac{\epsilon^{3}}{3888(3+2N)^{5}s^{3}}\Big{(}-227691(3+s)+4N(81(2170s-128871)+N(27(-2238507+554816s)$$

$$+2N(585(2414s-16143)+N(4N(5233698+N(1383001+16N(3832+54N-27s)-24986s)-371936s)$$

$$-3(8117973+761116s)))))+288(3+2N)s^2(2N(81+N(1917+4N(450+N(153+4N))))$$

$$\displaystyle-N(3+4N)(99+4N(21+N))s+81(3+s))\zeta_{3}\Big{)}\,.$$

$${\eta }_{\psi }=$$

$$\displaystyle-\frac{68607+2099304N+1629828N^{2}+1505352N^{3}+89536N^{4}+1248N^{5}}{7776(3+2N)^{5}s}-\frac{6(1+N)\zeta_{3}}{(3+2N)^{4}}\Big{)}\epsilon^{3}\,,$$

$${\eta }_{\varphi }=$$

$$\displaystyle+\frac{8829-24192N-1603476N^{2}-292200N^{3}-300224N^{4}-6112N^{5}+s(118926N^{2}+115564N^{3}+17312N^{4})}{1296(3+2N)^{5}s}\Big{)}\epsilon^{3}\,.$$

Full text

Introduction.

Dirac and Weyl fermions are an abundant form of quasi-particle excitations in condensed matter physicswehling2014 ; 2013arXiv1306.2272V appearing in very different materials ranging from graphene, via $d$-wave superconductors, to the surface states of topological insulators, or even three-dimensional (3D) materials such as Na3Bi and Cd3As2Liu864 ; PhysRevLett.113.027603 . While the physical origin of the quasi-relativistic energy dispersion can be quite different in various materials, it leads to universal low-energy properties shared by all these materials, such as, e.g., the density of states (DOS) and the concomitant thermodynamic properties or various response functions. Dirac and Weyl systems can also undergo transitions from their natural semi-metallic phase to a variety of broken symmetry phases as some parameter is varied. This includes continuous quantum transitions to interaction-induced, ordered, many-body ground statessorella1992 ; Khveshchenko2001 ; herbut2006 ; honerkamp2008 ; herbut2009 ; raghu2008 and a disorder-driven transition to a diffusive state with finite DOSPhysRevB.33.3263 ; PhysRevB.33.3257 ; PhysRevLett.107.196803 ; PhysRevLett.112.016402 ; PhysRevB.94.220201 ; PhysRevB.93.155113 . Depending on the broken symmetry of the ordered phase and the fermionic content, the corresponding transition represents a universality class with the concomitant critical behavior. Various of these transitions have been suggested to belong to the universality class defined by the chiral transition appearing in the 3D Gross-Neveu (GN) modelherbut2006 well-known in the context of high-energy physics and conformal field theoriesPhysRevD.10.3235 ; doi:10.1093/ptep/ptw120 . This applies, in particular, to the interaction-induced transition toward a charge density wave of electrons on the 2D honeycomb lattice that breaks the (Ising) sublattice symmetryherbut2006 , or the disorder-driven transition toward a diffusive metal in a 3D Weyl semi-metalPhysRevB.94.220201 .

Model.

The bosonized GN model – the Gross-Neveu-Yukawa (GNY) model – is represented by the Lagrangian

[TABLE]

which is defined in Euclidean space and is renormalizable in $D=4-\epsilon$ dimensions. It includes a real scalar field $\phi$ resulting from a Hubbard-Stratonovich decoupling of the four-Fermi interaction and lies in the same universality class as the GN model for (space-time) dimensions $2<D<4$HASENFRATZ199179 ; ZINNJUSTIN1991105 . Here, $\not{\partial}=\gamma_{\mu}\partial_{\mu}$ and we use a four-dimensional representation of the Clifford algebra $\{\gamma_{\mu},\gamma_{\nu}\}=2\delta_{\mu\nu}\mathbbm{1}_{4}$, with $\mu,\nu,=0,1,...D-1$. The conjugate of the Dirac field is given by $\bar{\psi}=\psi^{\dagger}\gamma_{0}$. The scalar field couples to the fermions with the Yukawa coupling $g$ and has a self-interaction (quartic term coupling) $\lambda$. For generality, we allow for an arbitrary number $N$ of four-component fermion species.

The precise determination of the critical exponents of the GN model is a formidable task and has been attempted by various methods, e.g., perturbativeRosenstein:1993zf ; Zerf:2016fti ; GRACEY1990403 ; Gracey:2016mio and non-perturbativePhysRevB.94.245102 ; PhysRevB.89.205403 renormalization group (RG) approaches, Monte Carlo (MC) simulationsKARKKAINEN1994781 ; PhysRevD.88.021701 ; 1367-2630-16-10-103008 ; 1367-2630-17-8-085003 ; PhysRevD.88.021701 ; PhysRevB.93.155157 and the conformal bootstrap, see Refs. Poland:2016chs ; Bashkirov:2013vya ; Iliesiu:2015akf . For the case of the purely bosonic $\phi^{4}$ theory with Ising symmetry in 3D, the development and comparison of these different methods has led to an impressive convergence across different theoretical approaches, settling beyond a three digit agreement, for example, for the correlation length exponent $\nu_{\mathrm{Ising}}\approx 0.630$0305-4470-31-40-006 ; PhysRevE.65.066127 ; PhysRevE.68.036125 ; PhysRevB.82.174433 ; Kos2016 . For the 3D Gross-Neveu, or “chiral Ising” universality class, the calculation of critical exponents has also been attempted by different methods; here, however, the situation is less settled, cf. Tab. 2. Recent progress in MC simulationsPhysRevD.88.021701 ; 1367-2630-16-10-103008 ; 1367-2630-17-8-085003 ; PhysRevD.88.021701 and the application of field-theoretical methodsGracey:2016mio ; PhysRevB.94.245102 ; PhysRevB.89.205403 could not resolve the existing discrepancies, and the difference between the results still shows up in the first digits. In fact, the paradigmatic role of the GNY model notwithstanding, its critical exponents are currently only known to two-loop order in the epsilon expansion near four dimensionsRosenstein:1993zf .

We have calculated the beta-functions and the critical exponents for the GNY model with the general number of fermion flavors $N$ at three-loop order, providing an important further step towards a more quantitative understanding of the fermionic universality classes. The full analytical expressions for general $N$ are presented below. Numerical values for the physically relevant cases $N=2,1,1/4$ and [math] are presented in Tab. 1, where we display the first three terms in the expansion in $\epsilon$ for the correlation length exponent and the anomalous dimensions. Furthermore, we use Padé approximants to extract estimates for the critical exponents at $\epsilon=1$.

The rest of the paper is organized as follows. First, we specify the RG procedure, the employed computer algebraical tools, and present the full set of three-loop RG functions. We then determine the fixed point solutions to order $\epsilon^{3}$ and calculate the anomalous dimensions and the correlation length exponent to that order. For the specific cases of physical interest $N=2,1,1/4,0$ we also calculate the Padé approximants for the universal critical exponents and discuss the physical applications for quantum phase transitions in Dirac and Weyl semimetals. Finally, we draw our conclusions.

RG procedure and tools.

We first explain how we set up the three-loop RG analysis in $D=4-\epsilon$ space-time dimensions. The bare Lagrangian is defined by Eq. (1) upon replacing the fields and couplings by their bare counterparts $\psi\to\psi_{0},\phi\to\phi_{0},g\to g_{0},\lambda\to\lambda_{0}$. The renormalized Lagrangian is then introduced as

[TABLE]

where $\mu$ defines the energy scale which parametrizes the RG flow of the couplings. We have defined the wave function renormalization constants $Z_{\psi}$ and $Z_{\phi}$ which relate the bare and the renormalized Lagranigan upon the field recaling $\psi_{0}=\sqrt{Z_{\psi}}\psi$ and $\phi_{0}=\sqrt{Z_{\phi}}\phi$. The explicit $\mu$ dependencies in $\mathcal{L}$ reflect that after introducing the integration over $D=4-\epsilon$ dimensional spacetime we shift $g^{2}\,\to g^{2}\mu^{\epsilon}$ and $\lambda\to\lambda\,\mu^{\epsilon}$. To simplify the notation in the following we introduce $y=g^{2}$. The renormalization constants for the mass term, the Yukawa coupling and the quartic coupling are further introduced by

[TABLE]

These relations provide the RG scale dependence of the renormalized quantities by taking into account the RG invariance of the bare quantities. We calculate the renormalization group constants $Z_{x}$ where $x\in\{\psi,\phi,\phi^{2},\phi\bar{\psi}\psi,\phi^{4}\}$ up to three-loop order employing dimensional regularization and the modified minimal subtraction scheme ($\overline{\text{MS}}$). To that end, we use sophisticated computer algebra which was established for higher-loop calculations in high-energy physics, in particular in the context of Standard Model of Particle Physics calculations: The complete set of Feynman diagrams is generated with the program QGRAFNOGUEIRA1993279 and further processed with the programs q2e and expHarlander1998125 ; Seidensticker:1999bb . Traces over matrix structures coming from the Clifford algebra and the tensor reduction of Feynman intergrals is then achieved with FORMVermaseren:2000nd ; Kuipers20131453 . The calculation of Feynman integrals is performed after reduction to master integrals via integration-by-parts identities. We evaluate the beta functions and the anomalous dimensions using two independent setups. In one of them we computed the vertex functions setting one or two external momenta to zero. In this case, the loop integrals are mapped to massless two-point functions that are implemented up to three loops in the code MINCERLarin:1991fz . In the second setup we introduce an infrared regulator for all the propagators as described in Ref. Chetyrkin:1997fm . Here, the loop integrals can be reduced to three-loop tadpole integrals that can be processed with the help of MATADSTEINHAUSER2001335 .

Beta functions.

The beta functions for the squared Yukawa coupling $y$ and the quartic scalar coupling $\lambda$ are defined as $\beta_{y}=dy/d\ln\mu$ and $\beta_{\lambda}=d\lambda/d\ln\mu$. The relation to the renormalization constants is derived from Eqs. (3)-(4). We work with rescaled couplings $y/(8\pi^{2})\to y$ and $\lambda/(8\pi^{2})\to\lambda$. The beta functions for the Yukawa and the quartic scalar coupling at three-loop order read

[TABLE]

where $\zeta_{z}$ is the Riemann zeta function. Our expressions fully agree up to two loops with the ones fromRosenstein:1993zf . Upon setting $y=0$, the beta function for the quartic coupling also agrees with the three-loop results for the real scalar $\phi^{4}$ theory with $\mathbbm{Z}_{2}$ or Ising symmetryKleinert:2001ax .

Fixed points.

The three-loop beta functions allow the determination of the RG fixed points of the system order by order in $\epsilon$ up to order $\epsilon^{3}$. Let us start the fixed-point analysis at the one-loop level, where the beta functions for $y$ and $\lambda$ give rise to four different fixed pointsPhysRevB.80.075432 : the unstable Gaußian fixed point with vanishing coordinates $(y_{\ast},\lambda_{\ast})_{0}=(0,0)$, the unstable bosonic Wilson-Fisher fixed point $(y_{\ast},\lambda_{\ast})_{\text{WF}}=(0,\epsilon/36)$, and a pair of fully non-Gaußian fixed points (NGFP)

[TABLE]

where $s=\sqrt{9+4N(33+N)}$. From the pair of NGFPs, the one with the negative solution is discarded as it has a negative quartic coupling. Here, we study the stable positive solution from Eq. (7), which we solve order by order in $\epsilon$. The expressions for general $N$ are lengthy at order $\epsilon^{3}$ and we refrain from fully displaying them here. Instead, we explicitly show the universal critical exponents at order $\epsilon^{3}$ as derived from the fixed-point solution.

Critical exponents.

The field renormalization constants are defined as the logarithmic derivatives of the wave function renormalizations of the fermion and the boson fields, $\gamma_{x}=d\ln Z_{x}/d\ln\mu$ for $x\in\{\psi,\phi\}$ and read

[TABLE]

Our expressions agree with the ones fromRosenstein:1993zf and Kleinert:2001ax in the corresponding limits. To obtain the fermion and boson anomalous dimensions characterizing the critical behavior, we evaluate these expressions at the NGFP and define $\eta_{\psi}=\gamma_{\psi}(y_{\ast},\lambda_{\ast}),\ \eta_{\phi}=\gamma_{\phi}(y_{\ast},\lambda_{\ast})$, see below.

Finally, we require yet another renormalization constant, i.e. the one related to the renormalized mass term. We define $\gamma_{\phi^{2}}=d\ln Z_{\phi^{2}}/d\ln\mu$. At three-loop order,

[TABLE]

The RG beta function of the dimensionless mass term $\tilde{m}^{2}=\mu^{-2}m^{2}$, then follows from Eq. (3) and reads $\beta_{\tilde{m}^{2}}=(-2+\gamma_{\phi}-\gamma_{\phi^{2}})\tilde{m}^{2}$. We extract the correlation length exponent at the stable NGFP $(y^{\ast},\lambda^{\ast})$ from the relation

[TABLE]

where $\eta_{\phi^{2}}=\gamma_{\phi^{2}}(y_{\ast},\lambda_{\ast})$. As our main result, we obtain at order $\epsilon^{3}$ for the inverse correlation length exponent and the anomalous dimensions

[TABLE]

Metal-insulator transition in graphene.

One of the motivations behind this study is the improvement of our understanding of metal-insulator transitions in graphene and related Dirac materials. The quantum phase transition from the semi-metallic state to the sublattice symmetry broken insulating state with charge order – the charge density wave (CDW) – belongs to the 3D Gross-Neveu(-Yukawa) universality classherbut2006 ; PhysRevB.89.205403 for fermion flavor number $N=2$, corresponding to an eight-component spinor $\psi$. These components reflect the presence of two sublattices of the underlying honeycomb lattice $A$ and $B$, two inequivalent Dirac points in the Brillouin zone at $K$ and $-K$ and two spin species $(\uparrow,\downarrow)$. The order parameter field is a spin singlet corresponding to a staggered density state with alternating charge densities on the different sublattices. Condensation to a non-vanishing vacuum expectation value $\langle\phi\rangle\neq 0$ spontaneously breaks the Lagrangian’s chiral symmetry and induces a finite gap in the energy dispersion. The quantum critical behavior of this transition for $N=2$ has previously been accessed by different approaches, most recently by higher-order perturbative and non-perturbative RG calculations and Majorana Quantum Monte Carlo simulations: For the purely fermionic Gross-Neveu model expanded in $D=2+\epsilon$ dimensions the RG functions are known to order $\epsilon^{4}$ and for $\epsilon=1$ yield an inverse correlation length exponent $1/\nu\approx 0.931$ after resummationGracey:2016mio , while the most sophisticated non-perturbative functional RG calculation gives a value of $1/\nu\approx 0.994(2)$PhysRevB.94.245102 . Novel lattice methods have managed to access the question avoiding the sign problem, but predicting a rather different value of $1/\nu\approx 1.20(1)$PhysRevD.88.021701 . We show a numerical evaluation of Eq. (11) for $N=2$ order $\epsilon^{3}$ in Tab. 1. We note that, e.g., for the inverse correlation length exponent, the prefactor of the order $\epsilon^{3}$ term is larger than the one for the order $\epsilon^{2}$ term. To obtain an estimate for the critical exponents at $\epsilon=1$, we first directly evaluate the series as given in Tab. 1 at $\epsilon=1$ which gives $\nu^{-1}(\epsilon=1)\approx 0.960,\eta_{\psi}(\epsilon=1)\approx 0.0404$ and $\eta_{\phi}(\epsilon=1)\approx 0.667$. Further, we employ the Padé approximant [2/1] to obtain for $\nu^{-1}_{[2/1]}\approx 1.048$ for $\epsilon=1$. The corresponding Padé approximants for the fermion and boson anomalous dimensions are $\eta_{\psi[2/1]}\approx 0.0740$ and $\eta_{\phi[2/1]}\approx 0.672$, respectively. Due to the considerable size of the order $\epsilon^{3}$ contribution the values from the direct substitution and the given Padé approximants are spread over a sizeable interval. We interpret this as a measure for the theoretical uncertainties.

Spinless fermions on the honeycomb lattice.

Another much studied case is the one of spinless fermions on the honeycomb lattice, which also undergo a metal-insulator transition for strong repulsive interactions. It corresponds to the universality class which is realized for $N=1$ in our model. Recent sign-problem free Monte Carlo studiesPhysRevD.88.021701 ; 1367-2630-16-10-103008 ; 1367-2630-17-8-085003 ; PhysRevD.88.021701 ; PhysRevB.93.155157 and functional RG approachesPhysRevB.89.205403 give improved estimates on $\nu^{-1},\eta_{\phi}$ and $\eta_{\psi}$, cf. Tab 2. Again, for our order $\epsilon^{3}$ estimates at $\epsilon=1$, we employ direct substitution of $\epsilon=1$ and obtain $\nu^{-1}(\epsilon=1)\approx 1.099,\eta_{\psi}(\epsilon=1)\approx 0.0773$ and $\eta_{\phi}(\epsilon=1)\approx 0.439$. The corresponding Padé approximants [2/1] for $\epsilon=1$ are $\nu^{-1}_{[2/1]}\approx 1.166,\eta_{\psi[2/1]}\approx 0.102$ and $\eta_{\phi[2/1]}\approx 0.463$.

Emergent SUSY in topological superconductors.

For $N=1/4$, the field content of the GNY model presented here is compatible with supersymmetry (SUSY) and an emergent SUSY scenario at the boundary of a topological phase that has been discussed inGrover:2013rc . In this case a superscaling relation $1/\nu=(D-\eta)/2$ with $\eta=\eta_{\psi}=\eta_{\phi}$ is satisfiedGies:2009az at the quantum critical point to all ordersHeilmann:2014iga . Estimates for the critical exponents have been calculated with SUSY functional RG methodsHeilmann:2014iga as well as with the conformal bootstrapBashkirov:2013vya ; Iliesiu:2015akf , cf. Tab. 2. Within our calculations, we find that the super-scaling relation is exactly satisfied order by order for $D=4-\epsilon$ up to three loops and for the Padé approximant [2/1] which gives for $\nu^{-1}_{[2/1]}\approx 1.419$ and for $\eta_{[2/1]}\approx 0.162$. For the Padé approximant [1/2], we observe a violation of the superscaling relation. We interpret this as a hint towards a superior behavior of the Padé approximant [2/1] over [1/2] for this model and order of the expansion, and then adopt this approximant also for the neighboring values of $N$, as listed in Tab. 2.

Disorder-driven transition in Weyl semimetals.

Another application concerns the replica limit $N\to 0$ which is believed to describe the transition from a relativistic semi-metallic state to a diffusive metallic phase in a 3D Weyl semi-metal. At one-loop order this unusual limit also gives rise to a NGFP which is non-trivial in both couplings $(y^{\ast},\lambda^{\ast})_{+}=(\epsilon/3,\epsilon/36)$. Beyond one-loop order, the limit $N\to 0$ suffers from the lack of multiplicative renormalizability of the model and contributions from evanescent operatorsPhysRevB.94.220201 ; Gracey:2016mio . A previous four-loop expansion of the purely fermionic Gross-Neveu model in $D=2+\epsilon$ has shown very large four loop coefficients in comparison with the three-loop termsGracey:2016mio . Here, we circumvent this problem by considering the GNY model in the same limit which allows us to provide three-loop estimates for the critical exponents of this disorder-induced transition. Our order $\epsilon^{3}$ estimates, with direct substitution of $\epsilon=1$ yields $\nu^{-1}(\epsilon=1)\approx 1.673$ and $\eta_{\psi}(\epsilon=1)\approx 0.167$. The Padé approximants are $\nu^{-1}_{[2/1]}\approx 1.610$ and $\eta_{\psi[2/1]}\approx 0.191$.

Conclusions.

We have studied the Gross-Neveu-Yukawa model at three-loop order in $D=4-\epsilon$ space-time dimensions and have extracted the solution of the stable non-Gaußian fixed point and the corresponding critical exponents to order $\epsilon^{3}$ for arbitrary number $N$ of fermion flavors. The model is believed to govern the universal critical behavior in a number of quantum phase transitions in interacting or disordered Dirac fermions that are of current interest, and we provided estimates for the correlation length exponent and anomalous dimensions based on Padé approximants for our three-loop series. Our calculation fully reproduces the previous two-loop results Rosenstein:1993zf at all $N$, and the three-loop Ising model results Kleinert:2001ax at $N=0$. For a single component fermion, coresponding to $N=1/4$ in our notation, our values of the exponents agree to all orders of the calculation with the scaling relation dictated by the supersymmetry that emerges in this limit. Finally, while the comparisons of the values of the critical exponents from different methods, i.e. Monte Carlo simulations, conformal bootstrap, perturbative and non-perturbative RG approaches, and our epsilon-expansion agree to a reasonable degree, they seem not to have fully converged yet.

In the future, it will also be interesting to study closely related field theories describing other patterns of symmetry breaking in presence of fermionsPhysRevB.80.075432 , e.g., the antiferromagnetic transition in the Hubbard model on honeycomb lattice, to compare the results with the large-scale Monte Carlo simulations PhysRevX.3.031010 ; PhysRevB.91.165108 ; PhysRevX.6.011029 .

Acknowledgements.

The authors are grateful to Holger Gies for discussions. M.M.S. is supported by DFG Grant No. SCHE 1855/1-1 and I. F. H. is supported by the NSERC of Canada.